Modeling four-dimensional metamaterials: A T-matrix approach to describe time-varying metasurfaces

: Exploring the interaction of light with materials periodically structured in space and time is intellectually rewarding and, simultaneously, a computational challenge. Appropriate computational tools are urgently needed to explore how such upcoming photonic materials can control light on demand. Here, we introduce a semi-analytical approach based on the transition matrix (also known as T-matrix) to analyze the optical response of a spatiotemporal metasurface. The metasurface consists of a periodic arrangement of time-varying scattering particles. In our approach, we depart from an individual scatterer’s T-matrix to construct the effective T-matrix of the metasurface. From that effective T-matrix, all observable properties can reliably be predicted. We verify our semi-analytical approach with full-wave numerical simulations. We demonstrate a speed-up with our approach by a factor of more than 500 compared to a finite-element simulation. Finally, we exemplify our approach by studying the effect of time modulation on a Huygens’ metasurface and discuss some emerging observable features.


Introduction
With the advent of nanofabrication technologies [1][2][3], metasurfaces have gathered much attention within the scientific community.They are essentially arrays of scatterers arranged in two dimensions (2D), which are highly tunable in terms of their optical response.While originally metasurfaces made from metallic constituents were at the focus of interest, it has been recognized that scatterers made from dielectrics or semi-conductors can also affect light propagation in a deterministic manner.Amplitude, phase, polarization, and spectral contents in reflection and transmission can be controlled.So far, a lot of fruitful research has been performed to tailor such metasurfaces to achieve desired functionalities for specific needs [4][5][6], with tunable metasurfaces playing a prominent role in that regard [7][8][9][10][11][12][13][14].Also, stacking metasurfaces to achieve metamaterials or, more general, photonic bulk materials that control light propagation is quite a mature research domain by now.However, all these studies were mostly restricted to time-invariant materials.
Whereas tunable metamaterials usually suffer from a low tuning speed in comparison to the operating frequency of light, there have been recent experiments that demonstrated a rather fast and considerably strong temporal modulation of the electric material properties of transparent conducting oxides (TCOs) [15][16][17].Consequently, time-varying materials and structures have recently emerged as a promising field of research [18][19][20][21][22][23].Our article is dedicated to this burgeoning topic.The time-variation of the scatterer's properties provides an additional degree of freedom to tune the optical response of the metasurfaces and metamaterials and opens the opportunity to reach four-dimensional metamaterials [24,25].Many interesting findings have been reported regarding the optical response of time-varying materials.These include anisotropic antireflection temporal coatings [26,27], photonic time crystals (PTCs) [28,29], nonresonant tunable PTC lasers [30], non-reciprocity [31,32], multifrequency perfect absorbers [33], power combining of waves [34], dynamic polarization converters [35], time-varying optical vortices [36], time-varying epsilon-near-zero materials [37,38], and temporal aiming [39].
In our contribution, we are specifically interested in developing a comprehensive theory to describe the optical response from a metasurface made from a periodic arrangement of time-varying scatterers.Two developments are primers for this contribution.On the one hand, the interaction of light with an individual, localized scatterer made from a time-varying material has already been solved semi-analytically [40][41][42].A prototypical example of such an object would be a sphere.Furthermore, efficient computation of the optical response of an array of time-invariant scatterers has also been demonstrated [43,44].However, the extension of the multiple-scattering formalism to the case of arrays of time-varying scatterers is still missing.Several full-wave numerical solvers exist [45] that can provide the scattering output for such clusters, but they require large computational resources and are considerably slow.These requirements limit their usability since most metasurface applications rely on repeated calculations to optimize for desired tailored functionality.Therefore, a semi-analytical tool is urgently needed for much more efficient computations, and it is developed and presented herein.
In this contribution, we semi-analytically solve the multiple-scattering problem of an infinite periodic array of time-varying scatterers arranged in 2D.We emphasize upfront that dispersive materials are considered.That is necessary because a material whose properties can tremendously be modulated in time can likely only be achieved by dispersive materials.Our semi-analytical approach is generally based on the T-matrix methodology.First, we introduce time modulation in our system by employing the Lorentz oscillator model with a time-varying electron density.Then, we use the existing knowledge of the T-matrix of a single time-varying scatterer [40] to construct an effective T-matrix of a spatiotemporal metasurface.We use the Ewald summation technique to compute the effective T-matrix of such a metasurface efficiently [43,46,47].It is an effective T-matrix because, in essence, the properties of the individual scatterer are renormalized by the interaction with all the other scatterers forming the infinite 2D arrays of particles.This technique has already been successfully applied for time-invariant metasurfaces [43].While our approach is general in its formulation, we assume spherical scatterers for which we know their T-matrices analytically.Once we know the effective T-matrix, all observable properties, such as reflection and transmission in each spatial diffraction order and at each frequency component generated by the time-varying metasurface, can be predicted reliably.
Next, we verify our theory against numerical simulations with a full-wave solver [45] based on the finite-element time-domain method.We observe excellent agreement between the numerical and analytical solutions.Then, we study the effect of time modulation on a Huygens' metasurface [48][49][50][51] and obtain insights into the richness of the physics involved in light interaction with spatiotemporal metasurfaces.We find that a Huygens' metasurface no longer exhibits zero backscattering under time modulation.Furthermore, we demonstrate the possibility of negative absorption already reported for time-varying scatterers [30,40,52].The theory developed in our paper can be straightforwardly extended to the case of three-dimensional (3D) time-varying multilayer structures and metamaterials.

Theory
We consider the problem of light scattering from an infinite periodic array of dispersive scatterers arranged in 2D.The electric material properties of each of these scatterers vary periodically in time with the same modulation frequency  m .Figure 1(a) shows an illustration of that scattering problem.First, we recapitulate the electrodynamics of time-varying bulk media.Then, we use the discrete translational symmetries of a spatiotemporal metasurface to introduce the representations for the incident and scattered fields.We define them in spherical and Cartesian basis for convenience.Then, starting from the T-matrix of an isolated time-varying scatterer [40], we construct the T-matrix of a metasurface in both representations.We use appropriate transformation rules between spherical and Cartesian representations for that purpose.Finally, we give expressions for relevant physically observable quantities.

Time-varying bulk media
Initially, we consider the propagation of an electromagnetic wave in a linear, isotropic, nonmagnetic, source-free, dispersive, and time-varying bulk medium.The constitutive relations that accompany Maxwell's equations inside such a medium, for the electric displacement D and the magnetic flux density B, are written as where  0 is the vacuum permittivity,  0 is the vacuum permeability, E is the electric field, P is the electric polarization density, and H is the magnetic field intensity.The electric polarization density P for such a time-varying medium is written as [53,54] where  e (, −) is the electric response function of the medium, which shall respect causality [40].
We use the Lorentz dispersion model to describe the electric response.A superposition of such Lorentz oscillators can describe the general electric response of the medium [55].However, in what follows, we consider a single oscillator.Specifically, we assume the electron number density of such a system to be an explicit function of time given by [40] where  0 is the electron density of the time-invariant medium,  s is the modulation strength, and  m is the modulation frequency.According to the Lorentzian model, the electric polarization density obeys the differential equation [40,56,57] where   and   are the damping parameter and resonance frequency of the Lorentz oscillator, respectively, and  and  e are the charge and mass of an electron, respectively.Equation ( 4) can be further solved to calculate the electric response function Re ( −  ′ ,  ′ ) of our system (see Eq. (S2) and [40]).Importantly, note that, in the frequency domain, due to the periodic time-variance of the electric response function of the system, a harmonic electric field at frequency  induces a polychromatic polarization density inside the medium at the harmonics  +  m , with  ∈ Z [40].Hence, Maxwell's equations inside such time-varying media are now coupled in frequency through the electric response.

Representations for the input and output fields of a time-varying metasurface and T-matrix calculations
Here, we want to study the electromagnetic response of a 2D periodic array of scatterers made from a time-varying medium.Such a system is periodic in space and time.Hence, it constitutes a spatiotemporal crystal.First, the spatiotemporal crystal is invariant upon discrete translations in space since it is composed of scatterers placed at Here,  1 ,  2 ∈ Z, and R 1 , R 2 are the two primitive lattice vectors of the 2D Bravais lattice of the metasurface (which we consider to lay in the -plane).Moreover, the spatiotemporal crystal is also invariant upon discrete translations in time since the electric response function of the scatterers is periodic in time, with a period of  m = 2/ m .Therefore, a time-varying metasurface comprises a three-dimensional (3D) spatiotemporal crystal.According to Floquet theory, our system is conveniently analyzed in terms of eigenstates of the respective discrete translation operators by introducing the temporal Floquet frequency Ω and the spatial Floquet frequency vector G (that also lays in the -plane).By letting those Floquet frequencies take values within the first Brillouin zone (BZ) of the reciprocal lattice of our spatiotemporal crystal [see Fig.
where BZ 1 denotes integration within the first BZ of the reciprocal lattice of the 2D array of scatterers.With Ẽ(r, ; Ω, G), we denote fields that are eigenstates of the following set of spatiotemporal discrete translation operators: where with  1 ,  2 ∈ Z, and  =  3  m , with  3 ∈ Z.We will represent those fields using two basis sets: a Cartesian and a spherical.That is to say, we will further expand those fields into a superposition of plane and spherical waves, respectively.The plane wave (PW) expansion is conveniently used to represent the input and output fields of our time-varying metasurface, as it constitutes a planar-like structure [43].However, the spherical wave (SW) expansion needs to be used to solve the scattering problem of such an array of scatterers based on the prior knowledge of the solution to the scattering problem of an individual scatterer.In what follows, we will use the superscripts "(c)" and "(s)" to refer to such Cartesian and spherical expansions of the fields, respectively.Specifically, we represent the incident field in the two bases as [58] Ẽinc (r, ; where we introduced the complex Cartesian and spherical incident amplitudes,  inc, (c) g (Ω, G) and  inc, (s)  (Ω, G; r 0 ), respectively.On the one hand, Ω  G g   represents a monochromatic PW with frequency , and We consider the square root to have a positive real part.Here, Δ =  ↓   ≥0 +  ↑  <0 (1 − Δ 0 ) +  ↓ Δ 0 , with Δ 0 taking the value of one, when  2  ∈ R and  2  < G g 2 , and zero otherwise.We use the index  to refer to the direction of propagation of the PW, taking the values ↑ (↓) to refer to propagation/decay towards the upper(lower) half-space, respectively.We introduce Δ to guarantee propagation/decay along the direction defined by the index  for all frequencies.Further,   is the Kronecker delta,   = (Ω  )Ω  / 0 is the wave number of the PW, g , = g 1 + g 2 , with ,  ∈ Z, and g 1 , g 2 are the two primitive reciprocal lattice vectors [43].The summation over g implies a summation over the indices , .Furthermore, (Ω  ) is the (generally dispersive) refractive index of the host medium that is supposed to be isotropic, and  0 is the speed of light in free space.Moreover, the index  takes the values M(N) to refer to TE(TM) polarization.
On the other hand, Ω r 0 represents a monochromatic vector spherical harmonic (VSH), that is either regular ( = 1) or radiating ( = 3), with total angular momentum  = 1, 2, . . ., angular momentum along the -axis  = −, . . ., , and frequency Ω  .The VSH is centered at an arbitrary position r = r 0 .Due to the discrete translation symmetry in space, the spherical incident amplitudes have the symmetry  inc, (s)  (Ω, G; r 0 + R) =  G•R  inc, (s)  (Ω, G; r 0 ).An incident field that has no singularities within the slab that bounds the array of scatterers at the planes  =  − and  =  + , is conveniently represented inside that slab through Eq. ( 8).Furthermore, from the Cartesian incident amplitudes  inc, (c) g (Ω, G), we can analytically get the spherical incident amplitudes  inc, (s)  (Ω, G; r 0 ).We refer to the supplementary material for the definitions and the spatiotemporal representations of Ω  G g   and Ω     (1) r 0 and for the aforementioned transformations between the two equivalent representations of the incident field.
The scattered field is similarly represented by , where we introduced the complex Cartesian and spherical scattered amplitudes,  sca, (c) g (Ω, G) and  sca, (s)  (Ω, G) respectively.Here, Ω  G    (3) represents an array of radiating VSHs distributed over the Bravais lattice, i.e., the origins of the scatterers: The summation over R implies a summation over the integers  1 ,  2 .From the spherical scattered amplitudes  sca, (s)  (Ω, G) we can analytically get the Cartesian scattered amplitudes  sca, (c) g (Ω, G) (see the supplementary material).Conclusively, by employing vectors that contain the incident and scattered amplitudes and by appropriately truncating the infinite sums over the indices of the eigenstates, we can represent the (linear) scattering system of a periodically time-varying, 2D array of scatterers, either in the Cartesian or the spherical representation, by where we introduced the Cartesian and the spherical T-matrices, T(c) (Ω, G), and T(s) (Ω, G).
They represent the scattering system in either of the two bases.The spherical representation of the T-matrix of our system, T(s) (Ω, G), is also known in the literature as the effective T-matrix of the scatterer.It is an effective T-matrix because it could be conceived as a renormalization of the T-matrix of the individual scatterer once it is placed inside the lattice (see Fig. 2).
We denote the T-matrix of an individual time-varying scatterer, placed at r = r 0 , as T(s) 0 (Ω; r 0 ).Its elements are given by (3)  r 0 Ω     T(s) 0 (Ω; r 0 ) Ω  ′  ′  ′  ′ (1) r 0 .They connect the scattered fields with the incident fields of the individual scatterer, once they are represented in the basis of a series of VSHs centered at r = r 0 , i.e., at the origin of the individual scatterer.In principle, T(s) 0 (Ω; r 0 ) can be calculated numerically for an individual scatterer of arbitrary geometry by generalizing an existing numerical method for static scatterers to the time-varying case [59].However, this, unfortunately, generally requires a rather big computational effort.Nevertheless, in Ref. [40], the T-matrix of an individual, homogeneous time-varying sphere is semi-analytically calculated.It can be conveniently used here to study the optical response of an array of time-varying spheres, even though our theoretical approach is generally valid for arbitrary scatterers.
In Fig. 2, we illustrate the different T-matrices representing the systems of static/time-varying isolated spheres/array of spheres in the spherical basis.The symmetries of each of the four scattering systems induce specific symmetry-protected zeros in the entries of the T-matrix.Specifically, on one side, time-invariance enforces zeros in the entries of the T-matrix that correspond to multipolar transitions between different frequencies, and, on the other side, the rotational symmetry of an isolated sphere enforces zeros in the entries of the T-matrix that correspond to multipolar transitions between multipoles with different multipolar indices.Note that the symmetries of the lattice in the cases of an array of spheres introduce extra symmetry-protected zeros in the elements of the T-matrix.They correspond to prohibited multipolar transitions between multipoles that belong to different irreducible representations of the symmetry group of the array [60].We neglect this for simplicity in our illustration.
As a next step, we can solve for T(s) (Ω, G) in Eq. ( 14) by writing [43] A sca, (s) (Ω, G) = T(s) 0 (Ω; R) This equation physically states that each scatterer of the metasurface, once excited by a field that has the same discrete translational symmetry as that of the lattice, 1) sees as an effective incident field the external incident field plus the scattered field from all the other scatterers in the 2D array, and 2) scatters the same field, upon a difference of phase, according to its individual T-matrix, T(s) 0 (Ω; R).The matrix Ĉ(3) (−R ′ ) represents the translation of radiating VSHs centered at R+R ′ , into a series of regular VSHs centered at r = R, Ω R (see the supplementary material).Combining Eqs. ( 14) and ( 15), we readily get the expression for the effective T-matrix of our system in the spherical representation: where Û is the identity matrix.The Ewald method can efficiently evaluate the infinite summation of the translation matrices over the lattice in real space in the above equation (see the supplementary material and [46]).Note that for the efficient calculation of the spherical T-matrix in Eq. ( 16), one should generally consider representations with a different number of multipoles for each frequency of a spectral comb characterized by some Floquet frequency Ω.Higher frequencies require a larger number of multipoles since they correspond to interactions with optically larger scatterers.
As a next step, we can analytically transform the T-matrices from the spherical to the Cartesian representation and get T(c) (Ω, G).This requires the transformation of the input and the output vectors of the T-matrix: where analytical expressions for the input and output transformation matrices, Î(Ω, G; r 0 ), Ô(Ω, G), are given in the supplementary material.We readily get the following expression for the Cartesian T-matrix: Finally, we can write the Cartesian S-matrix, that connects the Cartesian amplitudes of the outgoing waves,  out, (c) g (Ω, G) =  inc, (c) g (Ω, G) +  sca, (c) g (Ω, G), in either of the two half-spaces of the hosting medium, with the Cartesian amplitudes of the incoming waves,  in, (c) g (Ω, G) =  inc, (c) g (Ω, G), as follows: Let us note that an extension to multilayer time-varying metasurfaces, based on the S-matrices of each layer, is also straightforward (see the supplementary material and [43]).

Physically observable quantities
We need to define some physically observable quantities to study the response of time-varying metasurfaces.We assume that our metasurface is embedded inside a lossless host medium.The total energy flux towards the positive -axis of an electromagnetic field (represented in the Cartesian basis by some amplitudes  (c) g (Ω, G) and which may contain evanescent waves decaying exclusively along a single direction) propagating in the host medium is given by: where H(r, ) is the magnetic field, and  g (Ω, G) denotes the energy spectral density flux that a PW, Ω  G g   , carries along the positive -axis, which is given by the formula: Here,   is the wave impedance of the host medium at the frequency Ω  .Assuming an excitation of the time-varying metasurface with a single monochromatic PW, Ω  G g   (with  ≥ 0 in what follows), we can define the physical observables of reflectivity [ℛ g (Ω, G)], transmissivity [ g (Ω, G)], and absorptivity [ g (Ω, G)] as: where d ≠ .We can see that those observables are directly related to the elements of the S-matrix of the metasurface.

Comparing analytical solutions with full-wave numerical simulations
We verify at first the theory developed in the last section by comparing the near-field calculations of a full-wave numerical solver [45] with our analytical solutions.These full-wave numerical simulations are based on the finite-element time-domain method.We consider a time-varying metasurface composed of a square lattice of dispersive spheres embedded in free space [see Fig. 1(a)].We make use of the Lorentz oscillator model of Eq. ( 4), and in what follows, for generality, we will consider the resonance frequency   as a free parameter with respect to which we define all the relevant quantities.Specifically, the damping parameter   is taken as   /8, and the electron number density  0 is set to 11 2      0 / 2 .Furthermore, the radius of the spheres is chosen as  = 3.77  0 /  , and the spatial periodicity in both  and  directions is chosen as |R 1 | = |R 2 | = 5.Moreover, the frequency of the modulation of the electron number density is set to  m = 0.1  , while its strength is taken as  s = 0.9 [see Eq. ( 3)].Next, we set the incident field to be a Gaussian pulse of an x-polarized PW given by where  0 = 1 V/m,  0 = 2.9 × 2/  ,  0 = 0.31  ,  0 = 8 0 , and  0 is the speed of light in free space.We provide in the supplementary material the incident amplitudes that expand such excitation in the Cartesian basis, which we use to calculate the input vector in our algorithm.
To set up these simulations in the full-wave reference solver, we exploit the mirror symmetries of our system about the -and -planes [see Fig. 1(a)] to reduce the simulation domain.We appropriately use perfect electric conductor (PEC) and perfect magnetic conductor (PMC) boundary conditions to mimic a periodic system.The mesh size inside the sphere was taken as 0.7 0 /  , while in the host medium, it was taken as 1.7 0 /  .
We compute the near fields at two spatial points (above and below the metasurface).These points have coordinates P(1.25, 1.25, 1.2) and Q(1.25, 1.25, −1.2).The comparison of near fields is shown in Fig. 3.We observe an excellent agreement between our analytical solutions with the fields calculated using the full wave numerical solver.Besides the excellent accuracy, it is worth mentioning that the numerical full-wave solver takes around 68 hours for the computation.In contrast, our analytical solutions could be obtained within only 0.128 hours, using multipoles up to the third order and a truncation of the spectrum at 1.5  .These calculations were performed on a computing node of a cluster.

The effect of time modulation on Huygens' metasurfaces
Huygens' metasurfaces form an important class of metamaterials.They offer near-zero backscattering at specific spectral regions.Therefore, to highlight the opportunities of what can be explored with the present code, we study the impact of a time variation of the material properties on the optical response of such a Huygens' metasurface, i.e., in essence, we discuss a time-varying Huygens' metasurface.
As in the previous section, we consider a square lattice of spherical scatterers with the same Lorentz dispersion model.Initially, we assume the time-invariant case (i.e.,  s = 0) and find a specific geometry for which our metasurface shows near-zero backscattering at some spectral region.We consider the normal incidence of a single TE-polarized, upwards-propagating, monochromatic PW, whose frequency,  inc , we vary.The incident PW is the same as shown in Fig. 1(a).The Huygens' condition is achieved for the radius  = 4.27 0 /  and the spatial periodicities |R 1 | = |R 2 | = 11.930 /  .The cyan curve in Fig. 4(a) shows the reflectivity (ℛ) for the time-invariant case.We observe the existence of near-zero backscattering at two spectral regions centered at  inc /  = 0.18 and 0.32.
Next, we introduce the temporal modulation of the metasurface.The modulation frequency is fixed at  m = 0.1  .We consider two different modulation strengths,  s = 0.5 and 0.9, and observe the change in reflectivity ℛ.Note that ℛ includes contributions from all the reflected frequencies [see Fig. 1(a) and Eq. ( 23)].We find that the near-zero backscattering gets spoiled by the temporal modulation of the array [see Fig. 4(a)].This is expected since the T-matrix of our metasurface undergoes a renormalization process that is induced by the temporal modulation [see Fig. 2(b) and (d)].As expected, we observe a stronger deviation in the reflected spectra with increasing modulation strength.In fact, we observe a blue shift in the resonances of ℛ.These observations also hold for  [see Fig. 4(b)].These effects are mainly observed in the subwavelength regime, where the material losses of the Lorentz dispersion model are low.
To explain these observations, we look at the electric (ED) and magnetic dipolar (MD) contributions to the scattered field.For simplicity, we only show the plots of the multipolar spectra for the zeroth spectral diffraction order, i.e., for   ′ =   =  inc .However, as we will see later, a significant spectral coupling is also happening to other spectral diffraction orders.We plot in Figs.4(d) and (e) the amplitude ratios and the phase differences of these multipolar contributions.We denote them as  MD =  sca, (s) ,1,−1,M (Ω, 0) and  ED =  sca, (s) ,1,−1,N (Ω, 0), respectively.It is well known that to achieve the Huygens' condition (1st Kerker condition) for a metasurface operating in the dipolar regime (i.e., with optically small scatterers), we need |  MD |/|  ED | = 1 and ∠  MD − ∠  ED = − (see supplementary material of [52] and [50]).We observe that those two conditions are simultaneously satisfied at the expected frequency of the first Huygens' condition for the time-invariant case.They seem not to be satisfied for the case of the second Huygens' condition at the larger frequency, but this is just because the metasurface ceases to be purely dipolar in that spectral region and, hence, an other, generalized Kerker condition applies [61].Indeed, the multipolar conditions for zero-backscattering are still satisfied for that case if we consider all the multipolar contributions.As expected, the temporal modulation of the metasurface alters the multipolar content of the radiating fields in amplitude and phase, spoiling the Huygens' condition and destroying the zero-backscattering effect.
It is important to note that for the time-modulated metasurface, satisfying the Huygens' condition only at the zeroth spectral diffraction order is insufficient to achieve zero reflection.The reason for this is the existence of finite contributions to the reflected energy from other spectral diffraction orders [see Fig. 1(a) and Eq. ( 23)].In Fig. 5, we analyze the contributions of the different PWs that take part in the spectra of the observables previously discussed in Fig. 4, for the particular case of  s = 0.9.There, the complexity of the optical system of a time-varying metasurface can be appreciated in view of its S-matrix.In Fig. 5, for simplicity, we only show the absolute values of the S-matrix elements associated with up to the first spectral and spatial diffraction orders.Higher spatiotemporal orders have vanishing contributions to the reflected and transmitted fields.Note that the metasurface has mirror symmetries about the -, -, and -planes.Those mirror symmetries of the metasurface introduce symmetries between elements of the S-matrix that correspond to symmetric spatial diffraction orders.Therefore, we omit to plot the spectra of those symmetric spatial diffraction orders.See the supplementary material for more details about the effect of the above mentioned mirror symmetries of the system on its S-matrix.For example, such a symmetry analysis predicts that, upon normal incidence, due to the symmetries of the system concerning the -and -planes, there cannot be reflection or transmission in the zeroth spatial diffraction order of the opposite polarization compared to that of the excitation.We can also observe this symmetry selection rule in the spectra in Figs.5(d) and (j).We can observe in Figs.5(a) and (g) the rather significant contributions to ℛ and  from the non-zero spectral diffraction orders (  ′ =  ± 1) in the zeroth spatial diffraction order (g ′ 00 ).This spectral coupling effect primarily spoils the Huygens' condition in the time-varying case.Also note the solid blue lines in the second column of Fig. 5 extending into the light-blue-shaded spectral region in Fig. 4. As expected, this signifies the presence of propagating spatial diffraction orders for the first blue shifted spectral diffraction order (  ′ =  + 1) in the subwavelength spectral region of the time-invariant metasurface.Therefore, time modulation allows the scattered field to get coupled to the spatial diffraction orders even though the metasurface is sub-wavelength for the incident field.
Last but not least, another important feature to be noted is the onset of negative absorption for the time-modulated case [see Fig. 4(c)].According to Noether's theorem, time-invariance symmetry is responsible for energy conservation within a system.Therefore, in time-modulated media, the photons can exchange energy with the externally modulated matter.This explains the observation of negative absorption in our spatiotemporal metasurface.In fact, we can adequately adjust the parameters of such time-varying systems to observe parametric amplification effects [52].

Conclusions
We have presented an approach to semi-analytically compute the scattering response of a spatiotemporal metasurface.We started with the dynamics of a time-varying bulk media having Lorentzian dispersion.Then, we constructed the effective T-matrix of a spatiotemporal metasurface using the T-matrix of a time-varying sphere.We emphasize that our developed multiple-scattering theoretical framework is quite general and does not depend on the shape of the constituent scattering object.Then, we verified our theory against full-wave numerical simulations, demonstrating an excellent agreement.Finally, we studied the effect of time modulation on Huygens' metasurfaces and explained it in terms of spatiotemporal diffraction orders.
As a next step, our multiple scattering approach can be extended to the 3D metamaterials to study the complex band structure of such four-dimensional (4D) spatiotemporal crystals.Moreover, using our semi-analytical formalism, it is interesting to study parametric amplification effects in such spatiotemporal metasurfaces.Finally, another exciting extension of the current research can be the study of homogenization techniques in spatiotemporal metamaterials.
Funding.Content in the funding section will be generated entirely from details submitted to Prism.Disclosures.The authors declare no conflicts of interest.

Fig. 1 .
Fig. 1. a) Illustration of light scattering from a spatiotemporal metasurface made of periodically time-varying scatterers (with frequency  m ) arranged on a 2D lattice in the -plane.The excitation of a spatiotemporal metasurface with a monochromatic plane wave induces a response that is comprised of a rich spectrum of diffraction orders both in the spectral and the momentum dimensions.b) Illustration of the reciprocal lattice of such a spatiotemporal metasurface, which constitutes a 3D crystal with two spatial and one temporal dimension, i.e., demonstrating discrete translation symmetries along two directions in space and, also, in time.The reciprocal lattice resides in a 3D space comprised by the  and  components of the linear momentum in the two dimensions (  ,   ), and the frequency of light in the third dimension ().The shaded magenta box represents the first Brillouin zone of the 3D crystal that is spanned by the Floquet frequency vectors (denoted with yellow color) characterizing the reciprocal lattice with some temporal Floquet frequency Ω and spatial Floquet frequency vector G.
1(b)], we can generally express the involved fields as a superposition of fields that are eigenstates of the respective discrete translation operators in space [ Ts (R)] and time [ Tt ()] by

Fig. 2 .
Fig. 2. Illustration of T-matrices in the spherical basis for different scattering systems: a) time-invariant sphere, b) time-invariant metasurface, c) time-varying sphere, d) time-varying metasurface.Here, the white shaded entries represent the symmetryprotected zeros of the elements of the T-matrices, Ω  denotes the  th frequency of a given spectral comb, and   denotes the  th spherical multipole characterized by , , .The rotational symmetry of the isolated sphere allows only for invariant transitions with invariant multipolar indices, whereas the time-variance introduces coupling among multipoles of different frequencies.Note that, in general, there are extra lattice-symmetry-induced-zeros in the elements of the T-matrices of the arrays, which correspond to particular multipolar transitions prohibited by the symmetries of the lattice.Also, note that the number of multipoles for each frequency shall generally vary: typically, a larger number of multipoles for higher frequencies is needed.

Fig. 3 .
Fig. 3. Comparison between analytically and numerically obtained near field signals, at point P (transmission): a) in the frequency domain and b) in the time domain; and at point Q (reflection): c) in the frequency domain and d) in the time domain.The light green shaded region represents the frequency (time) interval where 99% of the energy of the incident pulse resides in the frequency (time) domain.All specific details on the scattering setting are described in the main text.

Fig. 4 .
Fig. 4. The effect of time modulation on a Huygens' metasurface.Plots of a) reflectivity [ℛ 0↑M (Ω, 0)], b) transmissivity [ 0↑M (Ω, 0)], and c) absorptivity [ 0↑M (Ω, 0)] for different modulation strengths,  s , as a function of the excitation frequency,  inc = Ω  , of a normally incident, TE-polarized, PW.The modulation frequency considered is  m /  = 0.1.d) Ratio of amplitudes, and e) phase difference of magnetic and electric dipolar contributions for the zeroth spectral diffraction order.Here, we use the following notation  MD =  sca, (s) ,1,−1,M (Ω, 0) , and  ED =  sca, (s) ,1,−1,N (Ω, 0).Two Huygens' conditions with zero backscattering are observed at  inc /  = 0.18, 0.32 for the case of unmodulated metasurface (denoted with the two black vertical lines).In the lower frequency, where the response of the metasurface is mainly dipolar, the 1st Kerker condition of |  ED | = |  MD |, ∠  MD − ∠  ED = − is satisfied.Note that at low excitation frequencies there exists a spectral region with negative absorptivity in (c).The light-blue-shaded region represents the frequency interval up to which the lattice is subwavelength concerning the wavelength of the zeroth spectral diffraction order.

Fig. 5 .
Fig. 5. Plots of the absolute value of the elements of the S-matrix,Ŝ  ′ g ′  ′  ′ ;(c)0↑M(Ω, 0) , that are related with the spectra of the metasurface presented in Fig.4for the case of  s = 0.9, as a function of the incident frequency  inc = Ω  , (a)-(f): for reflection ( ′ =↓), and (g)-(l): for transmission ( ′ =↑).Here, the solid curves correspond to propagating, and the dotted curves correspond to evanescent waves.The color of the lines denotes the different spectral diffraction orders considered.There are symmetries between the spectra of different spatial diffraction orders induced by the symmetries of the square lattice (see text).
Acknowledgments.P.G., A.G.L., and C.R. are part of the Max Planck School of Photonics, supported by the Bundesministerium für Bildung und Forschung, the Max Planck Society, and the Fraunhofer Society.P.G. and A.G.L. acknowledge support from the Karlsruhe School of Optics and Photonics (KSOP).D.B. and C.R. acknowledge support by the German Research Foundation through Germany's Excellence Strategy via the Excellence Cluster 3D Matter Made to Order (EXC-2082/1 -390761711) and the Carl Zeiss Foundation through the "Carl-Zeiss-Focus@HEiKA".B.V. and C.R. acknowledge support by the German Research Foundation within the SFB 1173 (Project-ID No. 258734477).