⇦Chinese Physics Letters, 2019, Vol. 36, No. 1, Article code 010401 Finite Volume Time Domain with the Green Function Method for Electromagnetic Scattering in Schwarzschild Spacetime * Shou-Qing Jia (贾守卿)** Affiliations School of Computer Science and Engineering, Northeastern University, Shenyang 110819 Received 7 September 2018, online 25 December 2018 *Supported by the National Natural Science Foundation of China under Grant No 61601105.
**Corresponding author. Email: jiashouqing@neuq.edu.cn Citation Text: Jia S Q 2019 Chin. Phys. Lett. 36 010401    Abstract The finite volume time domain (FVTD) algorithm and Green function algorithm are extended to Schwarzschild spacetime for numerical simulation of electromagnetic scattering. The FVTD method in Schwarzschild spacetime is developed by filling the flat spacetime with an equivalent medium. The Green function in Schwarzschild spacetime is acquired by solving initial value problems. Both the FVTD code and the Green function code are validated by numerical results. Scattering in Schwarzschild spacetime is simulated with these methods. DOI:10.1088/0256-307X/36/1/010401 PACS:04.40.Nr, 41.20.Jb, 02.70.Dh, 02.60.Cb © 2019 Chinese Physics Society Article Text For electromagnetic-wave-based communication between spacecraft in deep space, it is important to understand the property of the electromagnetic waves in curved spacetime, for which numerical simulation has been exploited by researchers. The electromagnetic particle-in-cell (EMPIC) algorithm is used to study the physics of high-frequency electromagnetic waves in a general relativistic plasma with Schwarzschild metric.[1] In addition, the Kerr–Schild metric is incorporated into the EMPIC code for the simulation of charged particles in the region of a spinning black hole.[2] The finite difference time domain (FDTD) method and Green function method (GFM) were extended to Schwarzschild spacetime in Ref. [3]. They have also been applied to electromagnetic scattering problems (not by black holes).[4] The finite volume time domain (FVTD) method[5-7] is a numerical method that is used to simulate electromagnetic waves in flat spacetime. The advantage of the FVTD method is its flexibility when choosing a mesh. Although only a regular hexahedral mesh can be used for FDTD, many types of mesh (e.g., tetrahedral mesh) can be used for FVTD. This is beneficial when simulating fine structures. In this Letter, the FVTD method is extended to Schwarzschild spacetime in the same way as introduced in Ref. [3], which is developed by filling the flat spacetime with an equivalent medium. The FVTD and GFM are incorporated to simulate the electromagnetic scattering in Schwarzschild spacetime. When solving scattering problems by the FVTD method, the connection boundary and output boundary should be incorporated into the FVTD code.[5] The electromagnetic field on the connection boundary needs to be calculated with the Green function. Generally, the observer is far away from the scatterer. It is impossible to simulate the whole field distribution between the scatterer and observer by the FVTD method. To calculate the far field outside the FVTD domain, the electromagnetic current on output boundary should be integrated using the Green function. To calculate the Green function in curved spacetime,[8,9] we should numerically solve the differential equations.[10,11] This has been introduced in Ref. [3]. The connection boundary and output boundary in the FVTD method are introduced and the methods are validated in this work. Maxwell equations can be written in the exterior differential form $$\begin{align} &d*F=Z_0*J,~~ \tag {1} \end{align} $$ $$\begin{align} &dF=*\tilde{J},~~ \tag {2} \end{align} $$ where $F$ is the electromagnetic tensor, $J$ and $\tilde{J}$ are the electric current density and the magnetic current density, respectively, $Z_0$ is the vacuum wave impedance, $*$ is the Hodge star operator, and $d$ is the exterior differential operator. These equations can be written in the covariant derivative form $$\begin{align} &F^{\mu\nu}{}_{;\nu}=Z_0J^\mu,~~ \tag {3} \end{align} $$ $$\begin{align} &F_{\nu\sigma;\mu}+F_{\mu\nu;\sigma}+F_{\sigma\mu;\nu}=(*\tilde{J})_{\mu\nu\sigma},~~ \tag {4} \end{align} $$ and the partial derivative form $$\begin{align} \frac{1}{\sqrt{-g}}\frac{\partial}{\partial x^\nu}(\sqrt{-g}F^{\mu\nu})=\,&Z_0J^\mu,~~ \tag {5} \end{align} $$ $$\begin{align} \frac{\partial F_{\nu\sigma}}{\partial x^\mu}+\frac{\partial F_{\mu\nu}}{\partial x^\sigma}+\frac{\partial F_{\sigma\mu}}{\partial x^\nu}=\,&(*\tilde{J})_{\mu\nu\sigma},~~ \tag {6} \end{align} $$ where the indices in Greek letters run from 0 to 3, $F_{\mu\nu}$ is the electromagnetic tensor with the contravariant form $F^{\mu\nu}$, the semicolon in lower indices denote covariant derivative, and $g$ is the determinant of metric tensor $g_{\mu\nu}$. The metric tensor of Schwarzschild spacetime in Cartesian coordinate is[3,12] $g_{\mu\nu}={\rm diag}[-(\frac{1-\frac{R_{\rm S}}{4R}}{1+\frac{R_{\rm S}}{4R}})^2, (1+\frac{R_{\rm S}}{4R})^4,(1+\frac{R_{\rm S}}{4R})^4,(1+\frac{R_{\rm S}}{4R})^4]$, where $R_{\rm S}$ is the Schwarzschild radius and $R=\sqrt{(x^1)^2+(x^2)^2+(x^3)^2}$. Let us define $$\begin{align} {\boldsymbol D}/\epsilon_0=\,&\sqrt{-g}(F^{01},F^{02},F^{03}),~~ \tag {7} \end{align} $$ $$\begin{align} Z_0{\boldsymbol H}=\,&\sqrt{-g}(F^{23},F^{31},F^{12}),~~ \tag {8} \end{align} $$ $$\begin{align} {\boldsymbol E}=\,&(F_{10},F_{20},F_{30}),~~ \tag {9} \end{align} $$ $$\begin{align} c{\boldsymbol B}=\,&(F_{23},F_{31},F_{12}),~~ \tag {10} \end{align} $$ $$\begin{align} c\rho=\,&\sqrt{-g}J^0,~~ \tag {11} \end{align} $$ $$\begin{align} {\boldsymbol J}=\,&\sqrt{-g}(J^1,J^2,J^3),~~ \tag {12} \end{align} $$ $$\begin{align} c\tilde{\rho}=\,&\sqrt{-g}\tilde{J}^0,~~ \tag {13} \end{align} $$ $$\begin{align} \tilde{{\boldsymbol J}}=\,&\sqrt{-g}(\tilde{J}^1,\tilde{J}^2,\tilde{J}^3),~~ \tag {14} \end{align} $$ where $\epsilon_0$ is the vacuum permittivity, and $c$ is the vacuum light speed. Using this definition, the Maxwell Eqs. (5) and (6) can be written as the traditional 3D form with the constitutive equations $$\begin{align} {\boldsymbol D}=\,&\epsilon_0\Big(1-\frac{R_{\rm S}}{4R}\Big)^{-1}\Big(1+\frac{R_{\rm S}}{4R}\Big)^3{\boldsymbol E},~~ \tag {15} \end{align} $$ $$\begin{align} {\boldsymbol B}=\,&\mu_0\Big(1-\frac{R_{\rm S}}{4R}\Big)^{-1}\Big(1+\frac{R_{\rm S}}{4R}\Big)^3{\boldsymbol H},~~ \tag {16} \end{align} $$ where $\mu_0$ is the vacuum permeability. The Schwarzschild spacetime is equivalent to flat spacetime filled with certain medium, and it is relative permittivity and permeability are[12] $\epsilon_{\rm r}=\mu_{\rm r}=(1-\frac{R_{\rm S}}{4R})^{-1}(1+\frac{R_{\rm S}}{4R})^3$. The FVTD method for an inhomogeneous medium can be used directly in Schwarzschild spacetime. If $\tilde{J}=0$, then we have $dF=0$ according to Eq. (2). The electromagnetic tensor can be written as $$ F_{\alpha\beta}=c(dA)_{\alpha\beta}=c(A_{\beta;\alpha}-A_{\alpha;\beta}),~~ \tag {17} $$ where $A^\alpha$ is the electric potential, and $A_{\alpha;\beta}$ is the covariant derivative of $A_{\alpha}$ (covariant form of $A^\alpha$). By substituting the above equation into Eq. (3) and applying Lorentz gauge, we obtain the wave equation $$ \Box A^\alpha-R^\alpha{}_\beta A^\beta=-\mu_0J^\alpha ,~~ \tag {18} $$ where $\Box A^\alpha=A^\alpha{}_{;\beta}{}^{;\beta}$, the semicolon in upper indices denote contravariant derivative, and $R^\alpha{}_\beta$ is the Ricci tensor. The electric field and magnetic field are $$\begin{align} E_i=\,&F_{i0}=c(A_{0;1}-A_{1;0},A_{0;2}-A_{2;0},A_{0;3}-A_{3;0}),~~ \tag {19} \end{align} $$ $$\begin{align} Z_0H_i=\,&\sqrt{-g}\epsilon_{ijk}F^{jk}\\ =\,&c\sqrt{-g}(A^{3;2}-A^{2;3},A^{1;3}-A^{3;1},A^{2;1}-A^{1;2}),~~ \tag {20} \end{align} $$ where the indices $i,j,k$ run from 1 to 3, and $\epsilon_{ijk}$ is the Levi–Civita symbol. Let us define the dual electromagnetic tensor $\tilde{F}=*F$ (or $F=-*\tilde{F}$). Equations (1) and (2) can be rewritten as $$\begin{align} d*\tilde{F}=\,&-*\tilde{J},~~ \tag {21} \end{align} $$ $$\begin{align} d\tilde{F}=\,&Z_0*J.~~ \tag {22} \end{align} $$ Equations (3) and (4) can be rewritten as $$\begin{align} &\tilde{F}^{\mu\nu}{}_{;\nu}=-\tilde{J}^\mu,~~ \tag {23} \end{align} $$ $$\begin{align} &\tilde{F}_{\nu\sigma;\mu}+\tilde{F}_{\mu\nu;\sigma}+\tilde{F}_{\sigma\mu;\nu}=Z_0(*J)_{\mu\nu\sigma}.~~ \tag {24} \end{align} $$ If $J=0$, then we have $d\tilde{F}=0$ according to Eq. (22). The dual electromagnetic tensor can be written as $$ \tilde{F}_{\alpha\beta}=-\frac{1}{\epsilon_0}(d\tilde{A})_{\alpha\beta} =-\frac{1}{\epsilon_0}(\tilde{A}_{\beta;\alpha}-\tilde{A}_{\alpha;\beta}),~~ \tag {25} $$ where $\tilde{A}^\alpha$ is the magnetic potential. By substituting the above equation into Eq. (23) and applying Lorentz gauge, we obtain the wave equation $$ \Box \tilde{A}^\alpha-R^\alpha{}_\beta\tilde{A}^\beta=-\epsilon_0\tilde{J}^\alpha.~~ \tag {26} $$ The electric field and magnetic field are $$\begin{alignat}{1} \!\!\!\!\!\!\!\!\!\!E_i=\,&-(*\tilde{F})_{i0}=\sqrt{-g}(\tilde{F}^{23},\tilde{F}^{31},\tilde{F}^{12})\\ \!\!\!\!\!\!\!\!\!\!=\,&-\frac{\sqrt{-g}}{\epsilon_0}(\tilde{A}^{3;2}-\tilde{A}^{2;3},\\ &\tilde{A}^{1;3} -\tilde{A}^{3;1},\tilde{A}^{2;1}-\tilde{A}^{1;2}),~~ \tag {27} \end{alignat} $$ $$\begin{alignat}{1} \!\!\!\!\!\!\!\!\!\!Z_0H_i=\,&-\sqrt{-g}\epsilon_{ijk}(*\tilde{F})^{jk}=-(\tilde{F}_{10}, \tilde{F}_{20},\tilde{F}_{30})\\ \!\!\!\!\!\!\!\!\!\!=\,&\frac{1}{\epsilon_0}(\tilde{A}_{0;1}-\tilde{A}_{1;0},\tilde{A}_{0;2} -\tilde{A}_{2;0},\tilde{A}_{0;3}-\tilde{A}_{3;0}).~~ \tag {28} \end{alignat} $$ The total electric field is the summation of Eqs. (19) and (27), and the total magnetic field is the summation of Eqs. (20) and (28). We need to calculate the covariant derivative of potential $A_{\alpha;\gamma}$ and $\tilde{A}_{\alpha;\gamma}$. The solutions to wave equations (Eqs. (18) and (26)) are $$\begin{alignat}{1} \!\!\!\!\!\!A_\alpha(x)=\,&\frac{\mu_0}{4\pi}\int{G_{\alpha\beta'}(x,x')J^{\beta'}(x') \sqrt{-g(x')} d^4x'},~~ \tag {29} \end{alignat} $$ $$\begin{alignat}{1} \!\!\!\!\!\!\tilde{A}_\alpha(x)=\,&\frac{\epsilon_0}{4\pi}\int{G_{\alpha\beta'}(x,x') \tilde{J}^{\beta'}(x')\sqrt{-g(x')}d^4x'},~~ \tag {30} \end{alignat} $$ where $G_{\alpha\beta'}$ is the Green function with the Hadamard form[9] $$\begin{alignat}{1} G_{\alpha\beta'}(x,x')=\,&U_{\alpha\beta'}(x,x')\delta(\sigma(x,x'))\\ &+V_{\alpha\beta'}(x,x')\theta(-\sigma(x,x')),~~ \tag {31} \end{alignat} $$ where $U_{\alpha\beta'}$ and $V_{\alpha\beta'}$ are two bi-tensors, $\sigma$ is Synge's world function,[12] $\delta$ is the Dirac function, and $\theta$ is the step function $$ \theta(\sigma)=\begin{cases}\!\! 1,~~~~\sigma\geq 0,\\\!\! 0,~~~~\sigma < 0. \end{cases}~~ \tag {32} $$ The method to calculate $A_{\alpha;\gamma}$ using the Green function method has previously been introduced.[3] Compared Eqs. (29) and (30), the expression of $\tilde{A}_{\alpha;\gamma}$ can be obtained by replacing $\mu_0$ with $\epsilon_0$, $J$ with $\tilde{J}$. We use electric dipole as the excitation source.[4] The incident wave on connection boundary can be calculated using the Green function method. The current ($J^\alpha$, $\tilde{J}^\alpha$) and their time derivatives ($J^\alpha_{,0}$, $\tilde{J}^\alpha_{,0}$) on output boundary are calculated by Eqs. (7)-(10) in Ref. [4]. If we establish a local coordinate system on the output boundary and set the outer normal vector on output boundary as the $x^3$ axis, then Eq. (9) in Ref. [4] can be written as $$\begin{align} J^0_{,0}(k)=\,&-\frac{\partial J^1(k)}{\partial x^1}-\frac{\partial J^2(k)}{\partial x^2},~~ \tag {33} \end{align} $$ $$\begin{align} \tilde{J}^0_{,0}(k)=\,&-\frac{\partial\tilde{J}^1(k)}{\partial x^1}-\frac{\partial\tilde{J}^2(k)}{\partial x^2}.~~ \tag {34} \end{align} $$ By substituting Eq. (7) in Ref. [4] into the above equations, we obtain $$\begin{align} J^0_{,0}(k)=\,&\frac{\partial H_2(k)}{\partial x^1}-\frac{\partial H_1(k)}{\partial x^2},~~ \tag {35} \end{align} $$ $$\begin{align} \tilde{J}^0_{,0}(k)=\,&-\frac{\partial E_2(k)}{\partial x^1}+\frac{\partial E_1(k)}{\partial x^2}.~~ \tag {36} \end{align} $$ In FVTD, the gradient of electromagnetic field is computed by discretizing the equation $$ \iiint_V\nabla{\boldsymbol Q}d^3x={\LARGE ∯}_{\partial V}{\boldsymbol Q}{\boldsymbol n}d^2x ,~~ \tag {37} $$ where ${\boldsymbol n}$ is the unit outer normal vector on output boundary, and $$ {\boldsymbol Q}= \begin{bmatrix} {\boldsymbol E}\\ Z_0{\boldsymbol H} \end{bmatrix},~ \nabla{\boldsymbol Q}= \begin{bmatrix} \frac{\partial Q_1}{\partial x^1}&\frac{\partial Q_1}{\partial x^2}&\frac{\partial Q_1}{\partial x^3}\\ \vdots&\vdots&\vdots\\ \frac{\partial Q_6}{\partial x^1}&\frac{\partial Q_6}{\partial x^2}&\frac{\partial Q_6}{\partial x^3} \end{bmatrix}.~~ \tag {38} $$ By applying coordinate transform, we can obtain the quantities on the right side of Eqs. (35) and (36).
cpl-36-1-010401-fig1.png
Fig. 1. Validating the program and the connection boundary.
The first example is to validate the FVTD program. The Schwarzschild radius $R_{\rm S}$ is set to 1 m. There is no scatterer in FVTD domain. The PML region is a spherical shell with an inner radius of 0.3 m and an outer radius of 0.8 m. The connection boundary is a spherical surface with a radius of 0.15 m. The center of FVTD domain is (1.3 m, 0 m, 0 m). Figure 1 shows the configuration of the computational domain. The event horizon $R=R_{\rm S}/4$ is also shown to illustrate the position relative to event horizon. The FVTD domain is discretized by tetrahedron with the mean size of 0.05 m. A $z$-directed electric dipole is placed at (1 m, $-$1 m, 0 m) (the midpoint of the two charges). The waveform of charge is a Gaussian pulse $$ q(t)=A{\exp}\Big[-\frac{1}{2}\Big(\frac{t}{\sigma}-4\Big)^2\Big],~~ \tag {39} $$ where $A=1\times10^{-8}C$ and $\sigma=3.22$ ns. The $z$-components of the electric field at two positions are calculated by the FVTD method. These positions are Pt(1.3 m, 0 m, 0 m) and Ps(1.5 m, 0 m, 0 m). Pt is located at total field zone and Ps is located at scatter field zone. The electric fields at Pt are also calculated by GFM and FDTD. Figure 2 shows the result at Pt by the three methods. This demonstrates that the numerical results obtained with all the methods match perfectly. This example validates the FVTD program and it also validates the connection boundary. The results at Ps are shown in Fig. 3, in which the values of electric field are in dB: $20{\lg}(|E_z|/{\max}|E_{\rm zt}|)$, where $E_{\rm zt}$ is the electric field at Pt. This demonstrates that the numerical scatter fields by FVTD are less than $-$60 dB.
cpl-36-1-010401-fig2.png
Fig. 2. The values of $E_z$ at Pt.
cpl-36-1-010401-fig3.png
Fig. 3. The values of $E_z$ at Ps.
cpl-36-1-010401-fig4.png
Fig. 4. Scattered by a disk.
The second example is scattered by a disk. The radius of the perfectly electric conductor (PEC) disk is 0.5 m and the thickness is 0.05 m. A $z$-directed electric dipole is placed at (6 m, 0 m, 0 m) (Fig. 4). The center of the disk locates at (3.5 m, 0 m, 0 m) which is also the center of the FVTD domain. The Schwarzschild radius $R_{\rm S}$ is set to 1 m. The connection boundary is an ellipsoid surface with a semi-axis of (0.1 m, 0.6 m, 0.6 m). The output boundary is an ellipsoid surface with a semi-axis of (0.2 m, 0.7 m, 0.7 m). The PML region is a spherical shell with an inner radius of 0.9 m and an outer radius of 1.4 m (Fig. 5). The FVTD domain is discretized by tetrahedron with the mean size of 0.05 m. The waveform of charge is a Gaussian pulse with $A=1\times10^{-5}C$ and $\sigma=1.93$ ns (Eq. (39)).
cpl-36-1-010401-fig5.png
Fig. 5. Validating the output boundary.
cpl-36-1-010401-fig6.png
Fig. 6. The values of $E_z$ at P1.
cpl-36-1-010401-fig7.png
Fig. 7. The values of $E_z$ at P2 in time domain.
The $z$ components of the electric field at P1(4.3 m, 0 m, 0 m) are calculated by the FVTD method. The electric fields at this position are also calculated by integrating the output boundaries using Green function method. The results are shown in Fig. 6, which validate the output boundary. The $z$-component of the scattered electric field (in both time domain and frequency domain) at P2(7.5 m, 0 m, 0 m) is shown in Figs. 7 and 8. The scattered electric field in flat spacetime ($R_{\rm S}=0$) is also shown in these two figures. From Eqs. (15) and (16), we can see that the effective light speed is smaller than that in flat spacetime. This leads to a time delay, as shown in Fig. 7. The inhomogeneity leads to pulse broadening in time domain and red shift in frequency domain (Fig. 8).
cpl-36-1-010401-fig8.png
Fig. 8. The values of $E_z$ at P2 in frequency domain.
In summary, the FVTD method in curved spacetime has been realized by filling flat spacetime with equivalent medium. Green function in curved spacetime is calculated by solving differential equations. These two methods are incorporated to simulating electromagnetic scattering in Schwarzschild spacetime. We validate the FVTD program and the Green function program by two numerical examples. The scattering field by a disk is simulated by the developed methods. In the numerical simulation, the Green functions linking source point to every surface element on the connection boundary should be computed, and the Green functions linking field point to every surface elements on the output boundary should also be computed. The computation is very time consuming. Developing the fast algorithm will be our future work. References Electromagnetic waves in a strong Schwarzschild plasmaA method for incorporating the Kerr–Schild metric in electromagnetic particle-in-cell codeNumerical simulation of electromagnetic waves in Schwarzschild space–time by finite difference time domain method and Green function methodNew version of FDTD with GFM program in Schwarzschild space–timeSpherical Perfectly Matched Absorber for Finite-Volume 3-D Domain TruncationFinite Volume Time Domain modelling of microwave breakdown and plasma formation in a metallic apertureThe Motion of Point Particles in Curved SpacetimeTransport equation approach to calculations of Hadamard Green functions and non-coincident DeWitt coefficients
[1] Daniel J and Tajima T 1997 Phys. Rev. D 55 5193
[2] Watson M and Nishikawa K I 2010 Comput. Phys. Commun. 181 1750
[3] Jia S Q, La D S and Ma X L 2018 Comput. Phys. Commun. 225 166
[4] Jia S Q 2018 Comput. Phys. Commun. 232 264
[5]Rao S M 1999 Time Domain Electromagnetics (New York: Academic Press)
[6] Fumeaux C, Sankaran K et al 2007 IEEE Trans. Microwave Theory Tech. 55 2773
[7] Hamiaz A, Klein R et al 2012 Comput. Phys. Commun. 183 1634
[8]Friedlander F 1975 The Wave Equation on a Curved Space-Time (Cambridge: Cambridge University Press)
[9] Poisson E 2004 Living Rev. Relativ. 7 6
[10] Ottewill A C and Wardell B 2011 Phys. Rev. D 84 104039
[11]Wardell B 2012 PhD Dissertation (Dublin: University College Dublin)
[12]Synge J L 1960 Relativity: the General Theory (Amsterdam: North-Holland Publishing)