Open Access
Add to BibSonomy
Abstract
Electromagnetic scattering in moving structures is a fundamental topic in physics and engineering. Yet no general numerical solution to related problems has been reported to date. We introduce here a generalized FDTD scheme to remedy this deficiency. That scheme is an extension of the FDTD standard Yee cell and stencil that includes not only the usual, physical fields but also auxiliary, unphysical fields allowing a straightforward application of moving boundary conditions. The proposed scheme is illustrated by four examples – a moving interface, a moving slab, a moving crystal and a moving gradient – with systematic validation against exact solutions.
© 2023 Optica Publishing Group under the terms of the Optica Open Access Publishing Agreement
1. Introduction
Moving structures are pervasive in electromagnetics. They may involve both moving bodies, such as stars or vehicles, and moving perturbations, such as fluid or elastic waves. The category of moving-body structures has a nearly 300-year history and has led to the discovery of many fundamental physical phenomena, such as Bradley aberration [1], Doppler frequency shifting [2], Fizeau dragging [3], relativity [4, 5], magneto-electric coupling [6], medium bianisotropy [7] and gravity emulation [8]. The category of moving-perturbation structures is more recent and particularly amenable to practical applications; it produces effects such as parametric amplification [9–11], nonreciprocal transmission [12,13], deflected reflection [14,15], space-time frequency transitions [16,17], pseudo-Fizeau dragging [18,19], spatiotemporal bandgaps [20–22] and acceleration-based beam bending [23].
Surprisingly, no general numerical tool is currently available to simulate electromagnetic moving structures, whether of the moving-body or moving-perturbation type. The development of such a tool would be highly desirable to simulate the plethora of moving structures mentioned above as well as future ones, particularly the non-canonical – and typically more practical – structures that do not admit analytical solutions. This issue is well-known. There have been only two workarounds in the literature so far [24–28], both based on the FDTD technique [29,30], probably selected for its natural incarnation of both spatial and temporal variations in Maxwell’s equations. However, one of these approaches is restricted to non-penetrable objects [24,25], while the other one implies cumbersome Lorentz frame transformations [26–28].
We present here a novel, general, and efficient FDTD scheme for simulating electromagnetic scattering in moving structures. That scheme, contrarily to that in [24,25], also applies to penetrable media, allowing hence to handle gradient structures and metamaterials [18,31], without requiring problematic numerical transformations, contrarily to that in [26–28]. This scheme is based on a generalized Yee cell and stencil, which include not only the usual, physical fields but also auxiliary, unphysical fields that carry the velocity information, to automatically satisfy the moving boundary conditions.
2. Types of systems
The proposed method applies to both moving-matter and moving-perturbation media [19], the only difference between the two, in terms of modeling, being that in the former case, motion transforms media that are isotropic at rest into bianisotropic media [32,33], whereas in the latter case, motion alters the values of the permittivity or permeability parameters without affecting the isotropy of the medium [16,34,35]. Figure 1 depicts some of these structures [50,53], which will be used for illustration and validation purposes in Sec. 6: a simple interface [36] [Fig. 1(a)], a slab [37] [Fig. 1(b)], a bilayer crystal [Fig. 1(c)], which can be operated in the Bragg regime [10] or in the metamaterial (homogeneous) regime [18], and a gradient [38] [Fig. 1(d)], all moving at a velocity $v$ in a direction selected as being $z$ in a Cartesian coordinate system. The different parts of these structures can represent arbitrary bianisotropic media, with general tensor $\overline {\overline {\chi }}=\left [\overline {\overline {\epsilon }},\overline {\overline {\xi }};\overline {\overline {\zeta }},\overline {\overline {\mu }}\right ]$, which may be all different from each other and may be arbitrarily truncated in the $xy$-plane. Note that all these structures can be seen as a succession of interfaces, including the gradient upon decomposition into subwavelength slabs. For this reason, we shall focus on an interface in the next section.
Fig. 1. Types of moving structures, of either the moving-matter or moving-perturbation type, that can be simulated by the proposed generalized FDTD scheme. (a) Simple interface. (b) Slab. (c) Bilayer crystal. (d) Gradient.
3. Generalized Yee cell
The boundary conditions at a moving interface, whether associated with moving modulation or moving matter, read [39,40]
withThe fields in Eq. (2) may be decomposed into tangential and normal components with respect to the interface. The tangential components are
In the case where the interface is perpendicular to the direction of motion, i.e., $\mathbf {n}\|\mathbf {v}$, the right-most terms in the right-hand sides of Eqs. (3a) and (3b) vanish. This is the scenario that is represented in Fig. 1 and that will be assumed in the sequel of the paper [51,52,54]. The paper will further assume non-conducting interfaces, and hence $\mathbf {J}_\text {s}=\rho _\text {s}=0$ in Eq. (1).
The dynamic boundary conditions (1) differ from the static or stationary boundary equations only insofar as they involve the auxiliary fields $\mathbf {E}_{1,2}^*$ and $\mathbf {H}_{1,2}^*$, whose tangential components include velocity-dependent extra terms [Eq. (3a)], instead of the fields $\mathbf {E}_{1,2}$ and $\mathbf {H}_{1,2}$. This observation suggests the generalization of the standard Yee cell [37] to a Yee cell with correspondingly modified $\mathbf {E}_\text {tan}$ and $\mathbf {H}_\text {tan}$ fields, without alteration of the conventional electric and magnetic field staggering arrangement. The result is shown in Fig. 2, for a pair of generalized Yee cells separated by a moving interface. The so-defined generalized Yee cell straightforwardly allows the field continuity conditions at a moving interface, i.e., Eqs. (1a) and (1b), to be satisfied by having the field components $E_x^*$, $E_y^*$ and $H_z$ uniquely defined at the permittivity interface of Fig. 2 (or $H_x^*$, $H_y^*$ and $E_z$ in the case of a corresponding magnetic interface).
Fig. 2. Pair of generalized Yee cells, separated by a moving interface [as in Fig. 1(a)], with the generalization consisting in the substitution of the tangential physical fields $\mathbf {E}_\text {tan}$ and $\mathbf {H}_\text {tan}$ by the corresponding tangential auxiliary fields $\mathbf {E}^*_\text {tan}$ and $\mathbf {H}^*_\text {tan}$ in Eq. (3a).
This simple generalization of the Yee cell is the core idea of the paper. It implies a mixture of physical and unphysical fields, namely, as seen in Eqs. (1) and (3), the physical fields $\mathbf {D}$, $\mathbf {B}$, $\mathbf {E}_{\mathrm{norm}}$ and $\mathbf {H}_{\mathrm{norm}}$ and the unphysical fields $\mathbf {E}_{\mathrm{tan}}^*$ and $\mathbf {H}_{\mathrm{tan}}^*$, which will also imply unphysical mixed-field electromagnetic equations. This strange situation is somewhat reminiscent of the unphysical anisotropic media appearing in Bérenger’s Perfect Matched Layers (PMLs) [41] and, as for the PMLs, it does not pose any fundamental problem. The corresponding unphysical Maxwell’s equations and constitutive relations, which will be derived and discretized in Sec. 4 to establish the generalized FDTD update equations, are mathematical mappings of the physical Maxwell’s equations via the simple changes of variable in Eqs. (3), and the physical fields $\mathbf {E}$ and $\mathbf {H}$ can be found anytime from the auxiliary fields by simply inverting these relations.
4. Generalized scheme
In this section, we shall derive the generalized FDTD scheme corresponding to the generalized Yee cell established in Sec. 3. This scheme must, of course, model the physics of the problem, and hence follow the usual Maxwell-Faraday and Maxwell-Ampère equations,
In order to apply the generalized Yee cell established in Sec. 3, we must express the fields $\mathbf {E}$ and $\mathbf {H}$ in Eqs. (4) and (5) in terms of the fields $\mathbf {E}^*$ and $\mathbf {H}^*$ using Eq. (2). This results in the mixed-field Maxwell’s equations
From this point, we shall restrict our developments, for the sake of simplicity and without any loss of generality, to structures that are i) based on moving modulation (as opposed to moving matter), ii) one-dimensional (1D) and iii) illuminated by a plane wave propagating in the direction of the modulation (i.e., normal to the interface(s)). Using a Cartesian coordinate system with $z$ corresponding to the direction of the modulation, the non-zero field components reduce to $E_x$ and $H_y$ and Eqs. (6) and (7) become then
We can now proceed to the discretization of Eqs. (8) and (9). The static ($v=0$) parts can be discretized in the usual fashion [30], which yields
In contrast, the dynamic ($v\neq 0$) terms, which have not been discretized yet in Eqs. (10), require special treatment for the numerical stability and minimal dispersion of the overall scheme. We empirically [57] found that the discretization choices given in Table 1 are appropriate in these regards. Note that these schemes are different for positive velocities, i.e., $\mathbf {v}\|\mathbf {\hat {z}}$ or $v>0$ and negative velocities, $\mathbf {v}\|-\mathbf {\hat {z}}$ or $v<0$.
Table 1. Discretization of the dynamic terms in Eqs. (10).
Figure 3 shows the stencil of the overall proposed generalized FDTD scheme, which provides a handy pictorial representation of Eqs. (10) with Eqs. (11)–(14).
5. Stability
In this section, we demonstrate the stability of the generalized scheme presented in Sec. 4 using von Neumann’s approach [43]. This approach consists in inserting plane-wave test fields of the form
where $\xi$ and $\zeta$ are the space-dependent oscillatory and time-dependent magnitude-oscillatory functions with $\alpha$ being the attenuation/amplification rate and oscillation, into the update equations and determining the condition under which the $\zeta (\xi )$ solution of the resulting matrix system is smaller than one – but as close as possible to one to avoid numerical attenuation [58].Inserting the test fields of the above form into Eqs. (10) with Eqs. (11a)–(14)(a) or (11)(b)–(14)(b), setting the determinant of the matrix associated with the resulting homogeneous system of equations to zero, and finding the $\zeta$ roots of the corresponding characteristic polynomial yields (see details in Appendix A.)
The fact that the stability analysis admits a double solution, $\zeta _\pm$ [Eq. (16)], may a priori appear surprising, given that $\zeta$ solutions are usually unique. However, this can be understood by realizing that waves that are co-moving and contra-moving with respect to the modulation see different numerical “landscapes”. Specifically, if $v>0$, then the scheme uses Eq. (11a)-(14a) [first column in Table 1] for the dynamic parts, which correspond to the dotted segments (with the dashed segment removed) in Fig. 3; the resulting stencils for the first three equations [Eqs. (10a–10c) with Eqs. (11a)–(13a)] are clearly asymmetric with respect the $z$ direction and therefore, although the stencil for the fourth equation [Eq. (10d) with Eq. (14a)] is symmetric, the stencil is overall asymmetric with respect to $z$, so that co-moving and contra-moving waves will see different stencils. A similar conclusion holds for the case $v<0$, so we expect that, for a given choice of the modulation direction, $\zeta _+$ will correspond to one of the two (co-moving or contra-moving) regimes, while $\zeta _-$ will correspond to the other regime. Note that the scheme must be stable for both waves since typical problems involve both of them [59]. The stability criterion ($\zeta <1$) should then be taken as corresponding to the more restrictive (larger $\zeta$) of the two $\zeta$ solutions.
The two regimes may be identified by comparing the signs of the spatial and temporal phases of the test plane-wave in Eq. (15a), with the spatial phase and temporal phase residing in $\xi$ and $\zeta _\pm$, respectively. For this purpose, let us consider a specific numerical position and time, say for simplicity $k=1$ and $n=1$, and analyze the test plane wave for the two possible propagation directions. We have, upon writing $\alpha _\pm =\omega _{\text {i}\pm }\pm i\omega _{\text {r}\pm }$ in $\zeta _\pm$ (with $\omega _{\text {i}\pm }>0$ assuming stability),
In these relations, the first expression of the second line reveals that the upper sign, which is associated with $\zeta _+$, corresponds to a forward wave, propagating in the $+z$ direction, while the lower sign, which is associated with $\zeta _-$, corresponds to a backward wave, so that $\zeta _+$ and $\zeta _-$ respectively correspond to the co-moving and contra-moving regimes for $v>0$ (and conversely for $v<0$).
Unfortunately, it is not possible to generally assert whether the $\pm$ signs in Eq. (17) correspond to the $\pm$ signs in Eq. (16), or their opposite, because Eq. (16) does not admit a general decomposition in real and imaginary parts. Therefore, that sign correspondence has to be determined in a numerical fashion. For instance, for the set of parameters $S=0.5$, $\beta =0.3$, $\epsilon =4$, $\mu =1$ and $\Delta z k_z=2\pi /5$, we have $\zeta _+=0.925-i 0.33$ and $\zeta _-=0.917+i 0.23$, and therefore the $\pm$ signs in Eq. (17) correspond to $\pm$ signs in Eq. (16). The correct sign mapping has to be determined on a case-by-case basis for different simulation parameters.
Let us now illustrate and verify our stability analysis. Figure 4 plots the attenuation rate, $\alpha =\ln (\zeta )/\Delta t$, versus the meshing density, $N_\lambda =\lambda _0/\Delta z=2\pi /(k_z\Delta z)=i 2\pi /\ln (\xi )$, for Gaussian pulse propagation and for the case $v>0$ [corresponding to the stencils with Eqs. (11a)–(14a)]. The analytical results are obtained by Eq. (16) while the numerical results are obtained by FDTD simulations and taking ratios of the pulse maxima at different times. The insets show the FDTD-simulated space-time field evolutions for two meshing densities. The FDTD attenuation results closely match the attenuation results predicted by the analytical stability analysis, with distinct levels corresponding to the co-moving and contra-moving regimes. Note that in the present case ($v>0$ and chosen parameters) and throughout the considered $N_\lambda$ range ($10<N_\lambda <100$), $\zeta ^\pm$ in Eq. (16) corresponds to $\zeta ^\pm$ in Eq. (17), so that $\zeta _+$ and $\zeta _-$ correspond to co-moving and contra-moving waves, respectively, as in the above numerical example. We see from the plot that the appropriate (more restrictive) criterion corresponds to that for the co-moving wave ($\zeta _+$, lower attenuation), and its satisfaction will a fortiori ensure the stability of the contra-moving wave ($\zeta _-$, larger attenuation). Symmetrically identical results (not shown) are obtained for the case $v>0$, corresponding to the stencil with Eqs. (11b)–(14b).
Fig. 4. Attenuation rate versus meshing density obtained by Eq. (16) and by FDTD simulation, for $S=0.5$, $v=0.3c$ ($v>0$, first column in Table 1), $\epsilon =1$ and $\mu =1$. The insets show FDTD-simulated space-time evolution of a $3\lambda$-wide Gaussian pulse modulated by a harmonic wave of wavelength $\lambda =2\pi /k_z$ and period $T=2\pi /\omega$, over a distance of $20\lambda$ and time of $45T$, for the mesh densities $N_\lambda =10$ and $N_\lambda =80$.
6. Illustrative and validating examples
Figures 5, 6, 7 and 8 provide illustrative and validating examples of scattering from the structures corresponding to Figs. 1(a), (b), (c) and (d): a moving interface, a moving slab, a moving crystal, and a moving gradient. In all the cases, the incident wave is a modulated Gaussian pulse with (starred) electric field $E^*_\text {i}(t)=e^{-i\omega _\text {i} t}e^{-(t/\tau _\text {i})^2}$, and the computational domain is terminated by standard Mur boundary conditions, to which the generalized Yee cells are automatically matched, as shown in App. B. The panels (a) show the space-time evolution of the incident and scattered fields, with the input and output temporal waveforms at specific $z$ positions plotted on the sides, for a pulse with $\tau _i=T_i=2\pi /\omega _i$. The panels (b) plot the Fourier transforms of these waveforms. The panels (c) plot the scattering coefficient functions $\Gamma (\omega )$ and $T(\omega )$, obtained for an incident pulse $E_i$ that is temporally short enough to have a nonzero spectrum across the frequency range of interest, as done for instance in [44]. All the analytical solutions are given in Appendix C.
Fig. 5. Results for a contra-moving interface [Fig. 1(a)], for the physical parameters $v=-0.3c$, $\epsilon _1=1$, $\epsilon _2=4$ and $\mu =1$, and numerical parameters $S=0.2$ and $\Delta z=\lambda _\text {i}/150$.
Fig. 6. Same as in Fig. 5 but for a co-moving slab [Fig. 1(b)], with velocity $v=+0.3c$ and length $\ell =\frac {\lambda _\text {i}}{4n_2}\frac {1+n_2\beta }{1-n_1\beta }$, with $\lambda _\text {i}=2\pi c/(n_1\omega _\text {i})$ corresponding to a space-time quarter-wave slab [30].
Fig. 7. Same as in Fig. 5 but for a co-moving bilayer crystal [Fig. 1(c)], consisting of 5 unit-cells each made of two space-time quarter-wave slabs, with lengths $\ell _1=\frac {\lambda _\text {i}}{4n_1}\frac {1+n_1\beta }{1-n_2\beta }$ and $\ell _2=\frac {\lambda _\text {i}}{4n_2}\frac {1+n_2\beta }{1-n_1\beta }$.
Fig. 8. Same as in Fig. 5 but for a co-moving gradient [Fig. 1(d)], corresponding to a linear increase of the permittivity from $\epsilon _1$ to $\epsilon _2$ over a length of $\ell =\lambda _\text {i}/2$.
In all the cases (Figs. 5, 6, 7 and 8) and for all the quantities (incident, reflected and transmitted electric fields, and reflection and transmission coefficients), excellent agreement is observed between the numerical results of the proposed generalized FDTD scheme and the analytical results (given in Appendix C.), with minor observable discrepancies in the scattering coefficients (panels (c)), which are explained by the limited bandwidth of the test pulse.
The physics of the dynamic scattering observed in Figs. 5, 6, 7 and 8 involves reflection Doppler frequency shifting, transmission index contrast frequency shifting, co-moving attenuation, contra-moving amplification, multiple space-time scattering, and space-time stopbands. All these effects have been described elsewhere [31,35,45] and are therefore not discussed here.
7. Conclusion
We have presented a generalized FDTD scheme that can simulate electromagnetic scattering in a great diversity of moving-medium problems. This scheme, using a combination of physical and auxiliary fields to satisfy moving boundary conditions in a natural fashion, is both simple and powerful. It fills a fundamental and important gap in the toolbox of electromagnetic computational techniques and is hence expected to find wide applications, particularly in the emergent area of space-time metamaterials.
Appendices
A. Stability calculation
This section derives the stability condition for the generalized FDTD scheme, presented in Sec. 4, using von Neumann’s approach [43]. We will therefore replace everywhere the electromagnetic fields with their plane-wave test counterpart, $\Psi |_k^n=\Psi _0\zeta ^n e^{ik_z k\Delta z}$, with $\zeta =e^{-\alpha \Delta t}$. In addition, we will use the convenient parameter $S=\Delta t/\Delta z$ wherever appropriate.
We shall restrict our derivations to the positive-velocity regime ($v>0$), whose dynamic part corresponds to the first column of Table 1. The negative-velocity derivations ($v<0$), whose dynamic part corresponds to the second column of Table 1, can easily be obtained upon following similar steps.
Using the plane-wave test field substitution in Eqs. (10a) and (11a), and dividing the resulting expression by $\zeta ^n$ and $e^{ik_z k\Delta z}$ yields
Similarly, Eq. (10b) with Eq. (12a) becomes
Next, Eq. (10c) with (13a) becomes, after factoring out $e^{-ik_z\Delta z}$ from the last two terms to form a cosine,
Finally, Eq. (10d) with (14a) becomes, after factoring out $e^{-ik_z\Delta z/2}$ from the last terms to form a cosine,
Equations (19), (21), (23) and (25) form a system of equations that can be put in the matrix form
Setting the determinant of this homogeneous system to zero to find a nontrivial solution for $\zeta$ leads to the characteristic polynomial
This is a polynomial of the second order in $\zeta$ whose solution is given by Eq. (16). The corresponding polynomial and solutions turn out to be exactly identical for the negative-velocity regime, although the related derivation involves distinct equations.
B. Impedance matching between the standard and the modified schemes
Inserting the constitutive relation (9b) into (8a), we find
Assuming medium homogeneity [$\epsilon \neq \epsilon (z,t)$, $\mu \neq \mu (z,t)$], as at the boundaries of the computational domain, these relations simplify to
Inserting (8b) into this result, we obtain
For a monochromatic plane wave, this relation reduces to
Thus, the starred fields have the same impedance as their non-starred counterparts, and the former are therefore matched to the latter.
C. Analytical solutions
This section provides analytical solutions that are used for validation of the illustrative examples given in Sec. 6. The solutions in Secs. C.2. and C.3 are closed-form solutions for the canonical problems of an interface and a slab, corresponding to Figs. 1(a) and 1(b), respectively, while the solutions in Sec. C.4 are transfer-matrix method solutions for the more complex problems of a bilayer crystal and a gradient, corresponding to Figs. 1(c) 1(d), respectively.
C.1. General formulas
In all cases, the incident field is a modulated Gaussian pulse, with temporal and spectral profiles
Although scattering from modulated interfaces has been extensively studied in the literature (e.g. [34,35,46,47]), the formula (35c) has not been explicitly published anywhere, to the best of our knowledge. The reason is that the quasi-totality of studies have been restricted to monochromatic waves, whose reflected field spectrum particularizes to $a_\text {r,t}\delta (\omega -a_\text {r,t}\omega _\text {i})\tilde A^*_\text {r,t}=\delta (\omega /a_\text {r,t}-\omega _\text {i})\tilde{A}^*_\text {r,t}$, after using the identity $a\delta (at)=\delta (t)$, where the $a_\text {r,t}$ factor has disappeared. However, the factor $a_\text {r,t}$ in Eq. (35c) cannot be dropped for other wave regimes such as the pulse regime.
C.2. Interface
The scattering coefficients for a modulated interface are found by applying the continuity conditions (1a) and (1b). They read [31,34,45], assuming wave incidence from medium 1 towards medium 2,
C.3. Slab
The scattering coefficients for a modulated slab may be found, as for a stationary slab, either by the multiple-reflection method or the boundary-value method [48]. This results in the scattering coefficients [31]
These relations are written in terms of $\omega _\text {i}$ through the expressions of the frequencies for waves propagating forward ($\omega _2^+$) and backward ($\omega _2^-$), namely
C.4. Transfer-matrix method
The transfer-matrix method [20,31,34] allows to compute electromagnetic scattering for any traveling-wave modulation profile with uniform (space-time) slabs in the frequency domain [60]. This is even true for a smooth gradient modulation profile, which can be approximated as a succession of subwavelength interfaces with constant parameters.
The transfer-matrix method models the interface scattering by the interface matrix
The transfer matrix of a structure composed of $N$ slabs is then found by taking the product of alternating discontinuity and propagation matrices, viz.,
Finally, the scattering coefficients can be obtained by converting the transfer matrix into a scattering matrix [49], which yields
where $A$, $B$, $C$ and $D$ are the $(1,1)$, $(1,2)$, $(2,1)$ and $(2,2)$ components of the $[M]_{N1}$ matrix.Disclosures
The authors declare no conflicts of interest.
Data availability
Data underlying the results presented in this paper are not publicly available at this time but may be obtained from the authors upon request.
References
1. J. Bradley, “IV. A letter from the Reverend Mr. James Bradley Savilian Professor of Astronomy at Oxford, and FRS to Dr. Edmond Halley Astronom. Reg. &c. giving an account of a new discovered motion of the fix’d stars,” Phil. Trans. R. Soc. 35(406), 637–661 (1728). [CrossRef]
2. C. Doppler, Über das farbige Licht der Doppelsterne und einiger anderer Gestirne des Himmels (Calve, 1842).
3. H. Fizeau, “Sur les hypothèses relatives à l’éther lumineux, et sur une expérience qui parait démontrer que le mouvement des corps change la vitesse avec laquelle la lumière se propage dans leur intérieur,” C.R. Acad. Sci. 33, 349 (1851).
4. A. Einstein, “Zur elektrodynamik bewegter Körper,” Ann. Phys. 322(10), 891–921 (1905). Reprinted in [49]. [CrossRef]
5. H. A. Lorentz, A. Einstein, H. Minkowski, H. Weyl, and A. Sommerfeld, The Principle of Relativity: a Collection of Original Memoirs on the Special and General Theory of Relativity (Dover Publications, 1952).
6. W. C. Röntgen, “Ueber die durch Bewegung eines im homogenen electrischen Felde befindlichen Dielectricums hervorgerufene electrodynamische Kraft,” Ann. Phys. 271(10), 264–270 (1888). [CrossRef]
7. H. Minkowski, “Die grundgleichungen für die elektromagnetischen vorgänge in bewegten körpern,” Nachrichten von der Gesellschaft der Wissenschaften zu Göttingen, Mathematisch-Physikalische Klasse 1908, 53–111 (1908).
8. U. Leonhardt and P. Piwnicki, “Optics of nonuniformly moving media,” Phys. Rev. A 60(6), 4301–4312 (1999). [CrossRef]
9. P. K. Tien, “Parametric amplification and frequency mixing in propagating circuits,” J. Appl. Phys. 29(9), 1347–1357 (1958). [CrossRef]
10. E. S. Cassedy and A. A. Oliner, “Dispersion relations in time-space periodic media: Part I-Stable interactions,” Proc. IEEE 51(10), 1342–1359 (1963). [CrossRef]
11. E. S. Cassedy, “Dispersion relations in time-space periodic media: Part II-Unstable interactions,” Proc. IEEE 55(7), 1154–1168 (1967). [CrossRef]
12. Z. Yu and S. Fan, “Complete optical isolation created by indirect interband photonic transitions,” Nat. Photonics 3(2), 91–94 (2009). [CrossRef]
13. N. Chamanara, S. Taravati, Z.-L. Deck-Léger, and C. Caloz, “Optical isolation based on space-time engineered asymmetric photonic band gaps,” Phys. Rev. B 96(15), 155409 (2017). [CrossRef]
14. Z.-L. Deck-Léger, A. Akbarzadeh, and C. Caloz, “Wave deflection and shifted refocusing in a medium modulated by a superluminal rectangular pulse,” Phys. Rev. B 97(10), 104305 (2018). [CrossRef]
15. A. M. Shaltout, V. M. Shalaev, and M. L. Brongersma, “Spatiotemporal light control with active metasurfaces,” Science 364(6441), 1–10 (2019). [CrossRef]
16. M. Lampe, E. Ott, and J. H. Walker, “Interaction of electromagnetic waves with a moving ionization front,” Phys. Fluids 21(1), 42–54 (1978). [CrossRef]
17. N. Chamanara, Y. Vahabzadeh, and C. Caloz, “Simultaneous control of the spatial and temporal spectra of light with space-time varying metasurfaces,” IEEE Trans. Antennas Propag. 67(4), 2430–2441 (2019). [CrossRef]
18. P. A. Huidobro, E. Galiffi, S. Guenneau, R. V. Craster, and J. B. Pendry, “Fresnel drag in space-time modulated metamaterials,” Proc. Natl. Acad. Sci. U.S.A. 116(50), 24943–24948 (2019). [CrossRef]
19. C. Caloz, Z.-L. Deck-Léger, A. Bahrami, O. C. Vincente, and Z. Li, “Generalized space-time engineered modulation (GSTEM) metamaterials: A global and extended perspective,” IEEE Antennas Propag. Mag. (2022).
20. F. Biancalana, A. Amann, A. V. Uskov, and E. P. O’Reilly, “Dynamics of light propagation in spatiotemporal dielectric structures,” Phys. Rev. E 75(4), 046607 (2007). [CrossRef]
21. O. Mattei and G. W. Milton, “Field patterns without blow up,” New J. Phys. 19(9), 093022 (2017). [CrossRef]
22. Y. Sharabi, A. Dikopoltsev, E. Lustig, Y. Lumer, and M. Segev, “Spatiotemporal photonic crystals,” Optica 9(6), 585–592 (2022). [CrossRef]
23. A. Bahrami, Z.-L. Deck-Léger, and C. Caloz, “Electrodynamics of metamaterials formed by accelerated modulation,” arXiv, arXiv:2208.11506 (2022). [CrossRef]
24. F. Harfoush, A. Taflove, and G. A. Kriegsmann, “Numerical implementation of relativistic electromagnetic boundary conditions in a laboratory-frame grid,” J. Comput. Phys. 89(1), 80–94 (1990). [CrossRef]
25. F. Harfoush, A. Taflove, and G. A. Kriegsmann, “A numerical technique for analyzing electromagnetic wave scattering from moving surfaces in one and two dimensions,” IEEE Trans. Antennas Propagat. 37(1), 55–63 (1989). [CrossRef]
26. H. Iwamatsu and M. Kuroda, “Over set grid generation method coupled with FDTD method while considering the Doppler effect,” IEEJ Trans. Fundam. Mater. 129(10), 699–703 (2009). [CrossRef]
27. K. Zheng, Z. Mu, H. Luo, and G. Wei, “Electromagnetic properties from moving dielectric in high speed with Lorentz-FDTD,” Antennas Wirel. Propag. Lett. 15, 934–937 (2016). [CrossRef]
28. Y. Zhao and S. Chaimool, “Relativistic finite-difference time-domain analysis of high-speed moving metamaterials,” Sci. Rep. 8(1), 7686 (2018). [CrossRef]
29. W. F. Ames, Numerical Methods for Partial Differential Equations (Academic Press, 1977).
30. A. Taflove and S. C. Hagness, Computational Electrodynamics: The Finite-Domain Time-Domain Method (Artech House, 2000), 2nd ed.
31. Z.-L. Deck-Léger, N. Chamanara, M. Skorobogatiy, M. Silveirinha, and C. Caloz, “Uniform-velocity spacetime crystals,” Adv. Photonics 1(5), 056002 (2019). [CrossRef]
32. D. Cheng and J. A. Kong, “Time-harmonic fields in source-free bianisotropic media,” J. Appl. Phys 39(12), 5792–5796 (1968). [CrossRef]
33. B. M. Bolotovski and S. N. Stolyarov, “Current status of the electrodynamics of moving media (infinite media),” Sov. Phys. Usp. 17(6), 875–895 (1975). [CrossRef]
34. C. S. Tsai and B. A. Auld, “Wave interactions with moving boundaries,” J. Appl. Phys. 38(5), 2106–2115 (1967). [CrossRef]
35. C. Caloz and Z.-L. Deck-Léger, “Spacetime metamaterials, Part I: General concepts,” IEEE Trans. Antennas Propag. 68(3), 1569–1582 (2020). [CrossRef]
36. Z. Li, X. Ma, A. Bahrami, Z.-L. Deck-Léger, and C. Caloz, “Generalized total internal reflection at dynamic interfaces,” Phys. Rev. B 107(11), 115129 (2023). [CrossRef]
37. K. Yee, “Numerical solution of initial boundary value problems involving Maxwell’s equations in isotropic media,” IEEE Trans. Antennas Propag. 14(3), 302–307 (1966). [CrossRef]
38. A. B. Shvartsburg, “Optics of nonstationary media,” Phys.-Usp. 48(8), 797–823 (2005). [CrossRef]
39. W. Pauli, Theory of Relativity (Courier Corporation, 1981).
40. J. V. Bladel, Relativity and Engineering (Springer Science, 1984).
41. J.-P. Berenger, “A perfectly matched layer for the absorption of electromagnetic waves,” J. Comput. Phys. 114(2), 185–200 (1994). [CrossRef]
42. J. A. Kong, Electromagnetic Wave Theory (EMW Publication, 2008), 7th ed.
43. J. von Neumann and R. D. Richtmyer, “A method for the numerical calculation of hydrodynamic shocks,” J. Appl. Phys 21(3), 232–237 (1950). [CrossRef]
44. J. B. Schneider, Understanding the FDTD Method (www.eecs.wsu.edu/schneidj/ufdtd, 2010).
45. C. Caloz and Z.-L. Deck-Léger, “Spacetime metamaterials, Part II: Theory and applications,” IEEE Trans. Antennas Propag. 68(3), 1583–1598 (2020). [CrossRef]
46. K. S. Kunz, “Plane electromagnetic waves in moving media and reflections from moving interfaces,” J. Appl. Phys. 51(2), 873–884 (1980). [CrossRef]
47. L. A. Ostrovskiĭ, “Some “moving boundaries paradoxes” in electrodynamics,” Sov. Phys. Usp. 18(6), 452–458 (1975). [CrossRef]
48. A. Sommerfeld, Optics-Lectures on Theoretical Physics, Vol. 4 (Academic Press, New York, 1957). Sec.7.
49. M. Born and E. Wolf, Principles of Optics (Pergamon, 1980).
50. G. W. Milton and O. Mattei, “Field patterns: A new mathematical object,” Proc. R. Soc. A 473(2198), 20160819 (2017). [CrossRef]
51. T. Shiozawa and K. Hazama, “General solution to the problem of reflection and transmission by a moving dielectric medium,” Radio Sci. 3(6), 569–575 (1968). [CrossRef]
52. D. Censor, “Scattering of a plane wave at a plane interface separating two moving media,” Radio Sci. 4(11), 1079–1088 (1969). [CrossRef]
53. In the case of moving-perturbation, these structures are essentially traveling-wave modulation-type structures, while standing-wave modulation-type structures, such as those studied in [21, 50], which consist combinations of purely spatial and purely temporal interfaces, can be handled by the standard FDTD algorithm. Indeed, the FDTD algorithm automatically enforces continuity of the tangential E and H fields at stationary interfaces and continuity of D and B fields at temporal interfaces.
54. If we had $\mathbf{n}\smash{\|\llap{\big{/}}}\mathbf{v}$, the interface would be slanted with respect to the motion direction. Related problems were treated analytically in [51, 52].
55. In this case, for motion in the z direction ($\mathbf {v}=v\hat {\mathbf {z}}$), the tensors in Eq. (5) are [41] $\overline {\overline {\epsilon }}=\epsilon '\begin {bmatrix}\alpha & 0 & 0\\0 & \alpha & 0\\0 & 0 & 1\end {bmatrix}, \quad \overline {\overline {\mu }}=\mu '\begin {bmatrix}\alpha & 0 & 0\\0 & \alpha & 0\\0 & 0 & 1\end {bmatrix}, \quad \overline {\overline {\xi }}=\begin {bmatrix}0 & \chi /c & 0\\-\chi /c & 0 & 0\\0 & 0 & 0\end {bmatrix}, \quad \overline {\overline {\zeta }}=\overline {\overline {\xi }}^T,$ with $\alpha =\frac {1-\beta ^2}{1-\beta ^2{n'}^2}, \qquad \chi =\beta \frac {1-{n'}^2}{1-\beta ^2{n'}^2}, \qquad \beta =v/c\qquad \text {and} \qquad n'/c=\sqrt {\epsilon ' \mu '},$ where $\epsilon '$ and μ′ are constant scalars that represent the permittivity and permeability of the medium at rest (v = 0), which is assumed to be isotropic, nondispersive and linear.
56. To show this, we insert $\mathbf {B}=-\int \left (\nabla \times \mathbf {E}\right )\,\mathrm{d} t$, from (4), into (2), which yields, using the vectorial identity $\mathbf {A}\times (\nabla \times \mathbf {B})=\nabla (\mathbf {A}\cdot \mathbf {B})-(\mathbf {A}\cdot \nabla )\mathbf {B}-(\mathbf {B}\cdot \nabla )\mathbf {A}-\mathbf {B}\times (\nabla \times \mathbf {A})$, $\mathbf {E}^*=\mathbf {E}-\smallint \left (\mathbf {v}\times (\nabla \times \mathbf {E})\right )\,\mathrm{d} t=\mathbf {E}-\smallint \left (\nabla (\mathbf {v}\cdot \mathbf {E})-(\mathbf {v}\cdot \nabla )\mathbf {E}-(\mathbf {E}\cdot \nabla )\mathbf {v}-\mathbf {E}\times (\nabla \times \mathbf {v}) \right )\,\mathrm{d} t.$ The last two last terms in the last expression both vanish in the assumed regime of constant velocity. Given $\mathbf {E}=E_x\hat {\mathbf {x}}$ and $\mathbf {v}=v\hat {\mathbf {z}}$, the term with v · E also vanishes since v⊥E. Finally, the term $(\mathbf {v}\cdot \nabla )\mathbf {E}$ reduces to $\hat {\mathbf {x}}\partial E_x/\partial z$ and is hence parallel to E. Thus, $\mathbf {E}^*\|\mathbf {E}$. It may be similarly shown that $\mathbf {H}^*\parallel \mathbf {H}$.
57. To the best of our knowledge, there is no general method to synthesize a stencil with guaranteed stability in the discretization of a system of partial differential equations (PDEs). Such stencils are always found, as done here, from a systematic testing of different stencil structures (forward, backward or centered differences, and possible field averages) in the different PDEs forming the system [28, 29].
58. In Eq. (15b), one must be careful not to confuse kz, the z-component of the wavenumber, with k, the spatial index appearing in the FDTD update equations. No related ambiguity should exist in the paper insofar as the two quantities are always distinguishable by the presence or absence of the z subscript. Furthermore, there should be no possible confusion between the ξ and ζ functions in (15b) and the bianisotropic tensors $\overline {\overline {\xi }}$ and $\overline {\overline {\zeta }}$ in (7), since the latter always appear with double overbars and are not involved in the following stability analysis.
59. This is even true for a simple mismatched moving interface, where the incident and transmitted waves are co-moving and the reflected wave is contra-moving if the incident wave propagates in the direction of the modulation.
60. Note that, however, the uniform slab and frequency domain requirements of the transfer-matrix methods are very restrictive. The former requirement precludes the simulation of a space-time structure composed of slabs with finite temporal duration (and hence involving space-time corners), which are obviously a basic constraint in practical scenarios. The latter requirement essentially precludes the determination of the spatio-temporal features of the scattered waves because this would imply an extremely complex set of Fourier transformations. In contrast, the proposed generalized FDTD scheme can straightforwardly simulate truncated slabs and exact space-time waveforms.
