Discretization of Maxwell’s Equations for Non-inertial Observers
Using Space-Time Algebra
Mariusz Klimek*, 1 , Stefan Kurz¶, 1 , Sebastian Schöps†, 1, 2 and Thomas Weiland‡, 2
1

arXiv:1611.07368v1 [] 21 Nov 2016

2

Graduate School of Computational Engineering, Technische Universität Darmstadt,
Dolivostraße 15, D-64293, Darmstadt, Germany
Institut für Theorie Elektromagnetischer Felder, Technische Universität Darmstadt,
Schloßgartenstraße 8, D-64289 Darmstadt, Germany
* Corresponding author, E-mail: klimek@gsc.tu-darmstadt.de,
Phone: +49 6151 16-24391, Fax: +49 6151 16-24404
¶ kurz@gsc.tu-darmstadt.de
† schoeps@gsc.tu-darmstadt.de
‡ thomas.weiland@temf.tu-darmstadt.de

Abstract: We employ Maxwell’s equations formulated in Space-Time Algebra to perform
discretization of moving geometries directly in space-time. All the derivations are carried out
without any non-relativistic assumptions, thus the application area of the scheme is not restricted to low velocities. The 4D mesh construction is based on a 3D mesh stemming from
a conventional 3D mesh generator. The movement of the system is encoded in the 4D mesh
geometry, enabling an easy extension of well-known 3D approaches to the space-time setting.
As a research example, we study a manifestation of Sagnac’s effect in a rotating ring resonator.
In case of constant rotation, the space-time approach enhances the efficiency of the scheme,
as the material matrices are constant for every time step, without abandoning the relativistic
framework.

1

Introduction

Numerical simulation of electromagnetic field in a rotating medium has been studied in frequencydomain [23, and references therein], as well as in time-domain [19, 18]. In all these approaches
it is assumed from the beginning that the angular velocity Ω is small, i.e., rΩ
c  1 with r the
Ω
distance from the centre of rotation, and in [23] additionally ω  1 with ω the frequency of the
electromagnetic field. These assumptions entail several drawbacks. First, they limit the range
of applicability, e.g., simulation of circular accelerators with respect to a rotating observer cannot be carried out. Second, to estimate the error committed in a low velocity approximation,
one certainly needs a fully relativistic description. Third, the physical interpretation is not
clear. Quantities as time and electromagnetic fields are assumed to be the same in a rotating
and a stationary reference frame. One may try to justify this assumption by the fact that the
difference due to relativistic effects is small, thus negligible. However, the quantities of interest,
such as the rotation-induced frequency shift, are also small. Therefore, one needs to distinguish
between quantities in different reference frames in order to interpret the computed quantities
properly.
Since Maxwell’s equations are invariant with respect to Lorentz transformations, it calls for
the relativistic treatment. The coordinate and metric free description of relativistic electromagnetism is outlined in [1]. In this paper we develop a fully relativistic computational scheme for
non-inertial observers. The scheme provides physical underpinning for further simplifications
and enables the possibility to account for the error of such simplifications. In principle one
could apply the scheme proposed in [20] to simulate the fields in a rotating reference frame.
However, this would require a 4D simplicial mesh generator, would not allow simple, explicit
time-marching and might require more computational resources, and one would also need to
1

employ the formalism of homology/cohomology groups as they are not trivial even in the case
of a relatively simple structure as a ring resonator considered in this paper.
We employ Clifford’s Geometric Algebra (GA) associated with Minkowski space-time; this
is known as Space-Time Algebra (STA). By imposing Minkowski metric, we restrict ourselves
to the setting of special relativity theory. However, we expect that all the steps can be repeated
for a general non-flat metric which corresponds to solving the electromagnetic field equations
in a background gravitational field. We approach the discretization of Maxwell’s equations
directly in space-time. The choice of GA as a mathematical language comes from the fact that
it enables concise and abstract notation, which in turn lets us focus on introduced mathematical
structures rather than their particular implementations.

2

Space-Time Algebra (STA)

The overview: In this section we introduce the basic concepts from relativity theory and GA,
which are necessary for the paper. We discuss the concept of physical and reference time, as
well as what do we mean by an observer. This concepts are crucial for the proper physical
interpretation of the quantities involved. We then introduce GA, STA and some concepts from
geometric calculus, which are used later.
To properly define the concept of the observer we need some basic notions, which are presented next.
The family of 1D curves filling densely, continuously and without intersections the space-time
region is called a congruence. If the tangent vector u to a curve (belonging to the congruence)
passing through any point is timelike, then the congruence is called timelike. Time-like congruence may be interpreted as follows. Each member of a family represents a world-line of a
particle. Particles can be imagined to form a body made of elastic rubber, which can deform
freely but cannot tear apart. As a consequence particles whose distance is small have similar
velocity.
The congruence is a purely geometrical concept; there was no notion of time so far. Since
each particle of the congruence can move with different velocity, the time passes at different
rates for each of them. The time measured by a particle is called the proper time τ of that
particle; we will also refer to it as the physical time. In simple words, if a particle looks at its
watch carried on its hand it sees τ .
The proper time τ is a natural parametrisation of any member of the time-like congruence.
Of course, in the definition of the natural parameter there is a norm (associated with the
metric) involved. In our case, the metric is not positive definite (which is sometimes called
pseudometric), thus resulting in the pseudonorm (not always greater than zero). However,
since we have restricted ourselves to time-like curves, the pseudonorm of the tangent vector
u is always positive, thus we can treat it as a usual norm, and thus the proper time τ being
a natural parameter is well defined (up to a constant as in the case of the positive definite
metric and the norm associated with it). The reference time tref is a parameter (not necessarily
a natural one) running along each member of the congruence. It is required that if we join
the points of different curves associated with the same value of tref the obtained hypersurface,
called hypersurface of simultaneity or constant time, is continuous.
With a time-like congruence there is associated the four-velocity field u(x), with x the
position in space-time, of tangent vectors. The congruence can be reconstructed from the fourvelocity field by means of integration. The integration constants are fixed by the particular
choice of time synchronisation.
The whole structure: time synchronisation and the four-velocity field u(x) (thus also the
time-like congruence) is referred to as the observer or the frame of reference.
The GA can be defined as a quotient algebra obtained by the factorisation of the full tensor
algebra with respect to a certain ideal, see [6]. However, we follow the approach presented in [5,
Chapter 4.1]. Namely, the resulting properties of the product in quotient algebra are treated
as axioms defining the geometric product.

2

The geometric product is
AB 6= BA ,

(noncommutative)

(AB)C = A(BC) ,

(associative)

A(B + C) = AB + AC ,

(left-distributive)

(A + B)C = AC + BC ,

(right-distributive)

where A, B and C are generic elements of GA called multivectors. If a is a vector, then its
square is a real number
a2 = g(a, a) ∈ R ,
(1)
with g the metric tensor.
The geometric product of two vectors can be split into its symmetric and antisymmetric
part
1
1
ab = (ab + ba) + (ab − ba) := a · b + a ∧ b ,
(2)
2
2
where · and ∧ denote the scalar and exterior product.
From the axioms follows that the scalar product is a real number
a·b=

1
(ab + ba) = g(a, b)
2

(3)

and coincides with the scalar product in linear algebra.
The exterior product a ∧ b of two vectors is called bivector or 2-vector. It can be interpreted
geometrically as an oriented surface of undefined shape and position, spanned by the vectors a
and b. Thus a circle in one point and a rectangle at the other point, but with the same area,
given by (a ∧ b)2 , and parallel to each other are associated with the same bivector.
Any multivector that can be written as an exterior product of n vectors
Vn = v1 ∧ · · · ∧ vn

(4)

is called an n-blade.
Any sum of n-blades
An =

X

Vni

(5)

i

is called an n-vector.
Let d be the dimension of the space. There are only two unit d-vectors, say, I1 and I2 ,
related via I2 = −I1 . By unit we mean that Ii2 = 1, with i = 1, 2. I1 and I2 correspond to
two possible orientations of the space. The chosen unit d-vector is denoted by I and called the
pseudoscalar. Once an ordered set of basis vectors γk , k = 1, . . . , d is given, the pseudoscalar is
uniquely defined via
γ1 ∧ · · · ∧ γd
I := r
.
(6)
2
(γ1 ∧ · · · ∧ γd )
Since we are interested in 4-dimensional space-time, we introduce more concepts particular
for this setting.
Basis vectors are denoted by γa , with a = 0, 1, 2, 3. The reciprocal basis vectors are denoted
by γ a and defined by requiring that they satisfy γa · γ b = δab , with δab the Kronecker delta.
The basis vectors of orthogonal Cartesian coordinate system t, x, y, z are denoted by γi , with
i = t, x, y, z. The reciprocal basis is γ t = γt and γ α = −γα , with α = x, y, z. This basis is
orthonormal by definition, thus we will exploit basic properties such as γi ∧ γj = γi γj = −γj γi
without explicitly commenting on that.
The pseudoscalar is given by
I=r

γ0 ∧ γ1 ∧ γ2 ∧ γ3

= γt γx γy γz .
2

(γ0 ∧ γ1 ∧ γ2 ∧ γ3 )

3

(7)

The 3D basis
σα := γ̂α ∧ γ̂0
(8)
p
with γ̂i := γi / |γi2 |, is used to embed 3D vectors in the space-time setting. This is also called
space-time split (first introduced in [8, Chapter 7], but the name appeared later in publications
by the same author) with respect to the observer with the four-velocity γ̂0 .
By an arrow over a symbol we denote 3D vectors. For example, the electric field
~ :=
E

3
X

E α σα .

(9)

α=1

We would like to justify the name “3D vectors”. From the point of view of 4D algebra (8)
are bivectors. However, one can build a 3D (not necessarily geometric) algebra treating σα as
basis vectors. As a consequence, one can identify σα with 3D basis vectors from traditional
vector analysis. Therefore, traditional 3D vector calculus can be used in a relative space of any
observer, provided one uses (8) and (9) to embed 3D vectors in 4D algebra.
The electromagnetic field in space-time is described by the Faraday bivector
~ + IB
~.
F := E

(10)

Knowing F one can split it into electric and magnetic field observed in the reference frame with
four-velocity u
~ = 1 (F − uF u)
E
2
~ = 1 (F + uF u) .
B
2I

(11)
(12)

This split is determined solely by the four-velocity field of the observer.
The geometric derivative
∇ :=

3
X

γ i ∂i = γt ∂t − γx ∂x − γy ∂y − γz ∂z

(13)

i=0

is a vector from the algebraic point of view, and an operator acting on the object to the right of
it from the operational perspective. In the same manner we define the 3D geometric derivative
~ :=
∇

3
X

σα ∂α .

(14)

α=1

For the comprehensive derivations and definitions of geometric calculus see [22] or [10], [5,
Chapter 6], [11]. As with the geometric product, we focus rather on a technical definition, i.e.,
how to evaluate integrals.
The geometric differential is defined as
dn x := e1 ∧ · · · ∧ en dx1 . . . dxn ,

(15)

where dxi are scalar differentials and ei are unit vectors spanning the hypersurface, thus e1 ∧
· · · ∧ en is a pseudoscalar of the hypersurface K n , compare [22, Equation (3.4)].
We define the integral of a multivector field A
Z
Z
dn xA := (e1 ∧ · · · ∧ en )Adx1 . . . dxn
(16)
Kn

Kn

as the limit of a sum as is usually done in traditional setting 1 .
1

In general, the integrand (e1 ∧· · ·∧en )A is not a scalar, thus one needs to define how to transport non-scalar
quantities (in order to sum them up), which in turn depends on the connection employed. In this paper we work
with flat Minkowski space-time, therefore this issue is irrelevant (the connection is trivial).

4

If the integrand is a scalar then it can be written as
(e1 ∧ · · · ∧ en ) · An ,

(17)

and the integral
Z

(dn x) · An

(18)

Kn

can be related to the integral of a differential form.

3

Maxwell’s Equations in STA

The overview: Here we start from traditional 3D form of Maxwell’s equations and constitutive
laws and derive their space-time equivalents expressed using GA. This is our starting point for
performing space-time discretization.
The traditional form of Maxwell’s equations separating space and time reads


~ ×E
~ = −∂t B
~
~ ×H
~ = ∂t D
~ + J~
∇
∇
(19)
~ B
~ =0
~ D
~ = %,
∇·
∇·
~ E,
~ H,
~ D
~ and B
~ are the electric charge density, the electric current density, the
where %, J,
electric and magnetic field strengths, and the electric and magnetic flux densities, respectively.
It should be stressed that Maxwell’s equations take the form (19) only in an inertial reference
frame.
Maxwell’s equations (19) do not form a complete set of equations. They need to be supplemented with constitutive material laws. Since we focus on lossless wave propagation, the
material laws in the reference frame, in which the materials are at rest, read
~ = εE
~
D

~ = νB
~,
H

and

(20)

with ε the electric permittivity and ν the magnetic reluctivity. We assume that the medium
is local and linear, possibly inhomogeneous and anisotropic. From mathematical point of view
this is equivalent to assuming that ε and ν are positive-definite, second rank tensor fields.
Equations (19) and (20) (with given boundary conditions and initial values) form a complete
set of equations, which we want to solve. Since we aim at performing space-time discretization,
we must express (19)–(20) in the space-time setting.
To this end, we combine all 3D scalar and vector fields into space-time counterparts in the
following way
~ + cI B
~
F := E

~
~ + 1 IH
G := D
c

~ 0,
J := (c% + J)γ

(21)

where γ0 is the four-velocity of the observer and c is the speed of light in vacuum.
Thus, Maxwell’s equations in the space-time form [5] ensue, viz.
∇ ∧ F = 0,

∇·G=J.

(22)

The space-time integral form of Maxwell’s equations is obtained in a manner described in
[9]. Namely, we integrate (22) over an arbitrary bounded 3D star domain K3 ⊂ K4 , with K4 the
4D space-time, and apply the fundamental theorem of geometric calculus [9, 22]. This results
in (SI units in square brackets)
I
I
Z
h
mi
2
2
.
(23)
(d x) · F = 0 [V · m] ,
(d x) ∧ G = (d3 x) ∧ J
C·
s
∂K3

∂K3

K3

The material mapping ξ is implicitly defined through
G = ξ(F ) .

5

(24)

If the materials are isotropic, then ε and ν are scalars, and we may easily obtain an explicit
formula for ξ, viz. (for a discussion on ξ for anisotropic materials, the reader is referred to [14])
~ + 1 IH
~ = εE
~ + ν IB
~ = ε 1 (F − γ0 F γ0 ) + ν 1 (F + γ0 F γ0 ) =
G=D
c
c
2
c2 2

i
ν
ν
1 h
ε + 2 F − ε − 2 γ0 F γ0 , (25)
=
2
c
c
i.e.,

i
ν
ν
1 h
ε + 2 X − ε − 2 γ0 Xγ0 .
(26)
2
c
c
q
In vacuum, we have that ε = ε0 , ν = ν0 and c = νε00 , from which the following simple mapping
ensues
G = ξ(F ) = ε0 F .
(27)
ξ(X) =

4

Discretization

The overview: We first present how we construct the 4D mesh complex from 3D primal mesh.
It is convenient due to the fact that there is a wide variety of available 3D mesh generators and
the material distribution is constant in time with respect to the 3D mesh. Then the integral
form of Maxwell’s equations is applied to this complex in order to discretize them. We then
introduce the interpolating functions; namely, Whitney ones. They are later used in creation
of material matrix, that is the discrete equivalent of the material mapping ξ.
In this section we use GA as well as linear algebra. Therefore, we introduce the following
notation. The elements of linear algebra (vectors and matrices) are denoted with a bar under
the symbol
(28)
f := [f1 , f2 , . . . , fn ] ,
where n follows from the context (for example, in the case of f it is the number of 2D facets).
Analogously M is the matrix with entries Mij . If the quantities from linear algebra are multiplied, the product from that algebra is employed.

4.1

4D Mesh Construction

We introduce the placement map, which specifies the observer. We derive the explicit expressions for the case of rotation with constant angular velocity. This map is used to create the 4D
primal mesh from a 3D one. In the end we introduce two possible choices for the barycentric
dual. We want to stress that, one can use the barycentric dual mesh (with broken edges) in
FIT as well. We present, however, the FIT derivations for a dual with straight edges, as that
is implemented in the code used to simulate numerical examples.
4.1.1

Placement Map

We extend the Lagrangian description of motion to our space-time setting by defining the
placement map
p : T × K3 → K4 , (tref , ~rref ) 7→ r ,
(29)
where tref is the reference time, ~rref is the reference position, and r is the physical position
vector in space-time 2 .
Our description may be thought of as a map rref 7→ r between reference and physical 4D
positions. Nonetheless, we explicitly separate the reference time tref
ptref (~rref ) := p(tref , ~rref ) ,

(30)

because this enhances later explanations.
When comparing our placement map with the traditional Lagrangian description, we observe
that
2 The reference domain is meshed using any 3D mesh generator; we obtain coordinates of the nodes in the
form [rx , ry , rz ], then ~
rref := rx σx + ry σy + rz σz .

6

1. We do not restrict the reference position ~rref to be the position at time zero. One reason
for this is that there is no unique notion of time (as discussed in Section 2).
2. We do not restrict the time component t of the position vector r to be the same as
the reference time tref . This allows to impose various synchronization of the clocks, i.e.,
change hypersurfaces of simultaneity.
4.1.2

Rotation with Constant Angular Velocity Ω

In order to deal with rotation we must define what we mean by constant angular velocity Ω.
(This discussion can be performed analogously to the discussion of constant acceleration in the
relativistic setting, see [16, Chapter 6.2]) If one applies the constantly accelerated motion as
v = at, then after a long enough time, one gets v > c, which is obviously not consistent with
relativity. Therefore, one rather takes as definition that the velocity at a time instant t + dt is
adt greater, with a relativistic formula for the addition of velocities, i.e.,
v=

v1 + v2
1 + v1 v2 /c2

(31)

employed. This corresponds to stating that the acceleration is constant as observed by the
accelerated observer. We do the same for constant rotation, that is, we define that the velocity
v at the radius r + dr is Ωdr higher than at r. (The reader familiar with the derivation for
accelerated observer will see that from mathematical point of view we just replace t with r and
a with Ω.)
Then at the radius r + ∆r, we have that
dv
∆r + O(∆r)2
dr
v(r) + Ω∆r
v 2 (r)
v(r + ∆r) =
= v(r) + Ω∆r − 2 Ω∆r + O(∆r)2 ,
2
1 + v(r)Ω∆r/c
c
v(r + ∆r) = v(r) +

(32)
(33)

that is, in the limit ∆r → 0, the velocity v satisfies the following ordinary differential equation


v2
dv
(34)
=Ω 1− 2 .
dr
c
By imposing the condition that the velocity at the center of rotation is zero v(0) = 0, we arrive
at the solution
 
rΩ
v(r) = c tanh
.
(35)
c
For small radii, (35) gives v(r) ≈ rΩ, which corresponds to the Galilean setting.
The cylindrical coordinates of ~rref are denoted by [r, ϕ, z]. For a system rotating with
constant angular velocity Ω around z-axis we have
 
rΩ ct
T
.
(36)
p = pt (~rref ) = [t, r cos(θ), r sin(θ), z][γt , γx , γy , γz ] , with θ = ϕ + tanh
c
r
Above, pt is such that it satisfies tref = t and ~rref = p0 (~rref ) ∧ γt coinciding thus with 3D
Lagrangian description of motion. Moreover, one can verify that the four velocity u of a
particle with reference coordinates [r, ϕ, z] is consistent with (35), see Appendix B.
4.1.3

Extrusion of the 3D Mesh

We depict the mesh construction in Fig. 1. Let M 3 be a 3D reference mesh depicted as a toroid
in the leftmost of Fig. 1.
We proceed to create images of M 3 under the placement map pt (M 3 ) at times t = t1 , t2 , ..., ti .
This is readily carried out by applying pt to the position vector ~rref of the nodes in M 3 while
preserving the connectivity, i.e., the incidence relations between elements. The resulting spacetime nodes are depicted as black points in Fig. 1 and the elements inherited from M 3 as black
horizontal lines joining these points. We then connect, along the “vertical” time line (see Fig. 1),
the nodes corresponding to the same reference point in M 3 . The resulting 4D mesh is referred
to as M 4 , and we refer to any n-dimensional element of M 4 as Kni .
7

pti (~rref )

~rref

Figure 1: Sketch of a space-time mesh used in simulation.
e i4
K
0

e i3
K
0
T3
T4

b

T2

T1
e i1
K
0

e i2
K
0

e i is denoted
Figure 2: A dual facet of M̄ 4 with n = 4 edges. The barycenter of the dual nodes K
0
by b. The facet consists of n = 4 flat triangles Ti .
4.1.4

Two Choices for Barycentric Dual

We define two kinds of dual meshes; namely, one is M̄ 4 that will be used in a 4D FIT context,
f4 in a 4D FEM context. In both cases the connectivity and dual nodes are
and the other is M
the same. As depicted in Fig. 1 with red squares, the dual nodes are the barycenters of the
4D cells of M 4 . Therefore, it appears natural to label a dual node within a 4D cell K4i with
the same index i. Nodes are used to form higher dimensional elements in the following way. If
K4i ∩ K4j = K3k then the dual nodes i and j are connected, and the resulting connecting edge is
labelled with index k. In the same frame of ideas, the dual edge k belongs to the dual facet l
if K2l ⊂ K3k , and so on for dual 3D and 4D elements.
Although we have defined the connectivity of all elements and positions of the dual nodes,
the exact geometry of the dual elements is not settled by these definitions. In fact, there are
two cases in which we are interested and describe them next.
The barycentric dual with straight edges M̄ 4 is obtained by joining the dual nodes with a
straight line segment. The dual facet with n edges (see Fig. 2) is composed of n flat triangles
with the nodes at dual edge’s nodes plus the barycenter of all dual nodes belonging to the dual
facet. If the dual nodes are coplanar our definition reduces to the flat dual facet. The dual 3D
element is defined as bounded by the dual facets.
f4 is defined as the 4D extension of the 3D construction described in
The barycentric dual M
[2, page 99]. Briefly speaking, the dual edge is defined as a broken line segment consisting of the
dual nodes and the barycenter of a corresponding 3D element (see Fig. 3). This construction
is extended to higher dimensional dual elements in such a way that the intersection of the
corresponding dual and primal objects is always the barycenter of the primal one. Moreover,
each dual object is a sum of flat simplices. This construction in 2D and 3D case is depicted in

8

K2j xbj
xij

ej
K
2

f4 (dotted line). The
Figure 3: The difference between dual edges of M̄ 4 (dashed line) and M
j
j
j
e in M̄ 4 is not the barycenter xb of K as in M
f4 . Nevertheless,
intersection point xij of K2 and K
j
2
2
b
we later use the interpolated fields at xj in our extension of FIT material matrices.
[2, Fig. 4.3 and Fig. 4.4], respectively.
Now, let us assume that the reference mesh M 3 bears some orientation of facets and edges.
To uniquely define the orientation of the 4D mesh M 4 we adopt the convention “first time then
space”. Therefore, if the edge is represented by a vector p, then the timelike facet associated
with that edge has the orientation given by the product u ∧ p.
Remark: In our space-time extension of FIT material matrices we apply several heuristics
such as the one explained in the caption of Fig. 3. In brief, one may perceive 4D FIT as
potentially more efficient, but there is no underlying theory, e.g., error analysis, i.e., convergence
proof. On the other hand 4D FEM is built on the sound mathematical framework being,
however, computationally more expensive. Therefore, we find it interesting to present both
approaches to the discretization of material mapping.

4.2

Maxwell’s Grid Equations

The overview: Here we discretize Maxwell’s equations. After general discretization in spacetime, we introduce some further simplifications coming from our particular construction of 4D
mesh. The physical interpretation of DoFs coming from these simplifications is given, and
relation to the discrete boundary conditions is briefly discussed. Our goal was to stick to the
3D FIT with leapfrog time integration, as close as possible.
We discretize Maxwell’s equations by applying their integral form (23) to the primal-dual
mesh pair. The (23)-left is applied to the primal grid and (23)-right to the dual one. This yields
I
I
Z
(37)
(d2 x) · F = 0 ,
(d2 x) ∧ G = (d3 x) ∧ J .
∂K3i

ei
∂K
3

ei
K
3

e i are sums
Please note that integrals over the boundary of primal and dual 3-cells ∂K3i and ∂ K
3
of integrals over 2D facets, i.e.,
I
I
XZ
XZ
(d2 x) · F =
(d2 x) · F ,
(d2 x) ∧ G =
(d2 x) ∧ G .
(38)
j

∂K3i

K2j

ei
∂K
3

j

ej
K
2

Therefore, in the case of a general mesh, like a simplicial 4D one as used in [20], we would
introduce scalar DoFs
Z
Z
Z
2
−1
2
−1
fj := (d x) · F ,
gj := I
(d x) ∧ G ,
ji := I
(d3 x) ∧ J ,
(39)
K2j

ej
K
2

ei
K
3

obtaining thus the discrete Maxwell’s equations
e 2g = j ,
D

D2 f = 0 ,
2

(40)

e are incidence matrices describing adjacency of 2D to 3D elements of the
where D2 and D
primal and dual mesh, respectively. Matrices Dn have entries

j
i

+1 if Kn ⊂ Kn+1 with matching orientation
n
j
i
Dij := −1 if Kn ⊂ Kn+1
.
(41)
with opposite orientation


j
i
0 if Kn 6⊂ Kn+1
9

p[tn ,tn+1 ]

l
x, y, z

n+1/2

m

K2j

el
ptn

t

bnm

x, y, z

x, y, z

Figure 4: Illustration of the relation between indices in (42). The leftmost is the reference
mesh. The l-th edge is depicted as a fat dot, and m-th 2D facet as a thick line. The rightmost
is the space-time mesh. 2D facets are depicted as lines.
Nevertheless, we have chosen a structured 4D mesh, and we shall exploit this fact in order to
simplify the discrete form of Maxwell’s equations. Since in this paper we are interested in wave
propagation phenomena, we assume J = 0 in (21) to simplify our derivations. This assumption
does not come from the limitations of the method.
Next, we rename and renumber the DoFs on the primal mesh as follows
(
Z
n+1/2
el
if K2j is timelike ,
fj = (d2 x) · F =:
(42)
bnm
if K2j is spacelike ,
K2j

where the relation between the indices j, l, m and n is as follows (see also Fig. 4). If the facet
K2j is timelike, it is the edge l in the reference mesh M 3 extruded in time via the placement
map pt such that t ∈ [tn , tn+1 ]. Similarly, if the facet K2j is spacelike, then it is the image of a
facet with index m in M 3 under the map pt with t = tn .
We perform a similar decomposition of the dual DoFs
(
Z
n+1/2
e j is spacelike
dl
if K
−1
2
2
gj = I
(d x) ∧ G =:
.
(43)
e j is timelike
hnm
if K
2
ej
K
2

In simple words, instead of an arbitrary enumeration of the 2D facets in M 4 , we rather label
them by using the facet/edge numbering in M 3 and the time step to which they correspond.
This numbering scheme is always possible and unique due to our particular construction of the
4D mesh. Such a decomposition aids the physical interpretation of the DoFs, as well as enables
a simple time marching procedure, which is in general not possible on unstructured 4D meshes
as considered in [20].
We now move towards the physical interpretation of the previously introduced DoFs. By
~ s and B
~ s we denote the 3D fields obtained by performing the space-time split by a stationary
E
~r
observer, i.e., the one with the four-velocity u = γt in (11) and (12). Analogously, the fields E
~
and Br are associated with (the four-velocity u given by (127) of) the rotating observer.
Spacelike case: Let us consider a spacelike facet K2j , and let us suppose that it is spanned by
orthonormal vectors γ1 and γ2 and parametrized with local coordinates x1 and x2 . The “normal” vector γ3 is uniquely defined by requiring that γi , with i = t, 1, 2, 3 form an orthonormal
set with the orientation γt γ1 γ2 γ3 = +I. Then
bnm

Z

2

Z

(d x) · F =

=
K2j


~s + IB
~ s) =
(γ1 dx1 ) ∧ (γ2 dx2 ) · (E

K2j

Z
=

dx1 dx2 (γ1 γ2 ) · (Esi γi γt + IBsi γi γt ) =

K2j

Z
K2j

10

dx1 dx2 Bs3 . (44)

Now we bring the expression above to the form from traditional 3D vector analysis. First, we
~ := ~ndS, where the normal vector ~n is defined as a cross
note that vector surface element dS
product of vectors spanning the surface; that is
~n := σ1 × σ2 = −Iσ1 ∧ σ2 = σ3 .

(45)

Next, the magnetic field component may be written as
~ s = ~n · B
~s .
Bs3 = σ3 · B

(46)

Therefore, in general it holds
Z
Z
Z
n
1
2 3
~
~ ·B
~s .
bm =
dx dx Bs =
dS(~n · Bs ) =
dS
K2j

K2j

(47)

K2j

Timelike case: Now, let us consider a timelike facet K2j , and let us suppose that it is spanned
by orthonormal vectors u and p. The local coordinates of the facet are: the proper time τ and
the distance along the edge l. Then we get analogously
Z
Z
Z
n+1/2
~r + IB
~ r ) = − dτ dl p~ · E
~r ,
el
= (d2 x) · F =
((u dτ ) ∧ (p dl)) · (E
(48)
K2j

K2j

K2j

where p~ := p ∧ u is a 3D vector representing an edge.
The same reasoning leads to the physical interpretation of the dual DoFs d and h. The
~ r and H
~ r , respectively,
DoFs denoted by e and h are related to the electric and magnetic fields E
in the moving reference frame, While the DoFs denoted by d and b are related to the electric
~ s and B
~ s , respectively, in the stationary reference frame.
and magnetic flux densities D
The physical interpretation of the DoFs enhances the imposition of boundary conditions.
For example, if we consider stationary Perfect Electric Conductor (PEC) boundary conditions,
then in 3D formulation one sets e = 0 on the boundary, which corresponds to setting the
~ s to zero, i.e.,
tangential component of E
~s = 0 ,
~ns × E

(49)

where ~ns is the vector normal to PEC in a stationary reference frame.
~ r to zero at the
If the PEC is rotating, then one needs to set the tangential component of E
boundary, i.e.,
~r = 0 ,
~nr × E
(50)
where ~nr is the vector normal to PEC in a rotating reference frame. Nevertheless, due to
the particular construction of the mesh and the resulting physical interpretation of the DoFs,
it is clear that (50) is equivalent to setting e = 0 on the boundary. Therefore, from the
implementation point of view nothing has changed due to movement.
Using (42) gives another advantage, namely, allows for a straightforward time marching
~s = 0
procedure. The Maxwell’s grid equations (40) (with omitted discrete equivalents of divB
~
and divDs = %s ) now read
bn+1 = bn + C en+1/2 ,

e hn ,
dn+1/2 = dn−1/2 + C

(51)

e [25] are, respectively,
which are first solved for n = 1, then for n = 2, and so on. C and C
primal and dual facets to edges incidence matrices of the reference 3D mesh, that is discrete
curl operators.
The signs in (51) depend on our convention “first time then space” for orientations of elements. The convention affects the decomposition of, e.g., D2 , in (40) to the temporal difference
and 3D curl-matrix C. For relation of our DoFs to the FIT ones, see (95).
We want to point out the similarity with 3D FIT + leapfrog [25]. It was shown, e.g., in
[14, 15], that leapfrog comes from the topology of the primal/dual mesh pair. However, this
11

is only true if the material matrices are local in time, that is they relate only the fields at the
same time step, i.e.,
dn+1/2 = Mεe en+1/2 ,

hn = Mνb bn .

(52)

Nevertheless, Maxwell’s grid equations, i.e., topological laws [7, Equation (4)], have not changed,
because the movement is encoded only in the geometry of the mesh, not in the topology.
However, our scheme does not always reduce to leapfrog as we consider space-time setting.
Material matrices may relate DoFs at different time steps, which in turn results in changing
the resulting time integration scheme.

4.3

Whitney Interpolating Functions

The overview: The concept of Whitney elements is explained in GA notation. Then we present
them for a particular reference element. The transformation to the physical domain is performed
later. The special properties of the lowest-order basis functions are discussed in the reference
element as well as in the physical ones, that is the images of the reference element with respect
to the particular diffeomorphism we use.
4.3.1

General Introduction

With every n-vector field An there is associated its discrete equivalent n-cochain An , which is
a vector in linear algebra sense, whose components are given by
Z
Ain := (dn x) · An ,
(53)
i
Kn

where Kni is the i-th n-dimensional cell of the mesh.
The deRham complex
∇∧

∇∧

∇∧

∇∧

A0 −−−−→ A1 −−−−→ A2 −−−−→ A3 −−−−→ A4
x
x
x
x
x





W 0
W 1
W 2
W 3
W 4
D0

D1

D2

(54)

D3

A0 −−−−→ A1 −−−−→ A2 −−−−→ A3 −−−−→ A4
will only commute for a proper choice of basis functions, namely, those belonging to the class
of Whitney elements. In general there are infinitely many members of this class. However, we
restrict our attention to the lowest-order elements, and thus they are uniquely defined up to
a sign coming from the orientation of the facets in Ξ. Therefore, in the sequel by Whitney
elements, we mean lowest-order Whitney elements. For lowest-order elements one can define
Whitney reconstruction operator (Whitney map) as
X
W n An =
Ain Nin (x) ,
(55)
i

with Nin (x) the n-vector valued basis function associated with n-cell Kni .
4.3.2

Our Choice of Whitney Forms on the Reference Element

So far, we have not specified any particular form of basis functions Nin (x). Due to the fact
that we intend to obtain a computational scheme, the explicit formulas for Whitney functions
Nin (x) are necessary. For the class of meshes discussed in this paper, 4D elements are tensor
products of 3D elements and 1D temporal intervals. Therefore, the extension of any particular
3D Whitney elements to 4D ones is straightforward. Since in our implementation we work with
hexahedral 3D mesh, our 4D elements are “4D hexahedrons”. Thus we shall extend here the
3D reasoning described in [13, Chapter 8.3.] to 4D.
We start by constructing the interpolating functions, that is Whitney elements, for a single
tesseract (“4D reference cube”) Ξ = [−1, +1] ⊗ [−1, +1] ⊗ [−1, +1] ⊗ [−1, +1]. The vertices of
12

(i, j)
Field component

(x, t)
Bx

(y, t)
By

(z, t)
Bz

(x, y)
Ez

(z, x)
Ey

(y, z)
Ex

2
Table 1: The physical interpretation of DoFs associated with the basis functions Ni±,
j± .

this tesseract are at the points (t, x, y, z) = (±1, ±1, ±1, ±1). Let us suppose we are interested
in the basis function associated with the facet [−1, +1] ⊗ [−1, +1] ⊗ −1 ⊗ −1. Since this is a txfacet, we know that the DoF living on such a facet is related to the (γt γx ) · F = E x component
~ + I B.
~ Clearly the function
of the Faraday field F = E
N1f =

1 x t
γ γ (1 − y)(1 − z)
16

vanishes at the other tx-facets and is orthogonal to the remaining facets, thus
Z
(d2 x) · N1f = δ1i ,

(56)

(57)

∆i

where ∆i , with i = 1, . . . , 24 are 2D facets of the reference tesseract.
One may construct a basis functions for all 2D facets, such that they fulfil the so-called
interpolation property
Z
(d2 x) · Njf = δji ,

i, j = 1, . . . , 24,

(58)

∆i

and give an exact interpolation for constant fields. Moreover, these basis functions guarantee
that the tangential component of the interpolated field F is continuous across a 2D facet, viz.
f · F1 = W
f · F2 ,
W

(59)

where F1 and F2 are the limits of the field F as approached from any two directions. Mathematically speaking
Fi = lim F (x + δvi ) , i = 1, 2,
(60)
δ→0

where vi is a unit vector describing the direction.
The space-time split of (59) leads to the more familiar 3D continuity relations
~ 1 · ~l = E
~ 2 · ~l ,
E

~ 1 · ~n = B
~ 2 · ~n
B

(61)

with ~n and ~l, respectively, a normal and a tangent vector to an arbitrary surface.
2
Regarding the physical interpretation of Ni±,
j± (defined in (141)) in the case of Φ1 being
an identity map, Table 1 provides a complete description.
We observe that the chosen basis functions Ni2 , besides being Whitney forms (see Appendix D), do also comply with the requirements of an energetic approach (as introduced in
[4]) in our space-time extension [3, equation iii)], that is
Z
fj ,
Nj2 |d4 x| = I −1 W
(62)
Ξ

fj is the bivector associated with the facet ∆
e j dual to
where Ξ is the reference tesseract and W
∆j , which is only the part contained in the reference tesseract Ξ. Please note that the dual
e j in the whole mesh complex are sums of the facets Φi (∆
e j ) in neighbouring 4D cells
facets K
2
j,i
i
e , being the part of the dual facet K
e j contained
K4 . Therefore, the bivectors associated with K
2
2
fj,i . This explanation is sketched schematically in
in the primal cell K4i , are later referred as W
Fig. 5.
The property (62) transfers to the physical domain as follows. Suppose it holds in the
reference domain. In this reasoning the quantities in the reference domain are denoted by a
bar.
13

Ξ
e l0
∆

el
∆

Φi

Φi 0
ej
K
2

e j,i
K
2

K4i

e j,i
K
2

0

0

K4i

e l0 and ∆
e l of the reference tesseract Ξ are mapped via Φi and Φi0 to
Figure 5: The facets ∆
j,i0
j,i
e
e j contained in K i and K i0 , respectively. Namely
e
K
theparts
4
4
2
 K2 and K2 of dual facets
0
0
e j,i is thus
e j,i ∪ K
el = K
e j,i . The bivector associated with K
ej = K
e l0 = K
e j,i and Φi0 ∆
Φi ∆
2

2

2

2

2

fj = W
fj,i + W
fj,i0 .
W
First, we note that
Z
Nj |d4 x| = lim
4

X

∆ x→0

K4i

Nj (xa )|∆4 xa | =

X

= lim
4

∆ x̄→0

∆ x̄→0



Φi N̄j (x̄a ) Φi |∆4 x̄a | =

(64)


Φi N̄j (x̄a )|∆4 x̄a | =

(65)

a, x̄a ∈Ξ

X

= lim
4

(63)

a, xa ∈K4i

a, x̄a ∈Ξ


= Φi  lim
4

∆ x̄→0



Z
N̄j (x̄a )|∆4 x̄a | = Φi  N̄j |d4 x̄| ,


X
a, x̄a ∈Ξ

(66)

Ξ

where |∆4 x| is a finite 4D volume element and we assumed that interchanging the order of
diffeomorphism and the limit of the sum is justified. Moreover, we have assumed that both
the reference and the physical domain are equipped with the flat Minkowski metric, thus the
connection is trivial and is not necessary in the sum.
So we calculate in physical domain


Z
Z


¯
f
=
Nj |d4 x| = Φi  N̄j |d4 x̄| = use (62) = Φi I¯−1 W
K4i

Ξ

 
 ¯
¯
f = I −1 Φi W
f
= Φi I¯−1 Φi W
, (67)
where Φ is the mapping to the physical domain, and we used the fact that the unit pseudoscalar
is invariant of the mapping (I was defined so). Therefore (62) will hold in physical domain if
 
¯
f
f,
Φ W
=W
(68)
which will hold if Φ maps
 the reference dual facet to the physical dual facet. Therefore, Φ ∈
¯
f
f ⇒ (62) . So if the mapping preserves (weighted) barycenters
affine mapping ⇒ Φ W = W
then (62) is implied, i.e., one cannot came up with a map that preserves barycenters and violates
(62).
14

K4i

Ξ

Φi

e j,i
K
2

el
∆

e l are mapped in general to Φp (∆
e l ), which does not necessarily coincide
Figure 6: The facets ∆
j
e
e l) =
with the part of the dual facet K2 . However, for the mapping discussed in this paper Φp (∆
j
e
K2 .
4.3.3

Basis Functions in the Physical Domain

The basis functions Nin satisfy the following set of properties:
1. they satisfy the interpolation property (58),
2. they fulfil (134), i.e., guarantee that (54) commutes (see Appendix D),
3. they are eligible for 4D energetic approach, i.e., meet (62).
All these properties hold for the reference tesseract; however, in many applications one is
interested in using meshes with non-cubic elements in order to approximate geometries more
accurately. Moreover, if one even uses a 3D Cartesian grid, the 4D elements are certainly
deformed by the movement, which is encoded in the geometry of the 4D mesh. Therefore, one
needs to construct analogous functions for a deformed tesseract. This is conveniently done by
first defining a mapping (diffeomorphism) from the reference to the physical tesseract and then
by transforming the basis functions accordingly.
We define a set of mappings
Φp : Ξ → K4p
(69)
from the reference domain to the physical 4D cells. One of important features of these mappings
is that they are affine maps, consequence of which is that the barycenters of any object in the
reference domain are mapped to the barycenters of a corresponding object in the physical
domain. For details see Appendix E.
Property 1. involves the integral of a differential form, which is invariant with respect to the
change of the metric properties; therefore, this property will hold for the transformed functions.
As explained in [12, text below Equation (2.13)], the transformed (pulled-back) functions will
satisfy the Property 2. on a deformed element. Property 3. is not guaranteed in general as the
dual elements of Ξ are not mapped to the barycentric dual elements of the deformed tesseract.
However, for the mapping we use, i.e., Φp , this is the case, that is barycenters are mapped to
barycenters, thus the reference barycentric dual is mapped to the physical barycentric dual; see
Fig. 6. Hence, all required properties are inherited by the deformed tesseract. Even if that was
not the case, one can still use the energetic approach, by simply stating that a different dual
(not the barycentric one) is used. Since in general, one does not need to calculate the dual
mesh in FEM, this bears no additional computational cost. The price to pay is that one cannot
~ and H
~ without calculating the dual.
easily relate d and h to D
In other words, properties 1. and 2. come from differential topology, thus change in the
metric, i.e., shape of an element, does not influence them. In contrast, the property 3. depends
on metric in a twofold way. First, the definition of the barycenter depends on the notion of
distance, thus the barycentric dual depends on metric. Second, the integrand in (62) is not a
scalar; therefore, the Levi-Civita connection is used, which depends on the metric.
There is no unique approach in the literature how to define the basis vectors in the physical
domain. We opt for the following definitions
γ̄i :=

∂x
,
∂ x̄i

γ̄ i := ∇x̄i (x) ,

15

(70)

where x is a position vector in physical space and x̄i (x) are reference coordinates (in γ̄i basis)
of x, i.e.,
x̄i = γ̄ i · Φ−1
(71)
p (x) .
Then we easily verify that they “transform” according to
∂xj ∂
∂xj
∂
x=
x=
γj ,
i
i
j
∂ x̄
∂ x̄i
∂ x̄ ∂x
∂
∂ x̄i j
γ̄ i = ∇x̄i = γ j j x̄i =
γ .
∂x
∂xj
γ̄i =

(72)
(73)

It means that basis vectors transform like vectors, and the reciprocal basis vectors as one-forms.
It is less surprising if one realises that γi , ∂i and γ i , (dxi )# , where the quantities on the right
hand side are defined in exterior algebra of differential forms (the canonical musical isomorphism
is used, as we assume the existence of some metric). The technical details of transforming Ni2
are in the Appendix F.

4.4

Material Matrices

The overview: As we discussed before, the discretized Maxwell’s equations are the same in FIT
and Whitney FEM on the primal/dual cell complex. From computational point of view, the
only difference is in the discretization of material laws. In this section we present our extension
of FIT material matrices and the energetic approach generalised to the space-time setting.
In both approaches the error is introduced in this step. In the following we assume that, the
material is allocated in primal grid, i.e., in each 4D cell K4i there is a constant material mapping
ξi .
4.4.1

4D FIT with Whitney Interpolation

The overview: Since we work with a non-orthogonal mesh pair, the traditional FIT approach for
creating diagonal material matrices is not possible. The construction of FIT material matrices
on non-orthogonal grids has been approached by several authors. The closest approach to
ours is described in [21], where the construction of material matrices for (non-orthogonal) 2D
quadrilateral meshes in FIT framework was proposed. However, it is not clear how this work can
be extended to space-time. Thus, we rather use our general STA discretization and Whitney
interpolation, which are ready to be applied in space-time. We would like to mention that
the sparsity pattern of the material matrices are exactly the same as in [21] 3 and are sparser
than in 4D energetic approach. The error is introduced in (75) by assuming that the fields are
e j . This assumption is needed, to recover the fields by knowing
constant on every dual facet K
2
e j . Then ξ can be used to relate the fields, resulting in the
only the DoFs, i.e., integrals over K
2
relation between DoFs. This is done with an approximation error.
We recall cf. (39) that the dual DoFs are defined as
Z
Z
2
(74)
Igj = (d x) ∧ G = (d2 x) ∧ ξ(F ) .
ej
K
2

ej
K
2

In order to construct the material matrices, we express F in terms of the primal DoFs fj and
use the formula above.
We follow the general FIT approach for approximating the integral in (74); namely, approximate it by taking the value of F at only one point xj . It may be thus perceived as an extension
of the one-point (midpoint?) quadrature rule to the geometric calculus. Thus, we approximate
(74) by


Z
Z


fj ∧ ξ(Fj ) ,
(d2 x) ∧ ξ(F ) ≈  (d2 x) ∧ ξ(Fj ) = W
(75)
ej
K
2

ej
K
2

3 Of course, to compare the 3D method [21] and our 4D one, one needs to either extend the first one to 4D
or more straightforwardly reduce our scheme to 3D in the case of no rotation Ω = 0.

16

e j with K j (or the barycenter of
where Fj := F (xj ) should be the field at the intersection of K
2
2
j
e according to midpoint rule). ξ is ξ averaged over the dual facet K
e j . However, due to the
K
2
2
fact that it is easy to interpolate the field at the barycenter of the primal facet K2j we have
used that point, see Fig. 3. Therefore, we are approximating the integral in (74) by the value
e j.
of F at the point which is not necessary lying on K
2
The mapping (143)–(146) is easy to calculate, i.e., having given (t̄, x̄, ȳ, z̄) it is easy to
calculate (t, x, y, z). However, it is difficult to invert, i.e., having given (t̄, x̄, ȳ, z̄) it is not
possible to obtain an analytic expression for (t, x, y, z); one needs to rather go for numerical
techniques to do that. Although it is difficult to invert it is easy to “guess” (t̄, x̄, ȳ, z̄) of
some interesting points, e.g., barycenters. Being more precise: barycenters are mapped to
barycenters. Therefore, if one is interested in the interpolated fields at the barycenter of K2j
which corresponds to, e.g., 2D facet ∆1 in the reference tesseract with vertices (±1, ±1, −1, −1),
then its barycenter is given by (t̄, x̄, ȳ, z̄) = (0, 0, −1, −1). The same holds true for barycenters
of edges, 2D facets, 3D volumes and the whole tesseract.
Using Whitney interpolation
X
F (x) =
fk Nk2 (x) ,
(76)
k

we calculate
Z
X Z
X
2
fj,i ∧ ξi (F (xj )) =
Igj = (d x) ∧ ξ(F ) =
(d2 x) ∧ ξi (F ) ≈
W
i

ej
K
2

i

e j,i
K
2

!
=

X
i

fj,i ∧ ξi
W

X

fk Nk2 (xj )

k

=

"
X X
k

#
fj,i ∧
W

ξi Nk2 (xj )



fk

(77)

i

and thus, we are led to the material matrix Mξ
X

fj,i ∧ ξi N 2 (xj ) .
[Mξ ]jk := I −1
W
k

(78)

i

e j,i is the part of K
e j contained in K i .
Above K
4
2
2
In the formula above the averaging is done explicitly. Therefore, it might be used to define
the averaged ξ j in (75) by requiring that it satisfies
fj ∧ ξ j (Fj ) =
W

X

fj,i ∧ ξi (Fj ) .
W

(79)

i

Due to the simplicity of predicting reference position of barycenters, we have used xj :=
barycenter of ∆j . The choice of xj ∈ K2j does not influence the fact that the obtained material
matrix Mξ is non-symmetric and thus endangers stability. However, this issue can be heuristically solved by taking only the symmetric part of Mξ as the material matrix. The numerical
experiments (without rotation) show that such operation does not influence the convergence of
the scheme, while making it stable.
4.4.2

4D FEM

The overview: We extend here to space-time the approach in [3]. The error is introduced in
(86) by assuming that the field G is constant in each 4D element K4i .
Maxwell’s equations can be obtained by applying variational principle δS = 0 to the action
integral, e.g., [24, Sections 2.4 and 2.5]
Z
1
S=
F · G|d4 x| .
(80)
2
K4

Since we have already discretized the operators, we will use S only to derive discrete material
equations.
17

We define the action integral with explicitly applied material laws
Z
1
SF =
F · ξ(F )|d4 x| .
2

(81)

K4

The core idea is to spot that
G=2

δSF
δS
=
= ξ(F ) ,
δF
δF

(82)

where δ/δF is partial functional derivative with respect to F . Indeed, the variation of S with
respect to F is
Z
G
δS =
· δF |d4 x| ,
(83)
2
K4

and the variation of SF
Z
Z
Z
1
1
4
4
δSF =
δF · ξ(F )|d x| +
F · ξ(δF )|d x| =
ξ(F ) · δF |d4 x| ,
2
2
K4

K4

(84)

K4

where we used the fact that the material mapping is self-conjugate, i.e., B1 · ξ(B2 ) = ξ(B1 ) · B2
for all bivectors B1 and B2 .
Here we extend the 3D reasoning found in [3, Section III] to our space-time setting. Approximating G by a piecewise constant field, i.e., in each element K4i it has the value Gi and
assuming the basis functions satisfy
Z
Z
X
fj,i ⇒
Nj2 |d4 x| = I −1 W
Gi · Nj2 |d4 x| = gj ,
(85)
i

K4i

K4i

we obtain (the approximation sign is due to using G(x) ≈ Gi if x ∈ K4i )
Z
Z
1X
1X
1X
2
4
fj Nj · G|d x| ≈
fj Gi · Nj2 |d4 x| =
fj gj ,
S=
2 i, j
2 i, j
2 j
K4i

(86)

K4i

from which follows
gi ≈ 2

∂S
.
∂fi

(87)

The action integral reads
Z
Z
1
1X
4
SF =
F · ξ(F )|d x| =
fi fj Ni2 · ξ(Nj2 )|d4 x| .
2
2 i, j
K4

(88)

K4

Using SF = S and comparing action integrals in (86) and (88) one gets


X Z


 Ni2 · ξ(Nj2 )|d4 x| fj =: M ξ f .
gi ≈
i
j

(89)

K4

By this construction, we have translated (82) to the discrete setting in an approximate sense
gi ≈ 2



∂S
∂SF
=
= M ξf i .
∂fi
∂fi

18

(90)

5

Resulting Scheme

The overview: In this section, we perform linear algebraic manipulations to obtain the numerical
scheme. In FIT case, (105) may be seen as an extension of leapfrog to the rotating case
or material equations, for which d and h depend on e and b as well. We have avoided any
extrapolation w.r.t. time as it may be the cause of instabilities as mentioned in [18, Appendix]
and confirmed by our numerical examples.
In general the decomposition of g = M ξ f according to e, b, d, h would read
−
−
hn = Mνb
(n)bn−1 + Mνe
(n)en−1/2 + Mνb (n)bn +
+
+
+Mνe
(n)en+1/2 + Mνb
(n)bn+1
n+1/2

d

=

−
Mεb
(n

+ 1/2)b + Mεe (n + 1/2)e

where the material matrices are defined
 

+
+
h0
Mνb Mνe
Mνb
0
−
+
 d1/2   Mεb
M
M
0
εe
εb
  −

−
+
 h1  Mνb
Mνe
Mνb Mνe
 3/2  
−
  0
 d
0
Mεb Mεe
 

−
−
 h2   0
0
M
νb Mνe
=

 ..   ..
..
..
..
 .   .
.
.
.
 

 hn−1   0
0
0
0
 

dn−1/2   0
0
0
0
hn

0

0

(91)

n

0

0

n+1/2

+

+
Mεb
(n

n+1

+ 1/2)b

,

(92)

according to
0
0
+
Mνb
+
Mεb
Mνb
..
.

...
...
...
...
...
..
.

0
0
0
0
0
..
.

0
0
0
0
0
..
.

0
0
0

...
...
...

Mνb
−
Mεb
−
Mνb

+
Mνe
Mεe
−
Mνe

0
0
0
0
0
..
.



b0
e1/2
b1
e3/2
b2
..
.

















.








+   n−1 
Mνb   b

+   n−1/2 
Mεb
e
bn
Mνb

(93)

+
Being more explicit, the block Mνe
(n) is given by
+
Mνe
(n) := Mξ (n(nf + ne ) + nf + 1 : (n + 1)(nf + ne ), n(nf + ne ) + 1 : n(nf + ne ) + nf ) (94)

where nf and ne are, respectively, the number of facets and edges in the reference mesh. The
other blocks are defined analogously.
In general all material matrices depend on the time step n. However, if the angular velocity
Ω is constant and we use a constant time step ∆t, then all the matrices are constant in time.
This comes from the fact that the 4D mesh has certain symmetry. Namely, the layer of cells
below (=before) and above (=after) time n has the same geometry as the layer corresponding
to n + 1. Thus in what follows we skip the dependency on n to simplify the notation.
We can recover 3D FIT with leapfrog time integration in our framework when the angular
velocity Ω = 0. The relation between our DoFs and the ones used in [25] is (the signs come
from the conventions employed in this paper)
_
_

d = −d

_
_

e = −∆t_
e

b= b

_

h = −∆t h .

(95)

Then the time step looks as follows
hn = Mνb bn
e n
dn+1/2 = dn−1/2 + Ch
−1 n+1/2
en+1/2 = Mεe
d

bn+1 = bn + Cen+1/2 ,

(96)

which is usually written in more compact form as
−1 e
en+1/2 = en−1/2 + Mεe
CMνb bn

b

n+1

n

n+1/2

= b + Ce

.

(97)
(98)

The Maxwell’s grid equations (51) do not change with the geometry of the mesh. Only
material matrices are affected by this change. Therefore, when Ω 6= 0 then we obtain (the first
19

and the third line are (91) and (92), respectively)
− n−1
+ n+1
− n−1/2
+ n+1/2
hn = Mνb
b
+ Mνe
e
+ Mνb bn + Mνe
e
+ Mνb
b
n+1/2
n−1/2
n
e
d
=d
+ Ch
h
i
− n
+ n+1
−1
en+1/2 = Mεe
dn+1/2 − Mεb
b − Mεb
b

bn+1 = bn + Cen+1/2 .

(99)
(100)
(101)
(102)

−
+
For 4D FIT Mνb
= Mνb
= 0. This comes from the fact that the basis functions associated
n+1
n−1
with b
and b
vanish at the time step n. Therefore, there is no contribution of bn+1 and
bn−1 to hn . This is the manifestation of the fact that our extension of FIT has narrower stencil.
For 4D FEM symmetry of Mξ translates to the following symmetries of M ’s



−
+
−
+ T
T
T
+ T
− T
Mνb = Mνb
Mεe = Mεe
Mεb
= Mνe
Mεb
= Mνe
Mνb
= Mνb
.
(103)

To keep the workflow, i.e., keep the time integrator scheme explicit, as in stationary case,
one would need to, e.g., extrapolate en+1/2 and bn+1 in (99) and (101). This however, leads to
instabilities as we will demonstrate later.
Therefore, we derive an implicit scheme. This can be done by, e.g., by plugging (102), (100)
and (99) into (101). After some algebra one arrives at
i
h
(104)
en+1/2 = M −1 γbn + βen−1/2 + αbn−1
bn+1 = bn + Cen+1/2 ,

(105)

where
+
e νe
e + C + M +C
M :=Mεe − CM
− CM
νb
εb
e −
α :=M − + CM
εb

νb
−
e
CMνe

β :=Mεe +
e νb + CM
e + − M−
γ :=CM
νb
εb

(106)
(107)
(108)
(109)

The good news are that we need to solve only one system of equations (corresponding to the
inversion of M ) at each time step, which was also already the case for 3D FIT with Ω = 0
(inversion of Mεe ).

6

Numeric Simulation of a Rotating Ring Resonator

As a research example we consider the rotating ring resonator depicted in Fig. 7. We assume
that it is filled with vacuum, and boundary is made of PEC. First we calculate the eigenmodes
of the non-rotating structure using 3D FIT in frequency domain [25]. These are collected in
Fig. 8. Afterwards they are used as initial values for the time-domain simulation of the rotating
structure.
According to the discussion in Section 4.2 the d and b values correspond to the stationary
observer. Therefore, they can be passed directly (without any processing) from the non-rotating
to the rotating solver. This is interpreted as exciting the mode field pattern in a stationary
frame, while the PEC boundary is rotating. One may perceive it that the stationary mode
~ s (x) and B
~ s (x) is transformed to the rotating reference frame via local Lorentz transpattern E
~ r (x) 6= E
~ s (x) and B
~ r (x) 6= B
~ s (x).
formations, yielding E
The other choice would be to insist that the mode pattern should be excited in rotating
~ r (x) = E
~ s (x) and B
~ r (x) = B
~ s (x). This
reference frame, i.e., should not be transformed, i.e., E
can be achieved by taking
Z
q
∆t
1/2
_
ei =
ei
(110)
ei := dτ _
e i (τ ) ≈ ∆τ _
e i (0) = 1 − vi2 ∆t_
cosh(ri Ω)
Z
q
∆t
_
_
_
_
h0i := dτ h i (τ ) ≈ ∆τ h i (0) = 1 − vi2 ∆t h i =
hi ,
(111)
cosh(ri Ω)
20

b
a

Figure 7: Geometry of the resonator. a = 5mm, b = 10mm. The height in z dimension is 2mm.

(a) m = 0

(b) m = 1

(c) m = 2

(d) m = 3

(e) m = 4

(f) m = 5

Figure 8: Considered eigenmodes

21

where vi and thus ri are values of 3D velocity and radial coordinate at, e.g., the middle of, the
_
edge corresponding to _
e i or h i .
Depending on the setting of a particular physical situation one may freely chose whether
the fields should be excited in the rotating or stationary reference frame by transferring d and
b or e and h as initial conditions.

6.1

Analytic Predictions

The eigenmodes have harmonic azimuthal and time dependency
z
Estat
= E z (r, θ, z, t) = E z (r, z) cos(kθ θ) cos(ωt) .

(112)

They may be perceived as two clockwise and counterclockwise modes
E z (r, z) cos(kθ θ) cos(ωt) =
1
+
−
+ Estat
= E z (r, z) [cos(kθ θ + ωt) + cos(kθ θ − ωt)] = Estat
, (113)
2
where
1 z
E (r, z) cos(kθ θ + ωt)
2
1
:= E z (r, z) cos(kθ θ − ωt) .
2

+
Estat
:=

(114)

−
Estat

(115)

According to [23] these two modes should decouple if the structure is rotating
1 z
E (r, z) cos(kθ θ + ω+ t) =
2
1
:= E z (r, z) cos(kθ θ − ω− t) =
2

+
Erot
:=
−
Erot

1 z
E (r, z) cos(kθ θ + ωt + δωt)
2
1 z
E (r, z) cos(kθ θ − ωt + δωt) .
2

(116)
(117)

Thus imposing initial value E z (r, θ, z, 0) is equivalent to exciting these modes with the same
amplitude. The solution therefore is
+
−
z
Erot
(r, θ, z, t) = Erot
+ Erot
=
1 z
= E (r, z) [cos(kθ θ + ωt + δωt) + cos(kθ θ − ωt + δωt)] =
2
= E z (r, z) cos(kθ θ + δωt) cos(ωt) . (118)

Therefore, if we chose a reference point (r, θ, z) = (r0 , 0, z0 ) and E z (r0 , z0 ) 6= 0, we should
observe the time signal
z
Erot
(r0 , 0, z0 , t)
= cos(δωt) cos(ωt) .
(119)
E z (r0 , z0 )
According to [23, Equation (4.5)] the rotation induced frequency shift should be
δω ≈ mΩ ,

(120)

where m is the number of wavelengths in azimuthal direction.

6.2

Comparison with Analytic Formula

According to (119) the field should look like in Fig. 9. At time T0 the amplitude is approximately
zero and this point correspond to the quarter of beating period Td . Therefore, the frequency
shift is
2π
π
δω =
=
.
(121)
Td
2T0
The frequencies of the discrete time-domain signal are extracted using one of the state of the
art algoriths described in [26].
22

Ez
10

Ez

5

0

t
Td 4

5

10
10

5

0

5

10

~ = [0, 0, E z ] [σx , σy , σz ]T at the sample point. Td is the
Figure 9: Behavior of the electric field E
beating period.

vmax /c

m

0
1
2
3
4
5

< 100%
Null
95.7%
95.8%
95.8%
96.%
96.2%

99.63% 30.42%
Null
Null
58.1%
2.2%
58.5%
2.8%
59.2%
2.6%
60.6%
4.6%
62.2%
5.%

3.14%
Null
0.
0.9%
2.%
1.3%
8.2%

0.31% 0.03%
Null
Null
0.2%
3.2%
0.9%
1.%
3.9%
3.4%
1.3%
2.4%
8.2%
2.6%

Table 2: FIT case: the relative error (122) vs. mode number and rotation rate/velocity of the
outer rim of the ring

vmax /c

m

0
1
2
3
4
5

< 100%
Null
95.6%
95.7%
95.7%
95.8%
95.8%

99.63% 30.42%
Null
Null
57.4%
3.5%
57.8%
0.
58.%
0.1%
58.4%
0.7%
58.7%
2.6%

3.14%
Null
0.4%
0.8%
5.2%
1.1%
13.%

0.31% 0.03%
Null
Null
6.7%
5.4%
3.5%
3.9%
4.9%
5.%
3.4%
7.9%
12.9% 2.4%

Table 3: FEM case: the relative error (122) vs. mode number and rotation rate/velocity of the
outer rim of the ring

23

Figure 10: Comparison of analytical and numerical frequency shifts in case of non-relativistic
velocities.
The obtained relative difference with the analytic formula (120)
|δω simulation − mΩ|
mΩ

(122)

for various m and Ω is presented in Tables 2 and 3.
As we can see from Tables 2 and 3 the numerical results are in close agreement with the
analytic formula (120) for velocities < 30% speed of light. This can also be seen from Fig. 10.
This is expected as this formula was derived assuming that the rotation rate Ω is low, i.e.,
rΩ/c  1. We want to stress that the two solutions: analytic and numerical one suffer from
two different types of error. The analytic solution suffers from the modelling error, namely,
absence of relativistic effects. The numerical solution was obtained without any non-relativistic
assumptions. Therefore, it suffers only from the error coming from the discretization.

6.3

Extrapolation Leads to Instabilities

We consider the mode with m = 3. The results for zeroth, first and second order extrapolation
in our 4D extension of FIT are depicted in Fig. 11. We can see that although first order
extrapolation is more stable than zeroth one, increasing the extrapolation order to two gives
less stable result. Therefore, suggestion in [19] to apply more accurate extrapolation is not an
option as it leads to more unstable scheme.
Applying an implicit scheme with FIT matrices leads to unstable scheme as we discussed
it before. However, it remains stable for a few millions of time steps as presented in Fig. 12.
Therefore, extrapolation increases the overall instability of the scheme.
Analogous results for FEM are in Fig. 13. We see that the zeroth order gets immediately
unstable. First order extrapolation “explodes” after 9000 time steps. Second order extrapolation although stable leads to the dissipation, which should not be present. Moreover, there is
no guarantee that for other geometries or initial values this extrapolation will lead to a stable
scheme.
The implicit FEM scheme remains stable up to 50 millions of time steps as shown in Fig. 14.

7

Conclusion and Outlook

Starting from the space-time setting we have proposed numerical scheme for simulation of
electromagnetic fields in arbitrary non-inertial reference frame. We have not made any lowΩ v
velocity assumptions, i.e., rΩ
c , ω , c  1, at any point, thus our approach can also be used for
a simulation of, e.g., linear accelerator or synchrotron radiation.
The application of a space-time method ensures that the obtained scheme is consistent with
relativity. Moreover, it enhances the efficiency of the scheme in the sense that the material
24

(a) Zeroth order

(b) First order

(c) Second order

Figure 11: Amplitude of the field for various extrapolations (FIT).

25

Figure 12: Energy conservation for FIT
matrices for a rotating observer are constant in time. The introduced DoFs allow to choose
whether the fields should be transformed or not when passing them from a stationary to a
rotating solver. Our DoFs are also convenient in terms of enforcing boundary conditions. Not
only PEC as discussed in the paper, but also, e.g., one may be interested in the implementation
of Absorbing Boundary Conditions (ABCs). Typically, they are derived by considering outgoing
plane waves, which are solutions defined in a stationary reference frame. Consequently, one
would prefer to relate to d and b as they relate to the stationary observer. For example,
extension of Mur’s ABCs to a 3D Cartesian grid rotating with arbitrary high rotation rate is
straightforward; namely, one needs to replace e by d in [17, Eq. 15-16].
Numerical experiments show the necessity of using an implicit scheme (instead of extrapolation of DoFs) in order to keep stability of the scheme. Our results both for 4D FEM and FIT
are in good agreement with the analytic predictions for low-velocities.
Acknowledgement: The work of the first, second and third author is supported by the
’Excellence Initiative’ of the German Federal and State Governments and the Graduate School
of CE at Technische Universitaet Darmstadt.

A

Spacelike, Lightlike and Timelike


 timelike if
The vector v is called lightlike if


spacelike if

v>0
v=0
v<0

(123)

This notion naturally extends to n-vectors. We will refer only to 2D facets in this context. The
2D facet ∆ is called timelike or lightlike if it contains a timelike or lightlike vector, respectively;
otherwise it is called spacelike. Mathematically, one can test the type of ∆ by squaring the
bivector W associated with it; namely


W2 > 0
 timelike if
The facet ∆ is called lightlike if
(124)
W2 = 0


2
spacelike if W < 0

26

(a) Zeroth order

(b) First order

(c) Second order

Figure 13: Amplitude of the field for various extrapolations (FEM).

27

Figure 14: Energy conservation for FEM

B

Velocity Field of a Rotating Observer

For simplicity we set c = 1, obtaining
dp
= [1, − sin(θ) tanh(rΩ), cos(θ) tanh(rΩ), 0][γt , γx , γy , γz ]T
dt
dp
1
=
dt
cosh(rΩ)

dp
dp
u=
=
dt
dt
=[cosh(rΩ), − sinh(rΩ) sin(θ), sinh(rΩ) cos(θ), 0][γt , γx , γy , γz ]T
ux
dx
=
= − tanh(rΩ) sin(θ)
vx =
dt
ut
dy
uy
vy =
=
= tanh(rΩ) cos(θ)
dt
ut
q
v = vx2 + vy2 = tanh(rΩ)

C

(125)
(126)

(127)
(128)
(129)
(130)

Whitney Interpolation

Thereby, (54) requires that

X

∇ ∧ (W n An ) = W n+1 (Dn An ) ,
X
∇ ∧ (Ain Nin (x)) =
([Dn ]ji Ain )Njn+1 (x) ,

i

(131)
(132)

i, j

X

Ain ∇ ∧ Nin (x) =

i

X

Ain [Dn ]ji Njn+1 (x) ,

(133)

[Dn ]ji Njn+1 (x) .

(134)

i, j

∇ ∧ Nin (x) =

X
j

One can indeed verify that Whitney interpolation is compatible with the approximationfree discretization introduced before; namely, they lead to the same system of equations as we
prove next. In order to obtain the discrete system of equations, we plug the reconstructed field
F = W 2 f in ∇ ∧ F = 0
!
∇∧F =

X
i

fi ∇ ∧

Ni2 (x)

=

X

X X 
 
fi D2 ji Nj3 (x) =
D2 ji fi

i, j

j

28

i

Nj3 (x) = 0 ,

(135)

and since the basis functions Nj3 (x) are linearly independent, it is equivalent to
X 
D2 ji fi = 0 ⇔ D2 f = 0 ,

(136)

i

which is exactly (40)-left.

D

Explicit Formulae for Whitney Elements

Let us now elaborate on the explicit construction of the basis functions. Since we consider one
4D element, there is only one basis function associated with a 4-vector, that is
(4)

N14 := I ,

(137)

I := γ t ∧ γ x ∧ γ y ∧ γ z .

(138)

where
(4)

Similarly, there are eight 3D facets and thus eight basis functions associated with 3-vectors,
viz.
(4)

3
Ni±
:= ±γi I (1 ± xi ), i = t, x, y, z .

(139)

We can check that (134) with n = 3 holds and the rightmost part of the deRham complex (54)
commutes, namely
!
(4)
(4)
(4)
X
X
3
k
i
∇ ∧ Ni± =
(γ ∂k ) ∧ ±γi I (1 ± x ) =
γ k ∧ (γi I )δki = I = N14 .
(140)
k

k

There are 24 2D facets and thus 24 basis functions associated with 2-vectors, viz.
2
Ni±
:=
1 , j±2

(4)
1
γi ∧ γj I (1 ±1 xi )(1 ±2 xj ) .
16

(141)

We can verify with (134) and n = 2 that (141) are indeed Whitney forms
i
(4) h
X 1
γ k ∧ (γi ∧ γj I ) ±1 δki (1 ±2 xj ) ±2 δkj (1 ±1 xi ) =
16
k
"
#
(4)
(4)
1
j
i
3
3
±1 γj I (1 ±2 x ) ∓2 γi I (1 ±1 x ) = ±1 ±2 Nj±
=
∓2 ±1 Ni±
. (142)
2
1
16

2
∇ ∧ Ni±
=
1 , j±2

In a similar fashion, we can define N 0 and N 1 . Nevertheless, since we never refer to them,
we do not state them explicitly.

E

Derivation of Φi

The reference tesseract described above has to be mapped to a deformed tesseract in the physical
domain. This is accomplished with the transformation
t=

16
X

Ti (t̄, x̄, ȳ, z̄)ti

(143)

Ti (t̄, x̄, ȳ, z̄)xi

(144)

Ti (t̄, x̄, ȳ, z̄)yi

(145)

Ti (t̄, x̄, ȳ, z̄)zi ,

(146)

i=1

x=

16
X
i=1

y=

16
X
i=1

z=

16
X
i=1

29

where (t̄, x̄, ȳ, z̄) are the coordinates in the reference tesseract, and (t, x, y, z) those in the physical one, and (ti , xi , yi , zi ) are the physical coordinates of the i-th vertex. The functions Ti are
given by
1
Ti (t̄, x̄, ȳ, z̄) =
(1 + t̄i t̄)(1 + x̄i x̄)(1 + ȳi ȳ)(1 + z̄i z̄) ,
(147)
16
where (t̄i , x̄i , ȳi , z̄i ) = (±1, ±1, ±1, ±1) are the reference coordinates of the i-th vertex.
This transformation comes from the following reasoning, compare [13, Section 8.1.3]. We
postulate the map between physical coordinates x̄i and the reference ones xi
x̄m = Am +

4
X
i

i
Am
i x +

4
X

4
X

i j
Am
ij x x +

i<j

i j k
Am
ijk x x x +

i<j<k

4
X

i j k l
Am
ijkl x x x x .

(148)

i<j<k<l

m
m
m
There are in total 4×(1+4+6+4+1) = 64 parameters associated with Am , Am
i , Aij , Aijk , Aijkl .
Then we impose that this mapping should map the nodes of the reference tesseract to the
physical ones. This accounts for 64 linear equations as there are 16 nodes with 4 coordinates
each. This linear system is composed of 64 equations with 64 unknowns, and can be solved for
A’s in terms of positions x̄ij of the nodes in the physical tesseract. By substituting the solution
for A’s in (148), we obtain the transformations stated in (143)–(146).
For illustration purposes we discuss here the properties of the linear interpolation in 1D.
Coordinate transformation is specified by the function x(x̄)

x=

2
X

Ti (x̄)xi ,

(149)

i=1

with

1
(1 + x̄i x̄) ,
2
where x̄i = ±1. Being more explicit, we write
Ti =

1
x̄
1
1
x = x1 (1 − x̄) + x2 (1 + x̄) = (x1 + x2 ) + (x2 − x1 ) .
2
2
2
2

(150)

(151)

The mapping (149) maps the reference domain x̄ ∈ [−1, +1] to the physical domain x ∈ [x1 , x2 ].
Moreover, the boundary of the reference domain, e.g., the point x̄ = −1, is mapped to the
boundary of the physical domain, e.g., x = x1 . Also the reference barycenter x̄ = 0 is mapped
to the physical barycenter x = 21 (x1 + x2 ).
The linear interpolation is extended to 4D by interpolating four coordinate functions xj (t̄, x̄, ȳ, z̄)
defined on a tensor product of the 1D reference interval, that is Ξ = [−1, +1] ⊗ [−1, +1] ⊗
[−1, +1] ⊗ [−1, +1]. We note that not only the barycenter of the reference tesseract 0 ⊗ 0 ⊗ 0 ⊗ 0
is mapped to the barycenter of the physical 4D cell, but also the barycenters of the lower
dimensional reference elements, e.g., the barycenter 1 ⊗ 1 ⊗ 0 ⊗ 0 of the 2D facet 1 ⊗ 1 ⊗
[−1, +1] ⊗ [−1, +1], are mapped to the barycenters of the corresponding physical elements Kni .

F

The Transformation to the Physical Domain

Manipulations below correspond to calculating the pullback of a differential form with respect
to a diffeomorphism. The Ni2 should transform as a 2-form. For example in our implementation,
for the function
1
2
Nx+,y+
=
γ̄x ∧ γ̄y I(1 + x)(1 + y)
(152)
16
we have used
γ̄x ∧ γ̄y I = γ̄x · (γ̄y · I) = γ̄x · (γ̄ t ∧ γ̄ x ∧ γ̄ z ) = −γ̄ t ∧ γ̄ z = (∇z̄) ∧ (∇t̄) .

(153)

We need to express ∇x̄α in terms of basis vectors γ α of physical space
γ̄ α =

∂ x̄α β
γ .
∂xβ
30

(154)

However, since we do not have an analytic expression for the inverse mapping x̄α (xβ ), we need
to use
  α
∂ x̄α
∂x
=
inv
,
(155)
β
∂x
∂ x̄ β

 α
∂xα
where inv ∂∂xx̄ is the inverse of the Jacobi matrix with the entries ∂x
∂ x̄ β = ∂ x̄β . After making
all this substitutions we obtain Ni2 expressed in terms of the basis vectors associated with
Minkowski space-time and its metric, but the reference coordinates are used.

References
[1] B. Auchmann and S. Kurz. Observers and splitting structures in relativistic electrodynamics. Journal of Physics A: Mathematical and Theoretical, 47(43):435202, 2014.
[2] A. Bossavit. Computational Electromagnetism . Academic Press, San Diego, 1998.
[3] L. Codecasa, R. Specogna, and F. Trevisan. Symmetric Positive-Definite Constitutive Matrices for Discrete Eddy-Current Problems. IEEE Transactions on Magnetics, 43(2):510–
515, Feb 2007.
[4] L. Codecasa and F. Trevisan. Piecewise uniform bases and energetic approach for discrete
constitutive matrices in electromagnetic problems. International Journal for Numerical
Methods in Engineering, 65(4):548–565, 2006.
[5] C. Doran and A. Lasenby. Geometric Algebra for Physicists. Cambridge University Press,
Cambridge, second edition, 2003.
[6] M. Fecko. Differential Geometry and Lie Groups for Physicists. Cambridge books online.
Cambridge University Press, 2006.
[7] F. W. Hehl. History Related Maxwell’s equations in Minkowski’s world: their premetric
generalization and the electromagnetic energy-momentum tensor. Annalen der Physik,
17(9-10):691–704, 2008.
[8] D. Hestenes. Space-time algebra. Documents on modern physics. Gordon and Breach, 1966.
[9] D. Hestenes. Differential Forms in Geometric Calculus. In F. Brackx, R. Delanghe, and
H. Serras, editors, Clifford Algebras and their Applications in Mathematical Physics, pages
269–285. Springer, 1993.
[10] D. Hestenes. The Shape of Differential Geometry in Geometric Calculus. In Leo Dorst and
Joan Lasenby, editors, Guide to Geometric Algebra in Practice, pages 393–410. Springer,
2011.
[11] D. Hestenes and G. Sobczyk. Clifford algebra to geometric calculus : a unified language
for mathematics and physics. D. Reidel ; Distributed in the U.S.A. and Canada by Kluwer
Academic Publishers, Dordrecht; Boston; Hingham, MA, U.S.A., 1984.
[12] R. Hiptmair. Finite elements in computational electromagnetism. Acta Numerica, 11:237–
339, 1 2002.
[13] J.-M. Jin. The finite element method in electromagnetics. Wiley, New York, 2002.
[14] M. Klimek, U. Roemer, S. Schoeps, and T. Weiland. Space-Time Discretization of
Maxwell’s Equations in the Setting of Geometric Algebra. In 2013 International Symposium on Electromagnetic Theory, editor, IEEE Xplore, 2013.
[15] C. Mattiussi. The Geometry of Time-Stepping. Progress In Electromagnetics Research,
32:123–149, 2001.
[16] C.W. Misner, K.S. Thorne, and J.A. Wheeler. Gravitation. W.H. Freeman and Company,
1973.
31

[17] G. Mur. Absorbing Boundary Conditions for the Finite-Difference Approximation of the
Time-Domain Electromagnetic-Field Equations. IEEE Transactions on Electromagnetic
Compatibility, EMC-23(4):377–382, Nov 1981.
[18] R. Novitski, J. Scheuer, and B. Z. Steinberg. Unconditionally stable finite-difference timedomain methods for modeling the Sagnac effect. Phys. Rev. E, 87:023303, Feb 2013.
[19] C. Peng, R. Hui, X. Luo, Z. Li, and A. Xu. Finite-difference time-domain algorithm for
modeling Sagnac effect in rotating optical elements. Opt. Express, 16(8):5227–5240, Apr
2008.
[20] J. Salamon, J. Moody, and M. Leok. Geometric Representations of Whitney Forms and
their Generalization to Minkowski Spacetime. arXiv, Feb. 2014.
[21] R. Schuhmann and T. Weiland. A Stable Interpolation Technique for FDTD on Nonorthogonal Grids. International Journal of Numerical Modelling: Electronic Networks, Devices
and Fields, 11(6):299–306, 1998.
[22] G. Sobczyk. Simplicial Calculus with Geometric Algebra. In Clifford Algebras and Their
Applications in Mathematical Physics, pages 279–292. Springer, 2011.
[23] B. Z. Steinberg, A. Shamir, and A. Boag. Two-dimensional Green’s function theory for
the electrodynamics of a rotating medium. Phys. Rev. E, 74:016608, Jul 2006.
[24] A. Stern, Y. Tong, M. Desbrun, and J. E. Marsden. Geometric computational electrodynamics with variational integrators and discrete differential forms. arXiv:0707.4470v2,
2009.
[25] T. Weiland. Time Domain Electromagnetic Field Computation with Finite Difference
Methods. International Journal of Numerical Modelling, 9:295–319, Jul. 1996.
[26] Y. Yang and M. Pesavento. A unified successive pseudo-convex approximation framework,
2015. cite arxiv:1506.04972v2.pdfComment: submitted to IEEE Transactions on Signal
Processing; original title: A Novel Iterative Convex Approximation Method.

32

