A Space-Time Approach for the Time-Domain Simulation in a Rotating Reference Frame

TL;DR

This paper develops a space-time discretisation method for Maxwell's equations in rotating frames, enabling accurate simulation of phenomena like the Sagnac effect, and compares the stability of FIT and FEM approaches.

Contribution

It introduces a direct space-time discretisation approach for Maxwell's equations in rotating frames, with a focus on stability and convergence analysis of FIT and FEM methods.

Findings

The scheme accurately models the Sagnac effect in a rotating resonator.

Comparison shows differences in stability between FIT and FEM approaches.

The method avoids non-relativistic assumptions in simulations.

Abstract

We approach the discretisation of Maxwell's equations directly in space-time without making any non-relativistic assumptions with the particular focus on simulations in rotating reference frames. As a research example we study Sagnac's effect in a rotating ring resonator. After the discretisation, we express the numerical scheme in a form resembling 3D FIT with leapfrog. We compare the stability and convergence properties of two 4D approaches, namely FIT and FEM, both using Whitney interpolation.

Figures

Tables

Table 1. TABLE I: FIT case: the relative difference η 𝜂 \eta vs. mode number and rotation rate/velocity of the outer rim of the ring

Equations

$$\nabla \wedge F=0,$$

$$G=\xi (F).$$

$$\vec{D}=\epsilon \vec{E}$$

$$\xi (F)=\frac{1}{2}\left[\left(\epsilon +\frac{\nu }{c^2}\right)F-\left(\epsilon -\frac{\nu }{c^2}\right)uFu\right],$$

$$p_t({\vec{r}}_{\text{ref}})=[t,r\cos (\theta ),r\sin (\theta ),z],$$

$$\underset{\partial K_3^i}{\oint }(d^{2}x)\cdot F=0,$$

$$f_j:=\underset{K_2^j}{\int }(d^{2}x)\cdot F,$$

$$f_{j}=:\begin{cases}e^{n+1/2}_{l}&\mbox{if $K_{2}^{j}$ is timelike}\,,\\ b^{n}_{m}&\mbox{if $K_{2}^{j}$ is spacelike}\,,\end{cases}$$

$$g_i=\sum _j{\left[M_{\xi }\right]}_{ij}{f}_j.$$

$${\left[M_{\xi }^{\text{FEM}}\right]}_{jk}:=\underset{K_4}{\int }{N}_j^2\cdot \xi (N_k^2)|d^{4}x|.$$

$${\left[M_{\xi }^{\text{FIT}}\right]}_{jk}:=I^{-1}{\tilde{W}}_j\wedge {\bar{\xi }}_j\left(N_k^2(x_j)\right),$$

$$M_{\nu b}^{\text{FIT,sym}}$$

$$M_{\epsilon e}^{\text{FIT,sym}}$$

$$M_{\epsilon b}^{\pm \text{FIT,sym}}$$

$$M_{\nu e}^{\pm \text{FIT,sym}}$$

$$M_{\nu b}^{\pm \text{FIT,sym}}$$

$$h^n$$

$$+M_{\nu e}^{+}{e}^{n+1/2}+M_{\nu b}^{+}{b}^{n+1}$$

$$d^{n+1/2}$$

$$e^{n+1/2}$$

$$b^{n+1}$$

$$\left[\begin{array}{c}{b}^{n+1}\\ e^{n+1/2}\\ b^n\end{array}\right]=U^n\left[\begin{array}{c}{b}^1\\ e^{1/2}\\ b^0\end{array}\right],$$

$$U:=\left[\begin{array}{ccc}1+C{M}^{-1}\gamma & C{M}^{-1}\beta & C{M}^{-1}\alpha \\ M^{-1}\gamma & M^{-1}\beta & M^{-1}\alpha \\ 1& 0& 0\end{array}\right],$$

$$M:=$$

$$\alpha :=$$

$$\beta :=$$

$$\gamma :=$$

$$\eta :=\frac{|\delta {\omega }^{\text{sim}}-\delta {\omega }^{\text{anal}}|}{|\delta {\omega }^{\text{anal}}|}$$

$$\Vert w{\Vert }_2:=\sqrt{\underset{t_{\text{min}}}{\overset{t_{\text{max}}}{\int }}dt{\left[w\left(t\right)\right]}^2},$$

$$\text{diff}\left(E_1^z,E_2^z\right):=\frac{\Vert E_2^z({\vec{r}}_{\text{sample}})-E_1^z({\vec{r}}_{\text{sample}}){\Vert }_2}{\Vert E_1^z({\vec{r}}_{\text{sample}}){\Vert }_2},$$

Full text

I INTRODUCTION

Concerning Clifford’s Geometric Algebra, we use the nomenclature and notation of [5], which we also perceive as a comprehensive introduction to this mathematical formalism. Maxwell’s equations [5, Eq. (7.39)] read

[TABLE]

where $F$ (or $G$) is the (dual) Faraday bivector, and $J$ the space-time current vector.

The bivectors $F$ and $G$ are related by constitutive material mapping $\xi$ defined implicitly by

[TABLE]

For linear, isotropic media $\xi$ is a space-time generalisation of the well known constitutive equations. To derive an explicit formula for $\xi$, we start with

[TABLE]

where $\varepsilon,\nu,\vec{E},\vec{H},\vec{D},\vec{B}$ are electric permittivity, magnetic reluctance, electric and magnetic field strengths, and fluxes, respectively. The resulting expression is given by

[TABLE]

with $u$ the four-velocity vector of the material, and $c$ the speed of light in vacuum.

II MOTION AS MESH’S GEOMETRY

Motion of the system is modelled by specifying Lagrangian placement map $p_{t}(\vec{r}_{\text{ref}})$, which gives a space-time position (at time $t$) of a particle with initial position $\vec{r}_{\text{ref}}$. The cylindrical coordinates of $\vec{r}_{\text{ref}}$ are denoted by $r,\varphi,z$. For a system rotating with constant angular velocity $\Omega$ around $z$-axis we use

[TABLE]

with $\theta=\varphi+\tanh\left(r\Omega/c\right)ct/r$.

We extrude a 3D mesh by applying $p_{t_{i}}$ to nodes’ positions for all time steps $t_{i}$ and connecting the 4D nodes stemming from the same 3D node; see Fig. 1. By $K_{n}^{i}$ (respectively $\widetilde{K}_{n}^{i}$) we denote the $i$-th $n$-dimensional element of the (barycentric dual) mesh.

III DISCRETIZATION

Maxwell’s equations (1) are discretised by applying their integral form to the primal/dual mesh pair, i.e.,

[TABLE]

and introducing scalar DoFs (Degrees of Freedom) on the primal mesh

[TABLE]

and analogously for the dual mesh. Next, we rename and renumber the DoFs on the primal mesh as follows

[TABLE]

where the relation between the indices $j$, $l$, $m$ and $n$ is as follows (see also Fig. 2). If the facet $K_{2}^{j}$ is timelike, it is the edge $l$ in the reference mesh $M^{3}$ extruded in time via the placement map $p_{t}$ such that $t\in\left[t_{n},t_{n+1}\right]$. Similarly, if the facet $K_{2}^{j}$ is spacelike, then it is the image of a facet with index $m$ in $M^{3}$ under the map $p_{t}$ with $t=t_{n}$. Analogously on the dual mesh, we split $g_{j}$ into $d^{n-1/2}$ and $h^{n}$.

III-A Material Equations

A discrete equivalent of $\xi$ in (4) is a material matrix $M_{\xi}$ relating discrete equivalents of $F$ and $G$, i.e.,

[TABLE]

In both methods described below, we use $n$-vector valued Whitney interpolating functions $N_{i}^{n}$ associated with $i$-th $n$-dimensional element of the mesh to derive explicit expressions for $M_{\xi}$.

III-A1 4D FEM

Applying energetic approach [4] also known as Galerkin Hodge star [3], we arrive at the material matrix with entries

[TABLE]

III-A2 4D FIT

The 4D FIT material matrix with Whitney interpolation written using Geometric Algebra reads

[TABLE]

with $\bar{\xi}_{j}$ is $\xi$ averaged over the $j$-th dual facet, and $\widetilde{W}_{j}$ the bivector associated with it, describing its magnitude and orientation in space-time. Although the collocation point $x_{j}$ may be chosen arbitrarily, we choose the barycenter of $K_{2}^{j}$ due to ease of calculation.

III-A3 Reduction to 3D Material Matrices and their Symmetrisation

The material matrix $M_{\xi}$ is split into 3D material matrices $M_{\varepsilon/\nu,e/b}^{\pm}$ to resemble a 3D time marching scheme. E.g., $M_{\varepsilon b}^{+}$ is the block of $M_{\xi}$ relating $d^{n+1/2}$ (role of $\varepsilon$) to the future (role of $+$) magnetic $b$ DoF (role of $b$), i.e., $b^{n+1}$.

While $M_{\xi}^{\text{FEM}}$ is symmetric by construction (10), the FIT material matrix $M_{\xi}^{\text{FIT}}$ is symmetrised in order to avoid instabilities, as explained in [6, especially Sec. II]. The symmetrisation of $M_{\xi}^{\text{FIT}}$ translates to the following redefinitions, where $M^{T}$ is the transpose of $M$,

[TABLE]

III-B Resulting System of Linear Equations

We have split space-time material matrix and DoFs into their 3D counterparts in a way, that the obtained numerical scheme resembles 3D FIT with leapfrog, i.e.,

[TABLE]

We would like to note, that if $\Omega=0$ then all material matrices except $M_{\varepsilon e}$ and $M_{\nu b}$ vanish, and on Cartesian grid the 4D FIT is naturally related to 3D FIT with leapfrog.

IV STABILITY

Starting from (17) we derive the recursive formula

[TABLE]

where $U$ is the update matrix

[TABLE]

where [math] and $1$ above are zero and identity matrices of proper dimensions and

[TABLE]

The solution vector $[b^{n},e^{n-1/2}]$ will stay bounded if the modulus of all eigenvalues $\lambda$ of $U$ is less or equal to unity, $|\lambda|\leq 1$. For both methods considered, the greatest $|\lambda|$ calculated on exemplary meshes used in simulations is at most $1$.

As another test of stability, we initialise the solver with randomly generated initial values and observe that the norm of the solution does not grow in time.

V CONVERGENCE

In order to verify the convergence of the scheme, we simulate the first six, $m=0\dots 5$, eigenmodes of the ring resonator structure studied in [1, Section IV], and compare the rotation induced frequency shift $\delta\omega^{\text{sim}}$ extracted from our numerical simulation with the non-relativistic semi-analytic approximation, [1, Eq. (4.5)], $\delta\omega^{\text{anal}}=m\Omega$. As depicted in Fig. 3 both approaches agree well in a non-relativistic regime.

However, in relativistic cases the relative difference

[TABLE]

between two approaches becomes significant, see Table I. We expect our solution to be correct as we have not made any non-relativistic assumptions as opposed to [1].

V-A 3D Wave Simulation without Rotation

The requirements of the convergence proof in [2] are neither fulfilled by $M_{\xi}^{\text{FIT,sym}}$ nor by $M_{\xi}^{\text{FIT}}$ (due to asymmetry). Since this is the feature of the proposed method itself, rather than its space-time extension, we focus now on 3D wave-simulation in the non-rotating ring resonator. We investigate the convergence by comparing results obtained using 4D FIT material matrix $M_{\xi}^{\text{FIT,sym}}$ with the ones using $M_{\xi}^{\text{FEM}}$, for which the proof [2] holds. We use temporal $L^{2}$ norm of a function of time $w:t\mapsto w(t)$

[TABLE]

to define the relative difference

[TABLE]

where $E^{z}_{1}(\vec{r}_{\text{sample}})$ is the time signal at the sample point $\vec{r}_{\text{sample}}$, obtained via interpolation of DoFs calculated using method and mesh indexed as $1$.

The results for a (non-)orthogonal mesh, left (right) in Fig. 4, are depicted in top (bottom) of Fig. 5.

The proof [2] can be applied to FEM independently of orthogonality of the mesh, therefore the convergence of FEM in both parts of Fig. 5 is in accordance with that theory. In case of orthogonal grid $M_{\xi}^{\text{FIT}}$ is diagonal, thus symmetric, and the proof [2] can be applied in this particular situation. Although, on non-orthogonal meshes the proof [2] cannot be repeated, we observe that $M_{\xi}^{\text{FIT,sym}}$ gives a solution that converges to the same solution as FEM (and with a similar rate of convergence).

VI CONCLUSION

We have applied space-time discretisation in two flavours, namely 4D FIT and FEM, without making any non-relativistic approximations. As a verification we recovered with good accuracy the rotation induced frequency shifts of the rotating ring resonator predicted by an alternative non-relativistic approach. The material matrix of our proposed extension of FIT is symmetrised to avoid instabilities. With taking only the symmetric part of the material matrix, the convergence is not guaranteed. However, we investigated numerically that FIT converges to the same solution as FEM, which is known to converge.

Acknowledgments

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.