Generalized Space-Time Engineered Modulation (GSTEM) Metamaterials

This article presents a global and generalized perspective of electrodynamic meta-materials formed by space and time engineered modulations, which we name Generalized Space-Time Engineered Modulation (GSTEM) Metamaterials, or GSTEMs. In this perspective, it describes metamaterials from a unified spacetime viewpoint and introduces accelerated metamaterials as an extra type of dynamic metamaterials. First, it positions GSTEMs in the even broader context of electrodynamic systems that include (non-modulated) moving sources in vacuum and moving bodies, explains the difference between the moving-matter nature of the latter and the moving-perturbation nature of GSTEMs, and enumerates the different types of GSTEMs considered, namely Space EMs (SEMs), Time EMs (TEMs), Uniform Space-Time EMs (USTEMs) and Accelerated Space-Time EMs (ASTEMs). Next, it establishes the physics of the related interfaces, which includes direct-spacetime scattering and inverse-spacetime transition transformations. Then, it exposes the physics of the GSTEM metamaterials formed by stacking these interfaces and homogenizing the resulting crystals; this includes an original explanation of light deflection by USTEMs as being a spacetime weighted averaging phenomenon and the demonstration of ASTEM light curving and black-hole light attraction. Finally, it discusses some future prospects. Useful complementary information and animations are provided in the Supplementary Material.


Introduction
Metamaterials are artificial structures consisting of supra-molecular but sub-wavelength particles that are engineered to provide medium properties beyond (µετά) those available in conventional materials [1,2,3]. Following rudimentary ancient nano-composites, medieval stained glasses and XX th century artificial dielectrics, they have experienced spectacular developments in the past two decades, where they have diversified and expanded to a point that they represent nowadays a powerful new paradigm in science and technology. This evolution has been largely facilitated by the advent of metasurfaces [4,5], which may be seen as the two-dimensional counterpart of voluminal metamaterials, with the benefits of easier fabrication, lower loss and greater flexibility, as well as drastic functional extensions of frequency-or polarization-selective surfaces, reflector transmit-arrays and spatial light modulators.
Most of the metamaterials and metasurfaces investigated until recently have been static, i.e., modulated only in space; we shall therefore refer to them as Space Engineered Modulation (SEM) metamaterials, or SEMs for short. A major advance in the field has been realized by making metamaterials dynamic, either by replacing the space modulation by a time modulation or by adding a time modulation to the space modulation. This introduction of the dimension of time as a new structural medium parameter has resulted in the metamaterial classes of Time Engineered Modulation (TEM) metamaterials, or TEMs [6,7,8,9,10,11], and Space-Time Engineered Modulation (STEM) metamaterials, or STEMs [12,13,14,15,16,17,18,19,20,21] 1 . Specifically, TEMs and STEMs are metamaterials that are formed by the variation (modulation) of a medium parameter in time and in both space and time, respectively, induced by an external drive. In the case of electromagnetic metamaterials, on which the paper focuses, the modulated parameter may be the refractive index, the permittivity, the permeability, or any of the bianisotropic and higher-order spatial-dispersion constitutive parameters and combination thereof, while the modulation drive may be acoustic (e.g., surface/bulk acoustic waves in a piezoelectric crystal), electronic (e.g., electric voltage variations in varactor chips), optical (e.g., laser pulses in semiconductor slabs), etc. [22,23]. TEMs and STEMs may thus be seen as medium -generally 3+1D, or 4D -extensions of electronic and optical active lumped element and circuit systems, such as parametric amplifiers [24,25] and acousto-electric/optic modulators [26,27].
This paper presents a global perspective of metamaterials and a generalization of dynamic metamaterials. The global perspective consists in describing all metamaterials, including SEMs and TEMs, in terms of space-time -or spacetime 2 -modulations, with various degrees of complexity, and in connection with the physics of moving bodies, while the generalization resides in the extension of STEMs with uniform (constant in both space and time) modulation velocity, i.e., Uniform STEMs (USTEM) metamaterials, or USTEMs, to STEMs with accelerated modulation, i.e., Accelerated STEM (ASTEM) metamaterials, or ASTEMs. We shall refer to these diverse possible types of metamaterials as Generalized STEMs (GSTEM) metamaterials, or GSTEMs, where it is noted that ASTEMs may feature different orders (derivatives) of acceleration and hence subdivide in further classes. Table 1 summarizes the terminology. 1 The terminology "time-modulated metamaterials" and "space-time modulated metamaterials" applies to metamaterials that have already a spatial modulation before being temporally modulated, but not to -equally relevant! -metamaterials whose dynamic structure is really formed (and engineered ) by a time or space-time modulation. Hence our introduction of the general terms TEMs and STEMs, and related terminology in Tab. 1. 2 The two spellings (with and without a hyphen) of this word are found in the literature on dynamic systems, whether for the noun or for the adjective. The one-word spelling is the universal standard when referring to the mathematical model that describes the merged nature of the space and time dimensions into a four-dimensional manifold in relativity physics (e.g., curved spacetime), while the spelling with a hyphen is preferable in reference to modulated structures, where the spatial and temporal features of the modulation are distinct and may exist independently of each other (e.g., space-time modulated metasurface). The present paper follows this convention.

Related Electrodynamic Systems
Electromagnetic GSTEMs are not the only electrodynamic systems. They only represent the category of moving-perturbation (or moving-modulation) electrodynamic systems. Two other fundamental types of electrodynamic systems should be considered here, vacuum moving-source systems and moving-matter (or moving-body) systems [28], because GSTEMs support physical effects that are inherited from them, although, as we shall see, in distinct embodiments. We shall next describe and compare the three categories, with the help of the illustrations provided in Fig. 1.

Vacuum
Matter Perturbation Vacuum moving-source systems, illustrated in Fig. 1(a), are systems involving objects (e.g., star or car) that emit or reflect light 3 while moving in vacuum relatively to the observer (e.g., Earth or road), with vacuum being defined as a portion of space that is essentially devoid of matter. The earliest reported related effect is Bradley aberration (top panel), whereby a terrestrial observer sees a star in a direction that is titled towards the direction of the motion of Earth in its orbit around the sun [29]. Another effect, which commonly manifests itself in daily life with sound sources, is the Doppler shift (bottom panel), whereby an observer of a moving source sees the frequency of the wave emitted or reflected by that source as depending on its velocity, with larger and lower frequency for approaching and receding motion, respectively [30]. Vacuum-moving-source systems are the simplest electrodynamic systems, since they are restricted to light propagation without light-matter interaction.
Moving-matter/body systems, illustrated in Fig. 1(b), are systems involving matter (e.g., water or dielectric) that moves relatively to the observer (e.g., laboratory experimenter) and that supports the propagation of light emitted from the reference frame of the observer, with matter motion defined as a collective translation or/and rotation of atoms and molecules over distances that are much larger than the molecular scale; these systems involve thus typically moving solids, fluids or gases. A related effect is the Fresnel-Fizeau drag (top panel), whereby the speed of light is reduced or increased for downstream or upstream propagation in a moving fluid [31,32,33]. Another effect that is of major importance in electrodynamics is the Röntgen magneto-electric coupling (bottom panel), whereby the motion (here, the rotation) of a solid submitted to an electric field induces a magnetic field in the frame of a rest observer, due to the creation of surface polarization currents [34,35]. Moving matter/body systems are more complex than vacuum moving-source systems because of the addition of their matter drag and magneto-electric coupling effects on top of the aberration and Doppler shift effects occurring in vacuum moving-source systems.
Finally, moving perturbation/modulation systems, illustrated in Fig. 1(c), are systems involving a perturbation (e.g., an acoustic wave in a piezoelectric crystals) that moves relatively to the observer (e.g., frame of an optical or microwave device) and that scatters light emitted from the reference frame of the observer, with perturbation motion defined as a traveling-wave (or standing-wave) modulation of some electromagnetic medium parameter, without any net transfer of matter, i.e., with motion restricted to oscillations of bound charges over submolecular distances (dielectric or magnetic polarization) 4 . A common example of such a system is the acousto-optic modulator (AOM) (top panel), whereby a periodic propagating perturbation ("spacetime modulation grating"), induced by variations of the molecular density of the medium from an electric signal (piezoelectricity), deflects the diffraction orders of the incident light in the direction of the perturbation via Bragg-Brillouin scattering [26,27]. GSTEMs (bottom panel), particularly USTEMs and ASTEMs, belong to this category of electrodynamic systems, where they generalize acousto-optic-modulator-type systems to multi-dimensional (2+1D = 3D and 3+1D = 4D), multi-velocity (uniform or nonuniform) [36], homogenized [21] and "new-physics" [19,20] electrodynamic systems. Figure 2 compares the electrodynamic structures of the moving-matter/body systems [ Fig. 1(b)] on one hand, and the moving perturbation/modulation systems in [Figs. 1(c)], which include GSTEMs, on the other hand. Figure 2(a) shows a moving-matter system (e.g., sliding curling stone), where the atoms and molecules (matter) move together with the body, along with the comoving frame, K , at a velocity v with respect to the (fixed) laboratory frame, K. Figure 2(b) shows a moving-perturbation system (e.g., dielectric slab excited by a laser-pump pulse or piezoelectric slab excited by a voltage source). In this case, the atoms and molecules oscillate about their bound position within the (solid) body, under the polarizing effect of the drive excitation, but do not experience any net motion in K. Only the related two perturbation interfaces move, inducing a STEM pulse modulation of the form n(z, t) = n 0 + ∆n · Π[(z − v m t)/D m ], where n 0 is the average refractive index, ∆n is the modulation depth, Π(·) is the pulse function, D m is the spatial extent of the pulse and v m = v is the velocity of the corresponding perturbation. Note that the atoms and molecules, while stationary in K, move with respect to K , in the opposite direction, with velocity v atom = −v. where v m = v mz = c/ sin θ, with even greater diversity attainable with a SEM or STEM mask at the output of the laser. Although the figures consider the specific example of laser-driven modulations, the principle holds for any type of modulation (optical, electronic, acoustic, mechanical, chemical, etc.) [22,23].

Moving Perturbation versus Moving Matter
Figure 2(c) shows a continuous (periodic) version of the pulse structure in Fig. 2(b) (e.g., using a periodically pulsed laser-pump or electro-acoustic drive), with the STEM medium function n(z, t) = n 0 + ∆n · sgn[(cos(k m z − ω m t) + 1]/2, where sgn(·) is the sign function, and k m and ω m , with v m = ω m /k m , are the spatial and temporal frequencies, respectively, of the modulation. Finally, Fig. 2(d) shows a variant of Fig. 2(c), where the drive excites the system from the top of the structure (e.g., oblique laser pump illumination or electro-acoustic source), under an angle θ with respect to the propagation axis, which provides superluminal modulation with velocity v m = v mz = c/ sin θ [19], reducing to instantaneous (pure time) modulation for θ = 0.
Moving perturbation/modulation systems [Figs. 2(b)-(d)], and particularly GSTEMs, are more promising than their moving matter/body counterparts [ Fig. 2(a)] towards reallife applications because i) they do not require cumbersome moving parts, ii) they easily attain relativistic velocities and accelerations, and iii) they posses richer functionality potential, resulting both from their dimensional extension of previous modulated systems and their capability to mimic and transcend cosmological systems (e.g., equivalent horizons and black holes; superluminality and negative mass equivalent) [37,38]. These are the reasons why GSTEMs are so attractive at this point of research in the field of metamaterials. We shall hereafter restrict our attention to GSTEMs, and refer to the moving matter/body dynamic systems only for the purpose of structural or property comparison.  [39] corresponding to the four main types of GSTEMs considered in the paper, with suggestive artistic illustrations (Supp. Mat. A). Such a global perspective offers multiple benefits, including i) an elegant classification, based on the natural concept of spacetime structuration, ii) a powerful unification, suggesting insightful comparisons and cross-fertilization concepts (e.g., time duals of space systems [40,10,41] or space-time extensions of pure space/time systems [42,18]), and iii) a generalization of the physics of special relativity [43,44,45] and general relativity [46,44,47] (Supp. Mat. B.1), which rather involves sources in vacuum [ Fig. 1(a)] or moving bodies [ Fig. 1

Perspective and Generalization
The concept of a continuous medium [48] is an idealization. In reality, all materials are formed by a more or less (crystal or amorphous) periodic collection of particlesatoms and molecules in the case of conventional materials and resonant scatterers in the case of metamaterials -which subtend the macroscopic response of the medium in terms of dipolar/multipolar responses at the microscopic scale. All materials can therefore be represented as periodically alternating regions of vacuum and particles, as suggested by the alternating gray-golden bands in the spacetime diagrams of Fig. 3. The case of GSTEMs seems a priori more complicated, because the modulation can take  . The subscript 'm' refers to 'modulation', n 1 and n 2 are the refractive indices of the constituent media, which are assumed to be isotropic and nondispersive, and z represents the spacetime hyperspace, which may include up to 3 spatial dimensions (x, y, z).
diverse and complex forms. Consider for instance, the sinusoidal modulation, n(z, t) = n 0 + ∆n cos(k m z − ω m t), which is common in practice. Even though this modulation continuously varies, the wave propagating in the medium is myopic to such fine detail in the metamaterial -sub-wavelength and sub-period -regime, corresponding to the twofold condition λ m λ and T m T (λ and T : smallest wavelength and period of the wave; λ m = 2π/k m and T m = 2π/ω m ); it only probes index extrema, with blurred transitions between them, and the modulation can therefore be safely approximated by the discrete binary function n(z, t) = n 0 + ∆n · sgn[(cos(k m z − ω m t) + 1]/2 [Sec. 3 and Fig. 2(c)] 5 . Thus, the spacetime diagrams of GSTEMs can generally be represented by the periodic bilayer spacetime patterns shown in Fig. 3.

Interface Physics
As illustrated by the spacetime diagrams in Fig. 3, GSTEMs can be modeled by alternating isotropic medium layers. The interfaces delimiting these layers are therefore the main discontinuities or non-uniformities of the structure and represent hence the entities that underpin the light-matter interaction of the metamaterial. For this reason, this section focuses on GSTEM interfaces, while Sec. 6 will reveal how the related principles extend to complete GSTEM media.
Let us start with the simplest cases of SEM and TEM interfaces. The electrodynamics of these interfaces is described in Fig. 4, with Fig. 4(a) and (b) representing the SEM and TEM cases, respectively [20]. When a wave hits a simple interface, or SEM interface, it splits into a reflected wave and a transmitted wave, which propagate in opposite directions over time, with well-known scattering (Fresnel) coefficients γ and τ , as shown in the left panel of Fig. 4(a). These coefficients are found by enforcing the conservation of the tangential E and H fields at the (spatial) interface discontinuity, which is required to avoid making the B and D fields singular (at the interface) through the spatial derivative (∇×) in Maxwell equations [20]. On the other hand, the transmitted wavelength is compressed (or the wavenumber increases) if the second medium is denser, as shown in the right panel of Fig. 4(a). Note that this transformation does not involve any change of temporal frequency (∆ω = 0, energy conservation) since the discontinuity is purely spatial.
The problem of a simple instantaneous interface, or TEM interface, is the perfect dual of that of the SEM interface. Now the incident wave splits into a later backward wave and a later forward wave, which also propagate in opposite directions, but in the same (later) medium and with different scattering coefficients [6], ζ and ξ, as shown in the left panel of Fig. 4(b). These coefficients are found by enforcing the conservation of B and D at the (temporal) interface discontinuity, which is required to avoid making E and H singular (at the interface) through the temporal derivative (∂/∂t) in Maxwell equations [20]. The two scattered waves are red-shifted (or their frequency decreases) if the second medium is denser, as shown in the right panel of Fig. 4(b). It is thus the temporal frequency that changes in this transformation, while the spatial frequency remains unchanged (∆k z = 0, momentum conservation), since the discontinuity is purely temporal. This space-to-time duality has several applications, including the inverse prism, a device that, instead of (a) (a) decomposing colors into angles as the Newton prism, maps angles into colors [40].
STEM interfaces are even more interesting than their SEM and TEM counterparts. The electrodynamics of USTEM interfaces is described in Fig [20]. Now the modulation occurs simultaneously in space and time, as in all the electrodynamics systems represented in Fig. 1. Scattering is space-like -with reflected and transmitted waves -for the subluminal case, and time-like -with later backward and later forward waves -in the superluminal case [20]. The corresponding scattering coefficients are found by enforcing the continuity of (E , H ) in the subluminal comoving frame (K , where ∆ω = 0) and of (D , B ) in the superluminal simultaneity frame (K , where ∆k z = 0) [42], which reveals, upon inverse-Lorentz transformation [53,45,54] [55,54,56,20]. On the other hand, the spectral transitions, whose reflective and transmissive parts are manifestations of the Doppler effect and of an index contrast effect, respectively, are oblique, since the discontinuity is both spatial and temporal, leading to simultaneous spatial and temporal frequency transformations [57]. In the case of oblique incidence, the dispersion curves, between which the oblique transitions occur, are altered by the change of the incident momentum on the interface and the scattered angles are deflected towards/against the direction of motion in the sub/super-luminal cases [42] (Supp. Mat. C).
Finally, an ASTEM interface may be seen as a generalization of a USTEM interface, where both the direct-spacetime interfaces and the normal-incidence dispersion lines change from straight to curved, as illustrated in Fig. 3(d) for the case of an ASTEM metamaterial with complex acceleration profile, including direction reversal, and hence jerk (∂a m /∂t = 0). While an accelerated system is locally uniform, so that special relativity and Lorentz transformation apply locally, it globally requires a much more complex treatment that belongs to the realm of general relativity and, hence, differential geometry [58,59,47]. In this case, the (linear) Lorentz transformations must be replaced by nonlinear transformations corresponding to the type of acceleration at hand [38,47], the simplest being the K -constant and rectilinear (or proper) acceleration, which is associated with Rindler transformations (Supp. Mat. B.3). Physically, the USTEM space-time transformation (diffraction-Doppler) effect is promoted to a space-time chirping effect [60,61], with still little explored physics and application potential. A first application is the recently reported ASTEM waveform generator, whose properly design acceleration trajectory allows virtually arbitrary pulse shaping [36].

Metamaterial Physics
Stacking the different GSTEM interfaces described in Sec. 5 leads to the formation of the corresponding GSTEM structures in Fig . Before homogenization, these structures are generally periodic spacetime structures (or, possibly, only locally quasiperiodic with a period gradient for the case of ASTEMs), or GSTEM crystals, and they involve the same direct (scattering) and inverse (transition) spacetime transformations as their interface counterparts. However, these crystals represent spacetime-extended structures, as opposed to spacetime-localized discontinuities, and support in addition multiple spacetime scattering and multiple spacetime transitions, which leads to specific spacetime crystal properties [18].
The GSTEM crystal can be spatially 1D, 2D or 3D, as shown in the bottom right inset of Fig. 3(a). If the wave of interest propagates in a space of dimension that is larger than the spatial dimension of the crystal [specifically, oblique incidence in a 1D crystal and off-plane incidence in a 2D crystal], the crystal is seen by the wave as anisotropic, with different tangential field components for different (e.g., p and s) wave polarizations. Since this represents the most general and most common configuration in early GSTEMs, let us consider henceforth this case of an anisotropic crystal, noting in passing that anisotropy, with an anisotropic medium being defined as a medium that exhibits different properties in different directions, is a purely spatial concept, since time is monodimentional.
In moving-matter dynamic systems [e.g., Fig. 2(a)], the problems are ideally solved in the comoving frame (K ), where everything is stationary, and hence the boundary conditions and the constitutive relations are both trivial. The situation is more compli-cated in moving-perturbation systems [e.g., Fig. 2(b)], because they involve motion in both frames [e.g., moving interfaces in K and (backward) moving atoms and molecules in K ]! Let us still choose the K frame, on the ground that moving matter with stationary boundaries is a known problem [62], whereas moving boundaries is a more difficult problem. To enter the metamaterial regime, we homogenize the GSTEM by prescribing subwavelength [max(λ m ) min(λ)] and subperiod [max(T m ) min(T )] operation. This operation leads to GSTEM-metamaterial anisotropic constitutive parameters with the same tensorial structure as that of their crystal parent [21,63]. However, another effect is present in K : the contradirectional motion of polarized atoms and molecules [e.g., focus on the moving K frame with K-stationary atoms and molecules in Fig. 2 The fundamental properties of a GSTEM metamaterial, in the frame of interest (K), may be inferred from these preliminary considerations with the help of the spectral graphs provided in Fig. 6, with Fig. 6(a) representing the problem from the K viewpoint and Fig. 6(b) from the K viewpoint. Since the Fresnel-Fizeau drag occurs in the −z (= −z) direction in K , as suggested in the panel at the extreme right bottom of the figure, the wave velocity, v gz = v pz = c/n (nondispersive medium assumption), is increased in the −k z (or −z ) direction and decreased in the +k z direction, i.e., 1/n − > 1/n + , so that the spectral cone is titled towards the positive k z direction [top panel of Fig. 6(a)]. As a result, the isofrequency curve at ω = ω 0 is a right-shifted ellipse [bottom panel of Fig. 6(a)], located between the n 1 (k ) and n 2 (k ) also right-shifted ellipses, which is a manifestation of the expected bianisotropy (off-centered ellipse → bianisotropy, centered ellipse → anisotropy, centered circle → isotropy) (Supp. Mat. E.2), and associated nonreciprocity [66]. The counterclockwise rotation, by the angle ϕ , of the group velocity, v g (k ) = ∇ k ω (k ), with respect to the phase velocity, v p = (ω /k )k , corresponds to the addition of a −z directed component to the velocity, and thus clearly shows the backward effect of the Fresnel-Fizeau drag.
We finally need to perform the required (Lorentz, Rindler, etc.) inverse transformation from K to K to complete the resolution of the electrodynamics problem at hand. Obviously, the elimination of the motion of the atoms and molecules in the transition from K to K eliminates motion-related bianisotropy, since v atom = 0. In contrast, the exact nature (bianisotropic, anisotropic, biisotropic, homoisotropic) of the K solution is difficult to predict and can be precisely found only by performing the exact K -to-K inverse transformation (inverse Lorentz transformation in the USTEM case of Fig. 6) (Supp. Mat. E.3). Figure 6(b) shows a typical result. Note that the n 1 (k ) and n 2 (k ) isofrequency ellipses have transformed to simple centered circles, corresponding to the expected curves for the assumed isotropic constituent media with (scalar) refractive indices n 1 and n 2 . Interestingly, the USTEM curve is still a right-shifted ellipse, with a deflection , compared with the cases of the bulk constituent media, with refractive indices n 1 (k ) and n 2 (k ) or n 1 and n 2 . (a) K viewpoint, with Fizeau-Fresnel drag, due to the motion of matter, in the −z direction. (b) K viewpoint, with "inverse drag", due to the motion of the interfaces.
of the group velocity by an angle of ϕ, still in the counterclockwise direction with respect to the phase velocity; this deflection is necessarily an effect of the moving-modulation interfaces, since matter motion does not exist anymore. We have here ϕ < ϕ , indicating that the deflection due to the modulation (K frame) is smaller than that due to the matter (K frame), but the effect is greater, and can reach ϕ > ϕ , if the constitutive media have less-than-one indices (n 1 , n 2 < 1) (plasma-type media). Most importantly, the observed light deflection is contradirectional to the modulation (+z), and hence opposite to the direction of the drag for a moving body, consistently with the finding in [67].
This GSTEM deflection effect is quite distinct from the Fresnel-Fizeau drag. Indeed, it does not involve any motion of matter that would "push" or "pull" -i.e., drag -light. It is rather an effect of spacetime weighted averaging, as first suggested in [18]. This effect is illustrated Fig. 7. In the SEM problem, represented in Fig. 7(a), light spends on average the same amount of time in medium 1 and in medium 2 in the forward and backward directions, so that the corresponding effective metamaterial indices are equal, i.e., n + = n − . In the K -USTEM problem, represented in Fig. 7(b), light is subjected to the conventional Fresnel-Fizeau drag, and propagates therefore faster in the backward (downstream) direction than in the forward (upstream) direction, so that n + > n − . In the K-USTEM problem, represented in Fig. 7(c), the spacetime slopes are the same as in the SEM problem, since no matter motion occurs, but the medium is spacetime-wise oblique, here tilting in the forward direction, and the following -rather subtle! -effect occurs due to this tilting. Consider, for simplicity, that n 2 n 1 , as in the figure [where the vertical axes have been denormalized to ensure a restricted graphical aspect ratio without introducing (unphysical) superluminal light curves]. Then, light spends much more time (see dashed lines) in medium 2, the slower medium, in the forward direction, where it propagates almost parallel to the medium trajectory, than in the backward direction, where the propagation is relatively more perpendicular to the medium trajectory, while the ratio of the forward to backward traveled distances increases in a much smaller ratio, as clearly apparent in the figure. As a result, n + > n − , and hence v + < v − , consistently with the observation in Fig. 6. This is really, as announced, a spacetime weighted averaging phenomenon; an exact mathematical formula for this phenomenon is provided in [18] (Supp. Mat. E.4).
As ASTEM interfaces are the curved-spacetime generalization of USTEM interfaces (Sec. 5), ASTEM metamaterials are the curved-spacetime generalization of USTEM metamaterials. Therefore, the principles exposed in conjunction with the USTEM graphs in Figs. 6 and 7 largely extend to ASTEMs metamaterials, although their rigorous treatment requires a quantum jump from the theory of special relativity [43,45], routinely applied to USTEMs, to the theory of general relativity [46,38]. Figure 8 presents two illustrative examples of ASTEM metamaterials [63]. Figure 8(a) shows a rectilinear ASTEM metamaterial, which exhibits the modulation-contradirectional group delay deflection [v g (r) S = E × H] and straight phase velocity [v p (r) k] propagation predicted in Fig. 6; the beam curvature can be here qualitatively inferred from piecewise (straight) USTEM deflection. Figure 8(b) shows a rectilinear ASTEM black hole, which attracts and absorbs light like a cosmic black hole, an effect that is unattainable in simple gradedindex lenses [68].

Future Prospects
Given their very fundamental nature and virtually unlimited diversity, GSTEMs have a formidable potential for scientific and technological innovation. The scientific prospects include 1) the study of the properties of higher dimensional (2+1D = 3D and 3+1D = 4D) GSTEM (unbounded) structures, 2) the analysis of the scattering and diffraction at the interfaces and wedges of spacetime-truncated [18] GSTEMs, 3) the exploration of new GSTEM geometries (e.g., Rindler, Schwartzshild, Kerr, jerk, snap, crackle, etc.) in both the homogeneous and Bragg regimes, 4) the extension of gravity analogs, currently restricted to interface horizons [72], such as white hole or Bose-Einstein condensate analogs, 5) the elaboration of new electrodynamic computational tools for spacetime nonuniformities, using for instance foliation decomposition [73], and for modulation-related multiphysics, and 6) the investigation of novel GSTEM physics (e.g., superluminal and interluminal scattering, spacetime reversal, generalized spacetime metrics, spacetime quantum and sub-cycle phenomena). The technological prospects include on the one hand the development of efficient modulation platforms and techniques (e.g., acoustic, electronic and optical) for the experimental implementation of the new GSTEM phenomena, and on the other hand the identification and demonstration of novel related applications. Many potential USTEM-related applications have been identified in [20]; we expect that these applications will generalize to ASTEM-type systems, with extra opportunities offered by various spacetime curvatures and generalized spacetime "chirping".

A GSTEM Spacetime Diagrams
The spacetime (or Minkowski) diagrams in Fig. 3, being (z, ct) patterns, include the complete information on the GSTEM structures. A useful complement to these representations are the time-evolution (t) of corresponding purely spatial -(z, x) -representations. Such representations are provided in the appended animation file GSTEM_Structures.mp4.

B Fundamentals of Relativity
A solid treatment of the electrodynamics of spacetime varying systems requires the utilization and adaptation of the principles and tools of the theory of relativity, or relativity. We shall therefore recall here the basic principles of this theory and establish the related essential tools.

B.1 Basic Principles
The fundamental principle of relativity is that physical phenomena measured in spacetime reference frames that are moving with respect to each other are different from each other, due to the different spatial and temporal perspectives of the frames. In this statement, spacetime is defined as a mathematical model that fuses the three dimensions of space and the one dimension of time into a single four-dimensional manifold, which is itself a geometrical topological space that locally resembles the Euclidean space [58]. Although not directly apparent in daily human experience, relativity effects play a major role in systems involving velocities approaching the speed of light, or relativistic velocities, and they underpin myriads of aspects of modern science and technology [38,54].
The theory of relativity, considered as a whole with Special Relativity (SR) and General Relativity (GR) combined, is based on the following four axioms: I First postulate of SR: The laws of physics are identical in all non-accelerated, or inertial, systems [43,44,45].
II Second postulate of SR: The speed of light in vacuum is the samec = 299, 792.458 m/s -for all inertial systems, regardless of the motion of the light source [43,44,45].
III Equivalence principle of GR: A uniform gravitational field is indistinguishable from a uniform acceleration field for a given reference frame [46,38,44,37].
IV Locality equivalence between SR and GR: The laws of physics in a gravitational or accelerated field (GR) are locally identical to their SR-counterparts at every point of spacetime [46,38,44,37].
Axioms I to III follow from experimental evidence, while Axiom IV is a consequence of the local flatness or differentiability of the spacetime manifold [59]. We shall deal here with the electrodynamics of both inertial -and hence force-lessand accelerated spacetime systems, which we shall hereafter refer to as uniform-velocity (UV) (as USTEM) and nonuniform-velocity (NUV) (as ASTEM) systems, respectively, in reference to an assumed inertial observation frame. UV and NUV systems respectively pertain to SR [43,44,45] and GR [46,38]. They are hence governed by Axioms I and II for the former case and by Axioms III and IV for the latter case. The first task in the study of the electrodynamics of such systems is to derive transformation relations between different reference frames for both the (direct and spectral) spacetime variables and the electromagnetic fields; this is accomplished in Sec. B.2 for the UV case and in Sec. B.3 for the NUV case.
The derivation of the variable and field transformations will be done with the help of the reference frame setup shown in Fig. 9. In this system, K (here, Earth) is the frame where the experiment is conducted, and is hence called the laboratory frame, while K (here, rocket) is a frame that moves at a velocity v = vẑ relatively to K, and is hence called the moving frame, or the rest frame because comoving entities are at rest in that frame, despite the motion relatively to K . The figure also shows a third element (here, spatial probe) in motion with velocities u and u relatively to K and K , respectively, and propagating light, which is measured as having the same speed c in K and K if these frames are inertial, according to Axiom II. Figure 9: Relative reference frame setup used to describe spacetime systems.

B.2 Uniform-Velocity (UV) Transformations
Historically, the UV [v = v(r, t) or v-constant] transformation relations between the physical quantities in K and K were first derived by Lorentz as the mathematical conditions for the form invariance of Maxwell's equations between inertial frames [53], and then shown by Einstein to follow directly from Axioms I and II [43]. The literature on the derivation of the these relations is rather obscure. We shall present here a clear and rigorous derivation, for both the spacetime variables and the electromagnetic fields, along with the velocity addition formula. A priori, the relations between the spacetime variables (z, t) of K and (z , t ) of K (Fig. 9) could be anything. However, the homogeneity and isotropy, or the 4D symmetry, of spacetime that prevails in inertial systems demands that the relation between K and K be the same at every point of spacetime. Therefore, where A, B, C and D are constants. This indicates that the sought after relations must be linear.
The constants A, B, C and D in Eqs. (1b) may be found by successively enforcing the following physical conditions: 1) motion of the origin of K (z = 0) at velocity v relatively to K, 2) motion of the origin of K (z = 0) at velocity −v relatively to K , 3) uniformity of the speed of light, c (Axiom II), 4) time reversal (t, t → −t, −t ) symmetry [28] upon successive direct (velocity v) and reverse (velocity −v) transformations, i.e.,

Solving this system of equations yields
and whose substitution into the last relations of Eqs. (1b) leads to the Lorentz transformation relations where we have added trivial relations in the motion-less directions x and y for generalization to the observation direction r = xx + yŷ + zẑ. Note that the direct and inverse relations [left and right in Eqs. (2b), resp.] are obtained from each other by exchanging the primed and unprimed quantities and changing the sign of v. Also note that −(cdt) 2 + dr 2 = −(cdt ) 2 + dr 2 , where r ( )2 = x ( )2 + y ( )2 + z ( )2 , which reveals that 6 E.g., ∂ ∂z ∂t ∂z = 0 because the relations between t and z cannot depend on z since space is assumed to everywhere be the same; similarly, ∂ ∂t ∂t ∂t = 0 because the relation between t and t cannot depend on t since time is assumed to be always the same. the quantity ds 2 = −d(ct) 2 + r 2 , which is called the Minkowski metric 7 , is an invariant of spacetime. Figure 10 plots the Lorentz reference frames corresponding to Eqs. (2b), from the K viewpoint in Fig. 10(a), where the axes (z , ct ) are obtained by setting (ct , z ) = (0, 0) in the direct relations, and from the K viewpoint in Fig. 10(b), where the axes (z, ct) are obtained by setting (ct, z) = (0, 0) in the inverse relations. Projecting the intervals between successive crests or troughs of the monochromatic plane wave drawn in the background onto the axes (z, z ; ct, ct ) immediately provides the corresponding wavelengths (or spatial periods) and (temporal) periods (λ, λ ; cT, cT ). The differences in the two frames illustrate the Doppler effect, whereby the wave, which propagates in the +z-direction in the first quadrant, is seen with larger frequency (ω = 2π/T , blue shift) from the −z-direction, relatively counter-moving frame K, and with smaller frequency (red shift) from the +z-direction, relatively co-moving frame K . It is useful to also establish the spectral (Fourier-transform) counterparts of the spatial Lorentz relations (2). These relations may be easily found from the latter by noting that K and K observers agree on the phase of a monochromatic plane wave at any point P , with spacetime coordinates (z, ct) ≡ (z , ct ) (Fig. 10), as corresponding to a crest, a trough, or any other level of the wave amplitude. Thus, we have φ (z , ct ) = φ(z, ct), where φ = k z z − ωt and φ = k z z − ω t are the phase of the wave at P in K and K , respectively, with k ( ) z = 2π/λ ( ) z and ω ( ) = 2π/T ( ) being the corresponding wavenumbers (or spatial frequencies) and (temporal) frequencies, respectively. Substituting the inverse relations in Eq. (2b) in the resulting relation, k z z − ωt = k z z − ω t , factoring out z and t , and noting that the last relation must hold for any z and t , yields then where the transformation between the direct and inverse relations is obtained, again, by exchanging the primed and unprimed quantities and by changing the sign of v. The last relations in Eq. (3) provide the exact formulas for the Doppler frequency shift; specifically, the shifted frequencies are ω = γω (1 ± |β| cos θ ) or ω = γω(1 ∓ |β| cos θ), where θ ( ) is the angle between the observation direction (r ( ) ) and the motion direction = ∓|β| cos θ. Spectral spacetime diagrams, described in [20], provide a powerful graphical technique to determine both the temporal and the spatial frequency transitions induced by moving systems.
The field transformation relations may now be simply found by applying Axiom I to electrodynamics, viz., by enforcing the form invariance of Mawell's equations in K and K . This is accomplished by expressing the spacetime derivatives versus K in terms of the spacetime derivatives versus K , obtained from Eq. (2b) as inserting these relations into the K-frame Mawell's equations, and identifying with the K -frame Mawell's equations. This gives the Lorentz field transformations with the usual connection between the direct and the inverse relations. Let us finally derive the velocity addition formula. For this purpose, consider an object moving at the speed u = dz /dt in K (probe in Fig. 9). Taking the ratio of the differentials of the inverse Lorentz transformation formulas in Eqs. (2b), namely dz = γ(dz + vdt ) and dt = γ(dt + β c dz ), gives the velocity of the object in K [43], which differs from the simple (Galilean) sum (u +v) by the factor βu /c in the denominator. In the case of a relatively small K -frame velocity (β 1) and moderate u velocity, the denominator of the last expression in (7) can be approximated by its first-order Taylor series, which, assuming an isotropic propagation medium with refractive index n so that u = c/n, simplifies the result to the original (approximate) formula of Fizeau [32].
Generalizations of the formulas derived in this section to arbitrary vectorial velocities, v, are available in many textbooks (E.g., [54,56]).

B.3 Rindler Transformations
There is a virtually unlimited diversity of possible NUV [v = v(r, t)] systems, as may be realized by considering their numerous gravitational counterparts (Axiom III) in cosmology, without even mentioning modulation-specific (e.g., equivalent negative mass or non-progressive motion) extensions. Examples include spacetime systems with spherical (Schwarzschild) symmetry, charged spherical (Reissner-Nordström) structure, or charged rotationally-moving spherical (Kerr-Newman) gravitation or acceleration features, which are commonly studied in GR [38,37,47]. Here, we shall restrict our attention to NUV systems, K , that involve rectilinear (locally identical to spherical) motion relatively to an inertial observer, K, i.e., v = v(r, t)ẑ, corresponding to the acceleration a = (dv/dt)ẑ (Fig. 9). This restriction, although considerable, will fortunately not prevent us from capturing the essence of the phenomena occurring in NUV systems.
In addition to the restriction of rectilinear motion, we shall assume that the acceleration of the considered NUV entities is independent of time in their reference frame 8 , i.e., a = a (t ), while allowing it to vary in space, i.e., a = a (z ). This is equivalent to assuming that the accelerated entities experience a force that may vary in space but that is time-wise constant any point of space, an assumption that is almost universally done in GR, both because no simple transformation relations could be found otherwise and because it corresponds to the physical reality of typical gravitational systems 9 . This assumption is also reasonable for an STM system, whose acceleration may depend on the position, e.g., via a nonuniform background medium, but will typically not depend on time at a given position.
Most textbooks on GR (e.g., [37,38,47]) derive the NUV transformation relations from metric tensors associated with differential geometry. This approach, while being 8 Time-wise constant acceleration implies ultimately exceeding the speed c. This is not at odds with Axiom II since this axiom applies only to inertial frames whereas K is noninertial. The speed of light (and objects) in noninertial frames is not limited to c. For instance, in the counter-rotating direction of a vacuum-based Sagnac system [33], it exceeds c in the rotating frame. Here, u is limited to c because K is inertial, but u is not since K is noninertial. 9 The gravitational force of a celestial body might depend on the position of the observer around it, if its mass is nonuniformly distributed (e.g., moon), but it does not depend on the time of observation at any given position.
the most general and the most powerful existing one, offers unfortunately little insight into the physics of NUV systems. To avoid this inconvenience, and to maintain the self-consistency spirit of this appendix, we shall present here an intuitive, yet rigorous, derivation of the NUV transformation relations, as in the case of uniform velocity.
The Lorentz transformations [Eqs. (2)] are linear because of the symmetry of spacetime that prevails in inertial systems. Such symmetry is broken in noninertial systems, as may be realized by thinking for instance of the gravitational potential of a massive object, which decays as 1/r from its center, from a high level at its surface to virtually zero farther in deep space. Therefore, the NUV relations between K and K must necessarily be nonlinear. Consequently, we must adopt a more elaborate derivation strategy than in the UV case.
We shall determine the NUV transformation relations by invoking Axiom IV, viz., by applying the UV transformation relations [Eqs. (2b)] locally in the NUV system of interest. This will be done with the help of Fig. 11. As we have seen in Fig. 10(a), the UV coordinate axis ct corresponds to the trajectory of a UV object that is co-moving with K in K. Similarly, Fig. 11(a) plots the trajectory of an accelerated object (e.g., spatial probe in Fig. 9) in K, cτ , which will ultimately become the NUV ct axis upon application of the differential forms of the UV relations for locality. This trajectory has a curvature, which corresponds to the assumed constant acceleration (and hence constant force) in K , a 0 . Consequently, the object will accelerate along its trajectory with acceleration a(t) until it reaches the speed of light, u(t) → c, in K, as indicated in the figure. According to Axiom IV, the laws of SR (Sec. B.2) may be applied locally to the GR problem at hand. Specifically, they may be applied at the origin of K in Fig. 11(a), where the accelerated (NUV) frame, (z , ct ), locally reduces to the Lorentz (UV) frame K (O ; z , ct ) shown in the figure. At that point, we have u = v (same velocity for the probe as for the rocket relatively to K), so that u = 0 (Fig. 9). Taking the derivative of the SR velocity addition formula (7) with respect to t leads then to the local relation (Sec. F.1) between the accelerations in K (a) and in K (a 0 ) at the origin.
Integrating the relation (8) from 0 at an arbitrary point P on the trajectory cτ , specifically integrating over the velocity from v to u, over time from 0 to t and over space from 0 to z, delocalizes the problem and leads to the global velocity and position functions along the entire trajectory cτ (z, t) in K (Sec. F.2), and The shape of the trajectory corresponding to Eq. (9b) may be better visualized upon moving the t-dependent right-hand side square root to the left-hand side and squaring the resulting equation. This yields which is the equation of a hyperbola centered at the point X ≡ (−c 2 γ/a 0 , −c 2 γβ/a 0 ) and with asymptotes (ct) = ±(z +c 2 γ/a 0 )−c 2 γβ/a 0 , in K, as shown in Fig. 11(a), where only the right branch of the hyperbola is considered given the assumption of motion in the +z-direction.
We can now globalize the problem from the hyperbola cτ to the entire spacetime plane. For this purpose, we consider a UV local frame (O ; z , ct ) that moves along the trajectory cτ , as indicated in Fig. 11(b), and describe the points Q(cτ ), noted Q n in the figure 10 , in the spatial neighborhood of O by the coordinates (z , ct = 0) in that local frame. We use then a cτ -parametrized version of the local Lorentz relations [Eqs. (2b)], z(z , ct ) =γ(τ )(z + vt ) and ct(z , ct ) =γ(τ )(ct +β(τ )z ), whereγ(τ ) andβ(τ ) are the τ -parametrized versions of γ and β in Eqs. (2a) (τ = 0) , with t = 0, and find the following form for the K-frame coordinates of Q: where z O (cτ ) and ct O (cτ ) are the coordinates of O and z is the distance of Q to O , as apparent in the figure. The cτ parametrization is performed, invoking again Axiom IV, by first integrating the local relation dt =γ(t)dt with v = u(t) over time from 0 to t for t and from 0 to τ for t in the left-and right-hand sides, respectively, which yields the t = τ mapping relation (Supp. Mat. F.3) and inserting then this relation into the cτ -trajectory velocity formula (9a) to obtain the functionsγ(cτ ) andβ(cτ where z 0 = c 2 /a 0 , z a = γz 0 , z b = γβz 0 and ξ = sinh −1 (γβ), and where γ and β are still given by Eq. (2a). Equations (13) are the Rindler-Kottler-Møller transformation relations [69] 11 , which are the sought after NUV counterparts of the Lorentz relations (2).
The following comments are in order about Fig. 11(b) and Eqs. (13). First, each hyperbolic coordinate curve corresponds to a time evolution (or aging) that is specifically experienced by a hypothetical object following it; therefore, the different times experienced along the different hyperbolas are called the proper times, with clocks ticking slower and faster on the left and right of the ct axis, respectively. Second, each oblique straight coordinate line (z ) corresponds to locations where all the points of K are synchronized; these lines are therefore called simultaneity planes. Third, the z = constant coordinate lines are the right branches of the hyperbolas described by Eq. (10) with a 0 replaced by the parameter a in the right-hand side of the relation and having the vertices (−c 2 γ/a 0 + c 2 /a 2 , −c 2 γβ/a 0 ) (same center and asymptotes but different vertices); inserting the inverse Lorentz relations into the new equation, setting t = 0 (time-invariant acceleration) gives the z -dependent acceleration where a (0) = a 0 , as expected. Finally, we have −(cdt) 2 +dr 2 = −[(z +z 0 )ζ] 2 (cdt ) 2 +dr 2 with ζ = z/z 0 (Sec. F.5), where ds 2 = −[(z + z 0 )ζ] 2 (cdt ) 2 + dr 2 is the Kottler-Møller metric, corresponding to the Minkowski metric ds 2 = −(cdt) 2 + dr 2 already encountered in Sec. B.2). Figure 12 shows the corresponding frame-wave representation, as the NUV counterpart of the UV representation in Fig. 10, from the K viewpoint in Fig. 12(a) and from the K viewpoint in Fig. 12(b). The curvature of spacetime, specific to GR [46,38,44,37], is immediately apparent in these graphs. The figure also shows some wavelengths (or spatial periods) and (temporal) periods (λ, λ ; cT, cT ). In contrast to those of the UV case ( Fig. 10), these quantities vary in spacetime, corresponding to the transformation of a monochromatic and monodirectional wave into a spacetime frequency rainbow 12 . Particularly interesting is the apparition of an event horizon, called the Rindler horizon for the present case of rectilinear, constant proper acceleration, which is a boundary of spacetime that can be crossed only one way, as the boundary of a black hole [38,37]. Given the aforementioned rainbow transformations, simple spectral relations such as the UV ones in Eqs. (3) do not exist. In contrast, closed-form expressions for the field transformation counterparts of Eqs. (16b) exist, and they may be easily derived, invoking again Axiom IV by locally applying the Lorentz field relations (16b). Consider for instance the first direct Lorentz relation in Eqs. (16b), which becomes, upon parametrizing γ as γ(t) = ∂t/∂t (boost) and v as v(t) = ∂z/∂t, E x = (∂t/∂t )E x − (∂t/∂t )(∂z/∂t)B y = (∂t/∂t )E x − (∂z/∂t )B y . In these relations, the derivatives are found from the Rindler-Kottler-Møller relations (13) as ∂t/∂t = (z + z 0 )ζ cosh(ζt + ξ)/c and ∂z/∂t = (z + z 0 )ζ sinh(ζt + ξ). The corresponding field transformation relation is then which, upon noting that (z + z 0 )ζ/c = (z + z 0 )(c/z 0 )/c = (z /z 0 + 1) = (1 + a 0 z /c 2 ), factoring out cosh(ζt + ξ), and further noting that β = v/c = (dz/dt)/c = tanh(ζt + ξ) [using the differentials the dz and cdt following from (13a)] and thus γ = 1/ 1 − β 2 = cosh(ζt + ξ), may also be written as The other relations may be found in a similar fashion. The results are, using the GRmetric notation In the case of an x-polarized transverse electromagnetic (TEM) wave propagating in free space, where B y = E x /c, and in the limit t → ∞, the relation (15b) reduces to This relation reveals that, for a given E x measured in K, E x increases (linearly) with z and decreases (exponentially) with t . These double effect may be explained from the Doppler effect, whereby the electromagnetic field compresses in spacetime and hence increases in magnitude with decreasing co-moving relative velocity between the source and the observer, which occurs for increasing z and decreasing t , as clearly visible in Fig. 12(b).
C Scattering Angles at a USTEM Interface under Oblique Incidence Figure 13 provides a graphical method to determine the scattering angles at a USTEM interface for oblique incidence in the subluminal regime. Figure 13(a) shows two cuts of the dispersion diagram: the ω/c − k z plane (top), corresponding to the dispersion diagram, and the k x − k z plane (bottom), corresponding to the wavevector diagram. The dispersion diagram is obtained by the following steps: 1. plot the hyperbolas corresponding to ω 2 n 2 1,2 /c 2 − k 2 z = k 2 xi , where k xi is the transverse (k x ) component of the incident wave wavevector; 2. plot the axes of the laboratory K frame and the primed K frame; 3. locate the incident wave spectrum coordinates ω i , k zi ; 4. apply the conservation of frequencies in the primed frame, ∆ω = 0, to deduce the ω, k z of the reflected and transmitted wave. This is a generalization of the frequency transition plot of Fig. 5(a) (right), with the only difference being that the dispersion curves are now hyperbolas instead of straight lines.
The wavevector diagram is obtained by the following steps: 1. draw the wavevector of the incident wave, which has a known angle; 2. project the results of the dispersion diagram onto the wavevector diagram; 3. apply the continuity of k x to deduce the scattering angles of the reflected and transmitted waves; 4. add the isofrequency circles corresponding to k 2 x + k 2 z = ω 2 n 2 /c 2 (optional, added for completeness).
Following this approach for a stationary interface and a modulated interface moving to the right (+z), as in Fig. 13(a), we see that modulation to the right leads to a clockwise rotation of the reflected and the transmitted waves. Figure 13(b) presents the results in the direct space-time, by representing two cuts of the space-time trajectory of the waves: the z−ct plane (top) and the x−z plane (bottom). The angles in the x − z plane are those from the bottom of Fig. 13(a), and the slopes at the top are calculated as n z = c/v z = c/(v cos θ) = n/ cos θ. This is a generalization of the space-time scattering of normal incidence, Fig. 5(a) (left), to oblique incidence.
Similar graphs can be found in [42] where only the reflected wave is derived. The deflection angles can be derived by inserting the Lorentz transformations of the wavenumber and the frequency (3) into the phase-matching conditions, k xi = k xr = k xt and ω i = ω r = ω t . After some manipulations, we find then cos θ r = (β 2 n 2 1 + 1) cos θ i − 2βn 1 β 2 n 2 1 − 2βn 1 cos θ i + 1 and Figure 14 shows the counterpart of Fig. 13 for the case of a superluminal modulation. Similar steps as those leading to Fig. 14(a) are applied, with some variations: 1. plot the hyperbolas of the two media for a given k xi ; 4. project onto the wavevector diagram, setting ∆k x = 0 as previously.
(19b) Full-wave (FDTD) animations related to this section are provided in the appended file GSTEM_Interface_scattering.mp4.

D.1 Constitutive Relations
We assume an isotropic and nondispersive medium in uniform motion with the constant rectilinear velocity v = vẑ. In the rest frame of the medium (K ), the constitutive relations are The corresponding relations in the laboratory frame (K) are found by applying the Lorentz transformations (16b) to the fields in (20) and rearranging the result so as to express D and B in terms of E and H. This results into the bianisotropic relations [56] Note that, according to (21), the response fields (B and D) are not parallel to the excitation fields (E and H) in a bianisotropic (and even homoanisotropic, i.e., without magnetoelectric coupling) medium. As a result, the wave propagating in such media are inhomogenous, with the phase velocity ( D × B) and group velocity ( S = E × H) being nonparallel.
Also note that if the medium at rest is a natural medium (as opposed to a metamaterial) and is not subjected to any external force (e.g., external magnetizing field) or operating in the nonlinear regime [66], then it is necessarily nonmagnetic, so that µ ≈ µ 0 .

D.2 Inference of First-Order Bianisotropy from the Lorentz Force
The (exact) bianisotropic relations (21) were derived mathematically, using Lorentz transformations, with little insight into the related physics. We provide here an original, first-order derivation of these relations, which shows bianisotropy (or Röntgen magnetoelectric coupling) is an effect that is dominated by the action of the Lorentz force.
The Lorentz force reads [28] where and induces the forces illustrated in Fig. 15 on the atoms and molecules of matter, with the electric part, F E , elongating the particle to produce the electric dipole moment response p, and the magnetic part, F H , forming a current loop to produce the manetic dipole moment response m, which are both represented in the inset of the figure. The different components of the so-induced electric and magnetic dipole moments, p Sv u and m Sv u (u, v = x, y, z), respectively corresponding to the response R = D, B to the excitation S = E, H, result at the macroscopic level into susceptibility components χ RS uv , which can be concisely deduced as follows: where the field pair (E, H) was rotated by π/2 about the z direction in the bottom relations. Ignoring in (24) the terms in O(β 2 ), whose effect is of second order and whose exact form is also too complex to establish without an exact relativistic treatment, we still, remarkably, retrieve the first-order expression of (21), corresponding to the approximation [α, χ] = [1, β]. This reveals that bianisotropy is, to the first-order, a direct result of the Lorentz force.

E.2 Dispersion Relation for a Crystal in the Moving Frame (K')
We shall homogenize a space-time modulated 1+1D crystal with modulation in the +z direction. In the frame that is co-moving with the modulation, K , the interfaces of the crystal are stationary, so that they will be subject to the standard boundary conditions. On the other hand, the media between the interfaces are moving, in the −z direction, and will therefore have the bianisotropic properties given in Supp. Mat D.1. In K , the (standard) continuity conditions for an s-polarized wave are and those for a p-polarized are We can find the constitutive relations for the homogenized structure by i) assuming that the quantities conserved at the interfaces [Eqs. (26a) and (26b)] do not depend on the position (z), but are rather constant, throughout the structure, as must be the case in the homogeneous regime, ii) forming weighted averages, in terms of spacetime travel lengths, m (n = 1, 2) and t = 1 + 2 (Fig. 7) [18], of the other quantities , and iii) using the bianisotropic relations (21) for the fields D and B.
For s-polarized waves, we have with n zav = av µ av (33a) and n xav = av µ zav or zav µ av , for s-and p-polarization, respectively.

E.4 Spacetime Weighted Averaging
To gain insight into the dependence of the effective refractive index on the angle, we find the effective refractive index for the pure-downstream and pure-upstream cases k = k zẑ and k = −k zẑ by a simple geometric argument, using Fig. 7. After some manipulations, given in [18], we find the expression n ± av = n 1 1 + n 2 2 ∓ vn 1 n 2 t /c 1 + 2 ∓ v(n 2 1 + n 1 2 )/c (34) which was numerically verified to correspond exactly to (32a), for the case k x = 0, when F Derivations for Sec. B.3