A Mathematical Journal
Vol. 7, No 2, (57 - 67). August 2005.

Quaternionic analysis and Maxwell’s
equations

Wolfgang Sprössig
Institute of Applied Analysis TU Bergakademie Freiberg

Prüfer-Strasse 9 09596 Freiberg Germany
sproessig@math.tu-freiberg.de

ABSTRACT

Methods of quaternionic analysis are used to obtain solutions of Maxwell’ s equa-
tions. By the help of time-discretisation Maxwell’s equations are reduced to an
equation of Yukawa type. Initial value and boundary value conditions are real-
ized by a representation formula in each time step. Approximation and stability
is proved.

RESUMEN

Se usan los métodos de análisis quaternionico para obtener soluciones de las ecua-
ciones de Maxwell. Con la ayuda de la discretización del tiempo, las ecuaciones
de Maxwell son reducidas a una ecuación del tipo Yukawa. Valores iniciales y
condiciones de valores en la frontera son realizadas por una fórmula de repre-
sentación en cada paso de tiempo. Se demuestra la aproximación y estabilidad.

Key words and phrases: Maxwell equations, quaternionic analysis, operator calculus.
Math. Subj. Class.: 35F10, 30G35



58 Wolfgang Sprössig
7, 2(2005)

1 Introduction

In 1873 J. C. Maxwell’ s fundamental paper A Treatise on Electricity and Magnetism
was published. Since this time generations of physicists and mathematicians felt
facsinated from the deepness and beauty of these equations. From the very beginning
scientists tried to give Maxwell’s equations a more simply algebraical structure maybe
in the form

Du + Au = F

with a suitable derivation operator D and a potential operator A. In this connec-
tion we should mention people like L. Silberstein (The theory of relativity, 1914), H.
Weyl (Raum -Zeit-Materie, 1921) and M. Mercier Expression des équations de elec-
tromagnetique au moyen des nombres de Clifford, 1946). New algebraical notions
were introduced and used for the description of Maxwell’s equations (for instance: C.
Lanczos (1929): spinors, A. Proca (1930): Clifford numbers, A. Einstein/A. Mayer
(1932): semi-vectors, F. Bolinder (1957): 4-d forms, G. Kron (1959): skew-symmetric
tensors and D. Hestenes (1968): multivector calculus).

We will use in our conception real and complex quaternionic analysis, more exactly
a quaternionic operator calculus, which also contains a corresponding quaternionic
numerical analysis.

We intent to apply a time-discretisation method (Rothe’ s method) in order to reduce
Maxwell’s equations to a disturbed Yukawa equation. The latter one is considered
under realization of initial and boundary values by means of a suited quaternionic
calculus. This paper belongs to a series of papers where we use Rothe’s method to
involve time-dependent problems in our quaternionic calculus. This paper can also
be seen as an alternative supplement to latest papers by V.V. Kravchenko et al (cf.
[3], [8] and [9].

2 Maxwell’ s equations in a chiral medium

Let G ⊂ IR3 a bounded domain with a sufficiently smooth boundary Γ. In MKS
Maxwell’s equations read as follows:

div D = ρ
ε0
, (1)

rot E = −∂tB, (2)
div B = 0, (3)
rot H = μ0J + ε0μ0∂tD, (4)

where μ0 is the permeability of the free space, ε0 the permittivity of the free space,
E the imposed electric field, B the magnetic field, ρ the (free) charge density, H
the effective magnetic field inside the dielectric medium, D the effective electric field



7, 2(2005)
Quaternionic analysis and Maxwell’s equations 59

inside the dielectric and J the charge density. In a homogeneous chiral medium G
the operators B and D have to fulfil the Drude-Born-Feodorov constitutive relations:

B : = μH + μβrot H, (5)
D : = εE + εβrot E, (6)

where β is the chirality measure of the medium. Further, we assume:
Initial value conditions:

E(x, 0) = E0(x) and H(x, 0) = H0(x). (7)

Boundary condition:

E(x,t) = g(x,t) (x ∈ Γ). (8)

Good references for Maxwell’s equations in chiral media are the book by A. Lakhtakia
[11] in 1994 as well as the article [1].

3 Preliminaries

Let us now introduce notations and operators needed in our approach. Let IH be the
set of real quaternions. Each quaternion permits the representation

a =
3∑

k=0

akek (ak ∈ IR; k = 0, 1, 2, 3),

where e0 = 1 and e1,e2,e3 are the so-called imaginary units. By definition these units
ek obey the following arithmetic rules:

e20 = 1,e1e2 = −e2e1 = e3,e2e3 = −e3e2 = e1 and e3e1 = −e1e3 = e2.

Addition and multiplication in IH turn it into a non-commutative number field. The
main-involution in IH is called quaternionic conjugation and defined by

e0 = e0, ek = −ek (k = 1, 2, 3),

which can be extented onto IH by IR-linearity. Therefore, we have

a = a0 −
3∑

k=1

akek.

Note that

aa = aa =
3∑

k=0

a2k =: |a|2IH.



60 Wolfgang Sprössig
7, 2(2005)

If a ∈ IH \{0} then the quaternion

a−1 :=
a

|a|2

is the inverse to a. For a,b ∈ IH we have ab = ba. The set of complex quaternions,
which we also need, is denoted by IH(C) and consist of all elements of the form

a =
3∑

k=0

akek (ak ∈ C; k = 0, 1, 2, 3).

By definition we state: iek = eki, k = 0, 1, 2, 3. Here i denotes the usual imaginary
unit in C. Elements of IH(C) can also be represented in the form

a = a1 + ia2 (ak ∈ IH; k = 1, 2).
Notice that the quaternionic conjugation acts only on the quaternionic units and not
on i.

Let now G ⊂ IR3 be a bounded domain with sufficient smooth boundary Γ. Assume
that all function spaces B(G, IH(C)) =: B(G) have the usual componentwise meaning.
Let u ∈ C1(G). The Dirac-Operator D is defined by

Du =
3∑

k=1

ek∂ku.

The operator D is right-linear with respect to complex numbers. On C2(G)) the
3-dimensional Laplacian permits the factorization

Δ = −D2.
We consider the disturbed Laplacian

Δ −ν2, (ν ∈ IR).
which is called Yukawa operator and acts on C2(G). This operator has the factorization
property

Δ −ν2 = (iν + D)(iν −D).
The factors iν + D and iν − D are called mutually generalized Dirac type operators.
Functions u ∈ ker(D + iν) are called left-(iν)-hyperholomorphic. The fundamental
solution of the Yukawa operator in IR3 is given by

Θν := −
1
4π

1
|x|e

−ν|x|.

Then the corresponding fundamental solution of the operator iν + D is given by

eiν (x) := (iν −D)Θν (x) =
(
iν +

x

|x|2 + ν
x

|x|
)

Θν (x)



7, 2(2005)
Quaternionic analysis and Maxwell’s equations 61

4 A time discretisation method

For simplicity we introduce the following abbreviations:

f :=
ρ

ε0ε
, a := −μ, b := −μβ, c := ε0μ0ε, d := ε0μ0εβ.

A simply calculation shows that Maxwell’s equations in a homogeneous chiral medium
transform into

DE = −f + ∂tH + b∂trot H, (9)
DH = μ0J + c∂tE + d∂trot E.

Let T > 0. The equations (9) are considered in the time-intervall [0,T ]. A decom-
position of [0,T )into n equal parts yields T = nτ, where τ is called the meshwidth of
the decomposition. We briefly write for k = 0, 1, ...,n :

Ek := E(kτ,x), Hk := H(kτ,x), fk := fk(kτ,x) and Jk := J(kτ,x).

We want to approximate the time derivatives ∂tE and ∂tH by the finite forward
differences:

Ek+1 −Ek
τ

and
Hk+1 −Hk

τ
.

More detailed we will consider the case β = 0. From (9) we obtain for k = 0, 1, ...,n−1:

DEk+1 = −fk +
a

τ
(Hk+1 −Hk), (10)

DHk+1 = μ0Jk +
c

τ
(Ek+1 −Ek). (11)

Setting now

ν2 = −ca
τ2

and L := −
√
−a
c
,

then we have

DEk+1 = −fk + Lν(Hk+1 −Hk), (12)
DHk+1 = μ0Jk −

1
L
ν(Ek+1 −Ek). (13)

Applying the Dirac operator D from the left we get

DDEk+1 = −Dfk + LνDHk+1 −LνDH,
= −Dfk −LνDHk + Lνμ0Jk −ν2Ek+1 + ν2Ek

and so

(Δ −ν2)Ek+1 = Dfk + LνDHk −Lνμ0Jk −ν2Ek. (14)



62 Wolfgang Sprössig
7, 2(2005)

With the same actions we obtain the dual relation

(Δ −ν2)Hk+1 = −μ0DJk −
1
L
νfk −

1
L
νDEk −ν2Hk. (15)

We intent to consider equations (14) and (15) in the hypercomplex setting in (cf.[4],
[12]). Main idea is to factorize the Yukawa operator on the left hand side. This leads
to

(D + iν)(D − iν)Ek+1 = −Dfk + Lνμ0Jk −LνDHk + ν2Ek. (16)

5 Borel-Pompeiu’s formula

Let G be a bounded domain in IR3 with the Liapunov boundary Γ and let n =
(n1,n2,n3) be the unit vector of the outward pointing normal at the point y ∈ Γ.
The kernel eiμ(x) function generates two important integrals: Teodorescu transforms,
which are defined by:

(T±iνu)(x) =
∫
G

e±iν (y −x)u(y)dy

as well as the Cauchy-Bizadse operators:

(F±iνu)(x) =
∫
Γ

e±iμ(x−y)n(y)u(y)dΓy.

These operators are well studied in several papers, see e.g. [6],[10] and [12]. In [7]
was obtained the following Borel-Pompeiu formula:

u = (D ± iν)T±iνu = T±iν (D ± iν)u + F±iνu in G, (17)
where u ∈ C1(G) ∪ C(G). Notice that Borel-Pompeiu’ s formula is also valid for
u ∈ W 12 (G). On the boundary Γ we have trΓu ∈ W 1/22 (Γ).

6 Representations

Applying Teodorescu transforms T±iν to formula (16) we get the iteration procedures:

Ek+1 = T−iνTiν
[
−Dfk + Lνμ0Jk −LνDHk + ν2Ek

]
+ T−iν Φ+ + Φ− (18)

and

Hk+1 = T−iνTiν

[
μ0DJk +

1
L
νμ0fk +

1
L
νDEk −ν2Hk

]
+ T−iν Ψ+ + Ψ−,

where Φ± and Ψ± belong to the kernel of the operators D±iν. Notice also that holds

Hk+1 = Hk +
1
L
νrot Ek.



7, 2(2005)
Quaternionic analysis and Maxwell’s equations 63

The unknown functions Φ+ and Φ− have to be determined now. In [12] is shown
that holds (D ± iμ)T±iμu = 0 and therefore also F±iμT±iμu = 0 for any function
u ∈ W 12 .This leads to

Φ− = F−iνEk+1 = F−iνtrΓEk+1 = F−iνgk+1

where gk := g(kτ,x). The determination of Φ+ is more complicated. Using Fiν Φ+ =
Φ+, which is a consequence of Borel-Pompeiu’s formula and trΓuk+1 = gk+1 we have

QΓ,−iνgk+1 = trΓT−iνTiν
[
−Dfk + Lνμ0Jk −LνDHk + ν2Ek

]
+ T−iνFiν Φ+,

where QΓ,−iν is one of the so-called Plemelj projections, which are defined by

n.− t.− lim
x′∈Ω±

x′→x∈Γ

F±iνu(x
′) =:

{
(PΓ,±iνu)(x), x′ ∈ Ω+ = G

(QΓ,±iνu)(x), x′ ∈ Ω− = IR3 \G
.

In [6] is shown that

trΓT−iνFiν : im PΓ,iν ∩W k−1/22 (Γ) → imQΓ,−iν ∩W
k+1/2
2 (Γ) (k > 1)

is an isomorphism. Notice that the pair of Plemelj projections act within correspond-
ing Hardy spaces (cf. [12],[6]). Further, we obtain

Φ+ = Fiν (trΓT−iνFiν )
−1 [trΓT−iνTiν (Dfk −Lνμ0Jk + LνDHk −ν2Ek)] +

+Fiν (trΓT−iνFiν )
−1QΓ,−iνgk+1.

Replacing Φ+ in (18) we get

Ek+1 = T−iν [Tiν
(
−Dfk + Lνμ0Jk −LνDHk + ν2Ek

)
+Fiν (trΓT−iνFiν )

−1trΓT−iνTiν (Dfk −Lνμ0Jk + νDHk −ν2Ek)
+Fiν (trΓT−iνFiν )

−1QΓ,−iνgk+1] + Fiνgk+1
= T−iν

[
I −Fiν (trΓT−iνFiν )−1trΓT−iν

]
Tiν (−Dfk + Lνμ0Jk −LνDHk + ν2Ek)

+T−iνFiν (trΓT−iνFiν )
−1trΓT−iν(D − iν)g̃k+1 + Fiνgk+1,

where g̃k+1 is a smooth extension of gk+1 into the domain G. We now introduce the
orthoprojections Qiν := I − IPiν , where IPiν := Fiν (trΓT−iνFiν )−1T−iν is a modified
Bergman projection onto the subspace ker(D + iν) and Qiν a projection onto the
subspace (D− iν)Ẇ 12 . Finally, we find

Ek+1 = T−iνQiνTiν
[
−Dfk + Lνμ0Jk −LνDHk + ν2Ek

]
(19)

+ T−iνP(D− iν)g̃k+1 + F−iνgk+1. (20)



64 Wolfgang Sprössig
7, 2(2005)

7 Realization of boundary conditions

The realization of boundary conditions for the imposed electric field makes use of the
follwing proposition:

Proposition 7.1 Let u ∈ L2(G, IH(C)). The condition

trΓT−iνu = 0

is sufficient and necessary that u ∈ imQiν ∩L2(G, IH(C)).
Proof. At first let u ∈ imQiν ∩L2(G, IH(C)). Then u permits the representation

u = (D − iν)w with w ∈ Ẇ 12 (G, IH(C))

and therefore trΓT−ivu = 0. On the other hand, if trΓT−iνu = 0 then immediately we
obtain from Hodge’s decomposition of the Hilbert space L2(G, IH(C))and the repre-
sentation of the generalized Bergman projection:

u = FΓ,iν (trΓT−iνFΓ,iν )
−1trΓT−iνu + Qiνu.

The first term vanishes and we have u = Qiνu. #

Let

Ẽk+1 := T−iνP(D− iν)g̃k+1 + F−iνgk+1,

Borel-Pompeiu’ s formula yields

(−Δ + ν2)Ẽk+1 = (D + iν)(D − iν)Ẽk+1 = 0.

Furthermore, we get from proposition (7.1)

trΓẼk+1 = trΓT−iνPiν (D − iν)g̃k+1 + PΓ,−iνgk+1
= trΓT−iν (D − iν)g̃k+1 − trΓT−iνQiν (D − iν)g̃k+1 + PΓ,−iνgk+1
= gk+1 −PΓ,−iνgk+1 + PΓ,−iνgk+1 − trΓT−iνQiν (D − iν)g̃k+1 = gk+1.

8 Approximation and Stability

From (15) we obtain by setting Fk := Ek + iLHk the representation

Fk+1 = TiνT−iν (iν−D)(fk−iμ0LJk)−iνTiν T−iν (iν −D)Fk + Tiνχ− + χ+,

where χ± belongs to the sets ker (D±iν). These functions can be defined by boundary
conditions. Further, we abbreviate Mk := fk−iμ0Jk. Using Borel-Pompeiu’ s formula

T−iν (iν −D)u = −T−iν (D − iν)u = FΓ,−iνu−u



7, 2(2005)
Quaternionic analysis and Maxwell’s equations 65

we get

Fk+1 = −T−iνMk + T−iνFΓ,−iνMk + iνTiνFk − iνT−iνFΓ,−iνFk + Tiνχ− + χ+.

Because the image of the Cauchy-Bizadse-type operators FΓ,±iν belongs to the kernels
ker (D ± iν) eventually we achieve the formula:

Fk+1 = iνT−iνFk −TiνMk + Tiνχ∗− + χ∗+.

The operator iνTiν is bounded in L2(G, IH(C)) which follows from ([2] Corollary 2.5).
There it is deduced the estimation for the generalized Teodorescu transform

‖T±iν‖[L2,L2] ≤
d

ν
,

where d depends on the diameter of the domain G and tends to zero for diam G →
0. Now it remains to analyze the approximation property of our semi-discretisation
procedure for convergence.

Put

L1 := c∂tE(t,x) −DH(t,x) and L2 := a∂tH(t,x) −DE(t,x).

Furthermore, we introduce the operators:

L1τ = c
E(t + τ,x) −E(t,x)

τ
+ DH(t + τ,x),

L2τ = a
H(t + τ,x) −H(t,x)

τ
+ DE(t + τ,x).

We have to estimate the differences Lj −Ljτ for j = 1, 2. With t = kτ it follows

|(L1 −L1τ )| ≤
c

τ
|(Ek+1 −Ek − τ∂tEk) + D(Hk+1 −Hk)|,

|(L2 −L2τ )| ≤
a

τ
|(Hk+1 −Hk − τ∂tHk) + DEk+1 −Ek|.

We intent to continue with the first estimate. It is easy to show that

|L1 −L1τ| ≤ cτ|∂ttE(kτ + θτ,x)| + μ0τ|∂tJk(kτ + θ′τ,x)| + cτ|∂ttE(kτ + θ′′τ,x)|
≤ τC1k (E,J)

with θ,θ′,θ′′ ∈ (kτ, (k + 1)τ). Analogously, we obtain for the second estimation

|L2 −L2τ| ≤ τC2k (H,f).

On this way the truncation error is estimated for sufficient smooth E, H, J and ρ



66 Wolfgang Sprössig
7, 2(2005)

Remark 8.1 With the same principle also the chiral case can be considered. One
obtains with similar calculations the following iteration procedure:

(
D + iν

1+βiν

)
Fk+1 =

Mk
1+βiν

+
(
D + iν

1+βiν

)(
βiν

1+βiν

)
Fk −

(
iν

1+βiν

)(
βiν

1+βiν

)
Fk.

Received: June 2004. Revised: November 2004.

References

[1] C. ATHANASIATIS, P.A. MARTIN AND I. G. STRATIS, Electromagnetic
scattering by homogeneous chiral obstacle: Boundary integral equations and
low-chirality approximations, SIAM J. Appl. Math. Vol. 59, No. 5, 1745 –
1762.

[2] H. BAHMANN, K. GÜRLEBECK, M. SHAPIRO AND W. SPRÖSSIG, On
a modified Teodorescu transform , Integral Transforms and Special Func-
tions, Vol. 12(2001), Number 3, 213-226.

[3] S. M. GRUDSKY, K. V. KHMELNYTSKAYA AND V. V. KRAVCHENKO,
On a quaternionic Maxwell equation for the time-dependent electromagnetic
field in a chiral medium J. Phys. A. Math.Gen. 37(2004), 4641–4647.

[4] K. GÜRLEBECK, Hypercomplex factorization of the Helmholtz equation,
ZAA, 5 (1986), 125-131.

[5] K. GÜRLEBECK AND W. SPRÖSSIG, Quaternionic Analysis and Bound-
ary Value Problems, Birkhäuser Verlag, Basel.(1990).

[6] K. GÜRLEBECK AND W. SPRÖSSIG, Quaternionic and Clifford Calculus
for Physicists and Engineers , John Wiley, Chichester, Sydney, New York.
(1998).

[7] K. GÜRLEBECK AND W. SPRÖSSIG, On a Teodorescu transform for a
class of metaharmonic functions. J. Nat. geopmetry, 21(2002), No.1-2, 17–
38.

[8] K. V. KHMELNYTSKAYA, V. V. KRAVCHENKO AND V. S. RABI-
NOVICH, Quaternionic Fundamental Solutions for Electromagnetic Scat-
tering Problems and Application, ZAA, Vol 22(2003), No. 1, 147–166.

[9] V. V. KRAVCHENKO, Applied Quaternionic Analysis, Research and Ex-
position in Mathematics Vol. 28(2003), Heldermann-Verlag,Lemgo.



7, 2(2005)
Quaternionic analysis and Maxwell’s equations 67

[10] V. KRAVCHENKO AND M. SHAPIRO, Integral representations for spatial
models of mathematical physics, Pitman research Notes in Mathematics 351,
Longman,, Harlow.(1996).

[11] A. LAKHTAKIA, Beltrami fields in chiral media, World Scientific, Singa-
pore. (1994).

[12] W. SPRÖSSIG, On the Decomposition of the Clifford Valued Hilbert Space
and their Applications to Boundary Value Problems, Advances in Applied
Clifford Algebras, 5(1995), 167 – 186.