– Space-Time Engineered-Modulation (STEM) structures are currently experiencing massive interest in the metamaterials community, with the manipulation of the additional tem- poral dimension leading to novel effects and applications. However, there is presently no general numerical tool to simulate such structures. This paper closes this gap by providing a modified FDTD scheme solution, which involves hybrid – auxiliary (non-physical) and physical – numer- ical fields. It describes the corresponding modified Yee cell and derives the update equations. Moreover, it validates the proposed methods with the illustrative examples of a STEM slab and STEM crystal. I. INTRODUCTION We are interested here in Space-Time Engineered-Modulation (STEM) structures, namely space-time varying structures formed by modulating a host medium with a traveling-wave perturbation. Modulated structures of that kind have been abundantly studied in the past in the context of parametric amplification [1] and acousto-optics [2]. However, recent metamaterial developments have led to more sophisticated modulation-based structures – STEM metamaterials – whose modulation can be engineered for novel physical effects and novel applications (e.g., [3–7]). No general numerical technique is currently available for properly analyzing STEM metamaterials and struc- tures. The best strategy towards filling up this gap seems to be generalizing the Finite-Difference Time-Domain (FDTD) scheme, since FDTD, being in a time-domain approach, automatically predicts the space-time frequency transitions that characterize STEMs. However, a naive application of the FDTD scheme, which would consist in discretizing space-time interfaces in a staircase manner, turns out to be fundamentally flawed and yield wrong scattering amplitudes! An correct FDTD scheme is presented in [8] for the particular case of a moving conducting wall, but no solution is available even for such canonical problems as STEM penetrable interfaces, slabs, gradi- ent sections or basic metamaterials. Frame hopping has been attempted as an alternative technique to sidestep the issue [9], but this approach, being indirect, requires large computer memory and involves awkward numerical transformations. We propose here a general solution to this long-lasting lack of computational tool for STEMs. That solution consists in a simple modification of the conventional FDTD scheme: it replaces the electric and magnetic fields by auxiliary (non-physical) fields in the Yee grid while keeping the usual (physical) electric and magnetic flux density fields unchanged, derives the corresponding (modified) update equations involving hybrid (physical and auxiliary) fields, and finally transforms the auxiliary fields back into their physical counterparts. We derive here the scheme, without loss of generality, for the 1D case, and demonstrate its validity with two illustrative examples. II. T YPICAL STEM S TRUCTURES STEM structures, whether operating in the Bragg regime or in the metamaterial (homogenized) regime, may be modeled by dielectric layers delimited by moving interfaces [10], as illustrated in Fig. 1(a)1. Note that, although typically considered infinite in analytical studies, STEMs must be practically finite in size, and the related trunca- tion, unless comoving with the STEM structure – a practically rare situation! – already necessitates a numerical method per se. 1Although the figure shows sharp (step) interfaces, smooth (graded) interfaces can always be modeled by a sequence of thin layered media of vanishingly small width. v v (a) (b)x y zE∗ x E∗ x E∗ x E∗ x Ex H∗ y H∗ y H∗ y Hy ǫ1ǫ1 ǫ2ǫ2Fig. 1: STEM concepts used in the paper. (a) Example of a STEM structure. (b) Zoom on three Yee cells across an interface, showing the auxiliary (non-physical) fields in Eq. (1). III. M ODIFIED FDTD S CHEME In the conventional FDTD scheme and for a stationary structure, the electric and magnetic fields, which reduce toExandHyin the 1D problem considered here [see Fig. 1(b)], are sampled on the Yee cell [11]. As a result, the fields in the Yee grid are automatically conserved at dielectric and magnetic interfaces2. In contrast, in the case of a STEM structure, the quantities that are conserved at the interfaces are Ex+vByandHy+vDx[13]. This suggests to use these quantities, which we denote here E∗ x=Ex+vByand H∗ y=Hy+vDx, (1) instead of their stationary counterparts, in the Yee cell. This what is done in Fig. 1(b) across a moving interface of the STEM in Fig. 1(a), which automatically ensures the continuity of E∗ xat the interface, since this quantity is uniquely defined there. We next derive the update equations for the hybrid (auxiliary and physical) fields. We first insert the auxiliary fields (1) into the appropriate Maxwell-Faraday, Maxwell-Amp `ere and constitutive equations, which yields ∂By ∂t=−∂E∗ x ∂z−v∂By ∂z,∂Dx ∂t=−∂H∗ y ∂z−v∂Dx ∂z, (2) Dx=ϵ(E∗ x+vBy), B y=µ H∗ y+vDx , (3) with the two systems of equations involving the auxiliary fields E∗ xandH∗ yand the physical fields DxandBy. We can now discretize the equations (2) and (3). For the terms without the factor v, we apply the conventional FDTD central difference in space and in time and the leapfrog method (electric and magnetic fields evaluated at different times). For the terms involving the factor v, the discretization depends on the sign of the modulation velocity. For positive (negative) velocities, we choose a backward (forward) difference in space to ensure stability. The resulting scheme for positive velocities is By|n k+1 2=By|n−1 k+1 2−S E∗ x|n−1 2 k+1−E∗ x|n−1 2 k −Sv By|n−1 k+1 2−By|n−1 k−1 2 , (4a) H∗ y|n k+1 2=By|n k+1 2 µ|n k+1 2−v(Dx|n−1/2 k+1+Dx|n−1/2 k)/2, (4b) Dx|n+1 2 k=Dx|n−1 2 k−S H∗ y|n k+1 2−H∗ y|n k−1 2 −Sv Dx|n−1/2 k−Dx|n−1/2 k−1 , (4c) E∗ x|n+1 2 k=Dx|n+1 2 k ϵ|n+1 2 k−v(By|n k−1/2+By|n k−3/2)/2, (4d) where the fields at times t=n∆tand positions z=k∆zare denoted ψ|n k, and S= ∆t/∆z(Courant factor). 2When there is a discontinuity in both the permittivity and permeability, a simple average ensures the correct boundary condition [12]. IV . I LLUSTRATIONS Figure 2 provides two illustrations of the proposed method and compares the numerical results with exact, analytical results for the scattering coefficients versus frequency [14,15]3. These results are obtained by sending a very short temporal pulse into the structure and Fourier transforming the resulting transmitted and reflected fields. Figure 2(a) presents the case of a single moving slab while Fig. 2(b) presents the case of a crystal made of 5 unit cells. In both cases, the numerical results agree well with the analytical results. T ΓTΓ ωv 00 0 .5 1 1.50.20.40.60.81 (a) T ΓTΓ ωv 00 0 .5 1 1.50.20.40.60.81 (b) Fig. 2: Comparison of numerical (dotted) and analytical (continuous line) results for scattering coefficients of a STEM with ϵ1= 1,ϵ2= 2.5,L=λ0/2/n2andv= 0.3c. For the numerical results, S= 0.2and∆z=λ0/150. (a) Single slab. (b) Crystal with 5 bilayer unit cells.