MANUSCRIPTDRAFT 1 Floquet-Mode Solutions of Space-Time Modulated Huygens’ Metasurfaces Shulabh Gupta, Tom. J. Smy and Scott A. Stewart Abstract—A rigorous Floquet mode analysis is proposed for exoticeffectssuchasgenerationofnewharmoniccomponents a zero thickness space-time modulated Huygens’ metasurface and Lorentz non-reciprocity [9][10]. Space-time modulated to model and determine the strengths of the new harmonic metasurfaces fall within the general framework of space-time componentsofthescatteredfields.Theproposedmethodisbased modulated mediums [11][12], which have found important on Generalized Sheet Transition Conditions (GSTCs) treating a metasurface as a spatial discontinuity. The metasurface is applications in acousto-optical systems for spectrum analysis, described in terms of Lorentzian electric and magnetic surface parametricoscillatorsandamplifiers,forinstance[13][14][15]. 7 susceptibilities, χee and χmm, respectively, and its resonant fre- Consequently, a combination of the wave-shaping capabil- 1 quencies are periodically modulated in both space and time. ities of Huygens’ metasurfaces with space-time modulation 0 The unknown scattered fields are then expressed in terms of principles, is an interesting avenue to explore for advanced 2 Floquet modes, which when used with the GSTCs, lead to a system of field matrix equations. The resulting set of linear electromagnetic wave control, in both space and time. To n equations are then solved numerically to determine the total investigate into the properties of space-time modulated Huy- a J scattered fields. Using a finite-difference time domain (FDTD) gens’ metasurfaces, a finite-difference time-domain (FDTD) 9 solver, the proposed method is validated and confirmed for technique has recently been proposed to analyze a zero several examples of modulation depths (∆ ) and frequencies 1 p thickness model of Huygens’ metasurfaces [16][17], based (ω ). Finally, the computed steady-state scattered fields are Fopurier propagated analytically, for visualization of refracted on Generalized Sheet Transition Conditions (GSTCs) [18]. ] s harmonics.Theproposedmethodissimpleandversatileandable In contrast to the other existing techniques for analyzing c todeterminethesteady-stateresponseofaspace-timemodulated static metasurfaces in frequency domain [19][20], the FDTD i Huygen’smetasurface,forarbitrarymodulationfrequenciesand t analysis is naturally applicable to the problem of space-time p depths. modulated metasurfaces, considered herein. o Index Terms—Electromagnetic Metasurfaces, Electromagnetic While such numerical methods are useful to determine the . s Propagation,FloquetAnalysis,ExplicitFinite-Difference,Gener- time-evolution of scattered waves from space-time modulated c alizedSheetTransitionConditions(GSTCs),LorentzDispersions, i metasurfaces for a given input wave, efficient determination s Parametric Systems. y of steady-state response of the metasurface is an equally h important problem to solve. This issue is addressed in this p I. INTRODUCTION work, whereby exploiting the periodic nature of the spatio- [ Recently, there has been a strong growing interest in temporal perturbation on the metasurface, the scattered fields 1 Huygen’s metasurfaces due to their impedance matching ca- areexpressedintermsofFloquetmodes.TheHuygens’meta- v pabilities with free-space and their versatile applications in surface is modelled using surface susceptibilities following a 1 wavefront shaping [1][2][3]. They are constructed using a physically motivated Lorentzian profile, whose resonant fre- 7 2-D array of electrically small Huygen’s sources, exhibiting quenciesareparametrizedtoemulateaspace-timemodulation 2 perfect cancellation of backscattered fields, due to optimal of the metasurface. Combined with GSTCs, the Floquet mode 5 0 interactionsoftheirelectricandmagneticdipolarmoments[4]. amplitudes are computed by solving a set of linear equations. . Some efficient implementations of Huygens’ metasurfaces are The proposed method thus efficiently computes the steady- 1 0 based on all-dielectric resonators [5][1][6] and orthogonally state response of a zero-thickness space-time modulated Huy- 7 collocated small electric and magnetic dipoles [7][8]. gens’ metasurface, consistent with the FDTD field solutions. 1 While the majority of the work on metasurfaces has been Furthermore, integrating the proposed method with analytical v: focussed on static (linear time invariant) metasurfaces, there Fourier methods [13], the propagation of scattered fields are i is also a growing interest in dynamic metasurfaces, where the conveniently visualized in free-space. X constitutive parameters of the metasurface unit cells are con- The paper is structured as follows. Section II describes ar trolled in real-time. An important class of such dynamic elec- the problem statement of this work, and develops the fields tromagnetic structures is space-time modulated metasurfaces, equations governing the scattering fields from a space-time wheretheirconstitutiveparametersareperiodicallymodulated modulated metasurface based on GSTCs and Lorentz surface in both space and time, at comparable frequency scales to susceptibilities. Section III presents the proposed method the input excitations. This leads to a complex interaction of based on Floquet mode expansions, forming the set of lin- the incident wavefronts with the metasurfaces resulting in ear equations to be solved numerically. Several results are then presented for both time-only and space-time modulated S. Gupta, T. J. Smy and S. A. Stewart are with the Department of metasurfaces.ConclusionsareprovidedinSec.IV,andabrief Electronics, Carleton University, Ottawa, Ontario, Canada. Email: shu- [email protected] summaryoftheFDTDmethodisprovidedinAppendixV,for MANUSCRIPTDRAFT 2 the sake of self-consistency and completeness of the paper. x t SPACE-TIME MODULATED METASURFACE II. SPACE-TIMEMODULATEDMETASURFACES A. Problem Statement χee(x,t),χmm(x,t) ψ2(r,t)sin{(ω0−ωp)t} Metasurfaces are zero thickness electromagnetic structures ψ(r,t)=ψ0(r,t)sin(ω0t) that act as discontinuities in space. The exact zero thickness natureofelectromagneticstructureswasdevelopedbyIdemen ψ1(r,t)sin{ω0t} t in terms of Generalized Sheet Transition Conditions (GSTCs) t [21], which were later applied to metasurfaces [18]. Consider a Huygens’ metasurface illustrated in Fig. 1, where the meta- ψ3(r,t)sin{(ω0+ωp)t} surfaceliesinthex y planeatz =0,withawaveincidence z on the left, normal t−o the surface. The wave interacts with the z=0 metasurface and produces a transmitted and a reflected wave, t along the forward and backward direction, respectively. This interaction of the metasurface with the electromagnetic waves Fig.1. AgeneralIllustrationofaspace-timemodulatedHuygens’metasur- faceundernormallyincidentCWplane-waveresultingingenerationofseveral is described by the GSTCs, using the transverse components frequency harmonics, refracted along different angles. ωp is the pumping of electric and magnetic surface polarizabilities P and M, as frequency. [22] dP (x,t) zˆ ∆H(x,t)= || (1a) where ωn = ω0 +nωp, and βn = βx0 +nβp with βx0 = 0 × dt due to assumed normal input incidence. Each harmonic term dM (x,t) ofthisexpansionwithatemporalfrequencyω ,representsan ∆E(x,t)×zˆ=µ0 |d|t , (1b) oblique forward propagating plane-wave in thne +z direction, where and making an angle θn measured from the normal of the metasurface, as illustrated in Fig. 1. These refraction angles are given by ∆E=(E E E ), ∆H=(H H H ). t− 0− r t− 0− r (cid:20) n β (cid:21) The surface polarizabilities on the Huygens’ metasurface are θ(ωn)=sin−1 p , (4) (1+nω /ω )k p 0 0 related to the average fields around the metasurface and can where k is the free-space wavenumber of the fundamental be described in terms of scalar electric and magnetic surface 0 susceptibilities χ and χ , respectively, as1. frequency. Based on this simple physical argument, the angle ee mm of refraction of these newly generated frequency components Q˜ (ω)=χ˜ E˜ (ω), (2a) fromaspace-timemodulatedmetasurface,canbedetermined. ee av || However, the strengths of these harmonics (weights qn) are M˜ (ω)=χ˜ H˜ (ω), (2b) still unknown and must be determined taking the exact elec- mm av || tromagneticinteractionofthewaveswiththemetasurfaceinto account using the GSTCs of (1). (cid:34) (cid:35) (cid:34) (cid:35) E˜ +E˜ +E˜ H˜ +H˜ +H˜ where E˜ = 0 t r , H˜ = 0 t r av 2 av 2 B. Field Equations The space-time modulated problem considered here is nat- and Q=P/(cid:15) is the normalized electric polarizability2. 0 urally treated in the time domain. Consequently, the surface Next, consider a space-time modulated metasurface, whose susceptibilities must be defined in time domain as well. The electric and magnetic susceptibilities are both a function of most common, and causal, description of the susceptibilities space and time, i.e. χ (x,t) and χ (x,t). Let us restrict ee mm is in terms of Lorentz dispersion, which is also typical of heretoaperiodicmodulationonly,withapumpingfrequency the Huygens’ sources used to construct the metasurfaces ω and the spatial frequency β . Due to this periodic spatio- p p [5][23][24]. Consequently, the electric and magnetic suscepti- temporal perturbation on the metasurface, assumed to be blities can be expressed as infiniteinsize,thescatteredtransmission(andreflection)fields can be expressed in terms of Floquet series as ω2 χ˜ (ω)= ep (5a) ee (ω2 ω2)+iα ω (cid:88)∞ e0− e Et(x,z =0+,t)= qnejωntejβnx, (3) ω2 n=0 χ˜mm(ω)= (ω2 ωm2)p+iα ω, (5b) m0− m 1Assumingnorotationofpolarizationandonlyonecomponentofthefields where (ω , ω ), (ω , ω ), (α ,α ) are the electric and alongtheprincipalaxis. e0 m0 ep mp e m 2The fields expressions with a tilde, ψ˜ denote the frequency domain magnetic resonant frequencies, plasma frequencies and loss quantities. coefficients, respectively. MANUSCRIPTDRAFT 3 The space-time modulation can now be introduced into III. PROPOSEDFLOQUETMODEEXPANSIONSOLUTIONS the metasurface, by sinusoidally modulating, for instance, the A. Field Matrix Equations resonant frequencies of the two Lorentzian susceptibilities, as Let us consider a metasurface with a periodic space-time [17] modulation of the resonant frequencies following (6), and expand the unknown reflected and transmitted fields using ω (x,t)=ω 1+∆ cos(ω t+β x) (6a) Floquet expansion as e0 e0 e p p { } (cid:88) (cid:88) ωm0(x,t)=ωm0{1+∆mcos(ωpt+βpx)}, (6b) Et(t)= qnejω0tejn(ωpt+βpx) = qnejω0tejnΩ, (12a) n n where ∆ and ∆ are the modulation depths and the spatial (cid:88) (cid:88) e m E (t)= p ejω0tejn(ωpt+βpx) = p ejω0tejnΩ, (12b) frequencies of the perturbation, for electric and magnetic r n n n n resonances, respectively. They can now be inserted into the where a new variable Ω = (ω t + β x) is introduced for GSTCs of (1), for a prescribed input fields, to determine the p p compact notation. Similarly, the corresponding unknown po- total scattered fields. larizabilities, can also be expanded in Floquet series as To illustrate the procedure, let us assume a normally in- cident plane-wave (E = const.) where the corresponding | 0| (cid:88) (cid:88) transmitted and reflected fields, given by Q (t)= a ejω0tejnΩ, M (t)= d ejω0tejnΩ, 0 n 0 n n n (cid:88) (cid:88) Q (t)= b ejω0tejnΩ, M (t)= e ejω0tejnΩ, ˆz E (z,t) t n t n E0(z,t)=E0ej(ω0t−k0z) yˆ, H0(z,t)= × 0 (7) n n η (cid:88) (cid:88) 0 Q (t)= c ejω0tejnΩ, M (t)= f ejω0tejnΩ. ˆz E (z,t) r n r n Et(z,t)=Et(x,t)e−jk0z yˆ, Ht(z,t)= × t n n η 0 Next, substituting the above expressions in (8) and (9) and E (z,t) ˆz E (z,t)=E (x,t)e+jk0z yˆ, H (z,t)= r × . re-arranging the terms, for incident, transmitted and reflected r r r η0 E- and H-fields, we get (10), where Each of these field components are related to their respective ∆2ω2 ∆ ω2 polarizabilities on the metasurface (at z = 0) through the ∆ = e e0, ∆ = e e0, Lorentz susceptibilities of (5) with (6), which can then be e1 4ωe2p e2 ωe2p expressed in the time domain as3 [16] (cid:26)ω2 (ω +nω )2 ∆2ω2 (cid:27) A = e0 0 p + e e0 . e,n ω2 − ω2 2ω2 ep ep ep dd2tQ2i +ωe20{1+∆ecos(ωpt+βpx)}2Qi =ωe2pEi(t) (8a) ∆m1 = ∆42mωωm2m2p0, ∆m2 = ∆ωmm2ωpm20, (cid:26)ω2 (ω +nω )2 ∆2 ω2 (cid:27) A = m0 0 p + m m0 d2M ω2 m,n ω2 − ω2 2ω2 i +ω2 1+∆ cos(ω t+β x) 2M = epE (t),4 mp mp mp dt2 m0{ m p p } i η i B =j(ω +jnω ). (13) 0 n 0 p (8b) Theseequationsareobtainedbyexpressingthecosinefunction where the subscript i = 0,t,r for incident, transmitted and using Euler’s form and then grouping the terms of common reflectedfields,respectively.Finally,allthescatteredfieldsare complex exponentials. Each of these series equations, can be related to their polarizabilities following the GSTCs of (1), as truncated to 2N + 1 harmonic terms5, and be written in a matrix form as shown in (11). In compact form, the matrix dQ dQ dQ 2 0 + t + r = (E E E ), (9a) equation (11) can be written as 0 r t dt dt dt η (cid:15) − − 0 0 [C ][V ]=[E ], (14) n n n dM dM dM 2 0 + t + r = (E E E ) (9b) t r 0 which represents 8 (2N + 1) linear equations with same dt dt dt µ − − 0 × number of unknowns. The sought field solutions can now be Equations (8) and (9) thus represent two sets of field numerically computed using this matrix equation as equations, that must be solved to determine the transmitted andreflectedfields,Et(x,z,t)andEr(x,z,t),insteadystate, [Vn]=inv{[Cn]}[En], (15) for a given input excitation E0 = E0ejω0t yˆ at the input of from which the corresponding p ’s and q ’s can be extracted n n the metasurface at z =0 . to construct the reflected and transmitted fields, respectively, − following (12). 3forthelosslesscase,forsimplicity. 4withpropersignoftherighthandsideterm,dependingonincident(−), 5assumingthattheharmonicamplitudesrapidlyfallstozerowithincreasing reflected(+)ortransmittedfields(−),accordingto(7). indexn,whichisthecaseinalltheforthcomingresults. MANUSCRIPTDRAFT 4 (cid:88) [∆ a +∆ a +A a +∆ a +∆ a ]ejnΩt =E (Incident E-Fields) (10) e1 n 2 e2 n 1 e,n n e2 n+1 e1 n+2 0 − − n (cid:88) [∆ b +∆ b +A b +∆ b +∆ b q ]ejnΩt =0 (Transmitted E-Fields) e1 n 2 e2 n 1 e,n n e2 n+1 e1 n+2 n − − − n (cid:88) [∆ c +∆ c +A c +∆ c +∆ c p ]ejnΩt =0 (Reflected E-Fields) e1 n 2 e2 n 1 e,n n e2 n+1 e1 n+2 n − − − n (cid:88)[∆ d +∆ d +A d +∆ d +∆ d ]ejnΩt = E0 (Incident H-Fields) m1 n 2 m2 n 1 m,n n m2 n+1 m1 n+2 − − −η0 n (cid:88) [∆ e +∆ e +A e +∆ e +∆ e b ]ejnΩt =0 (Transmitted H-Fields) m1 n 2 m2 n 1 m,n n m2 n+1 m1 n+2 n − − − n (cid:88) [∆ f +∆ f +A f +∆ f +∆ f f ]ejnΩt =0 (Reflected H-Fields) m1 n 2 m2 n 1 m,n n m2 n+1 m1 n+2 n − − − n (cid:88) [2cB a +2cB b +2cB c +p +q ]ejnΩt =E (GSTC Equation) n n n n n n n n 0 n (cid:20) (cid:21) (cid:88) µ0Bnd + µ0Bne + µ0Bnf +p q ejnΩt = E (GSTC Equation) n n n n n 0 2 2 2 − − n Cn Vn En (cid:122) (cid:125)(cid:124) (cid:123)(cid:122) (cid:125)(cid:124) (cid:123) (cid:122) (cid:125)(cid:124) (cid:123)  Se 0 0 0 0 0 0 0  an   I0   ...   0 Se 0 0 0 0 0 −I  bn   0   0   0 0 Se 0 0 0 −I 0  cn   0   0   00 00 00 S0m S0m 00 00 I/0η0  denn = −I00/η0 ,where I0 = 10   0 0 0 0 0 Sm −I/η0 0  fn   0   0   Sg/(2c0) Sg/(2c0) Sg/(2c0) 0 0 0 I I  pn   I0   .  . 0 0 0 S S S I I q I . g g g n 0 − −     A ∆ ∆ 0 0 A ∆ ∆ 0 0 e, N e2 e1 m, N m2 m1 .− . . . ··.· .− . . . ··.·  . . . . .   . . . . .   . . . . .   . . . . .   ···   ···   ∆e2 Ae, 1 ∆2 ∆e1 0   ∆m2 Am, 1 ∆2 ∆m1 0   − ···   − ···  Se = ∆e1 ∆e2 Ae,0 ∆e2 ∆e1 , Sm = ∆m1 ∆m2 Am,0 ∆m2 ∆m1   ···   ···   0 ∆e1 ∆e2 Ae,+1 ∆e2   0 ∆m1 ∆m2 Am,+1 ∆m2   ···   ···   .. .. .. .. ..   .. .. .. .. ..   . . . . .   . . . . .  ··· ··· 0 0 ∆ ∆ A 0 0 ∆ ∆ A e1 e2 e,+N m1 m2 m,+N ··· ··· S =diag B ,B , ,B ,B ,B , ,B ,B . (11) g N N+1 1 0 1 N 1 N { − − ··· − ··· − } B. Results FouriertransformsandforaclosercomparisonwiththeFDTD method, given by To validate the results of the proposed method, an FDTD solver is used, which has been recently proposed in [16][17] to directly simulate (8) and (9), for an arbitrary time-domain (cid:88)+N (cid:40) (cid:18) t (cid:19)2(cid:41) input. A brief summary of the method is described in the Et(x,t)= qnexp cos(ωnt+βnx), − T 0 appendix for self-consistency of the paper. In the results n= N described next, the modulation depths of the electric and (cid:88)+−N (cid:40) (cid:18) t (cid:19)2(cid:41) E (x,t)= p exp cos(ω t+β x). magnetic resonant frequencies are assumed to be equal for r n n n − T 0 simplicity, i.e. ∆ = ∆ = ∆ . Furthermore, the output n= N e m p − (16) fields from the proposed method are constructed using a Gaussian envelopes6, for better conditioning of the numerical Figure 2 shows several examples of the transmission and reflection spectrum for different modulation depths and fre- 6ThisprovidesafinitebandwidthinthetemporalFFTsasopposedtoideal quencies. While the modulation frequency ω only fixes the p deltafunctionswhicharedifficulttocaptureinnumericalcomputations.These location of the newly generated harmonic components, the Gaussianenvelopesshouldnotbeseenastheinputpulseshapes,butrather aswindowingfunction. modulationdepth∆p controlstherelativestrengthsofthehar- MANUSCRIPTDRAFT 5 monics, compared to the fundamental frequency of excitation. ∆ . For simplicity, and considering the practical point of p An excellent agreement is observed with FDTD solver, for view, a perfectly matched metasurface with χ˜ = χ˜ is ee mm all the cases considered, validating the proposed procedure. assumedexhibitingzeroreflections.Thespectralmapprovides Minor discrepancies are observed between the two, which interesting and useful information about the interaction of the are attributed to several factors. Firstly, the Floquet solution metasurface with the input fields. For example, at ω 2ω , p 0 ≈ assumed zero losses here7, while the FDTD time-domain a substantial amplification of the fundamental frequency ω 0 results consider finite losses. Secondly, the Floquet solution is is observed, which is also the dominant spectral component the steady state response of the metasurface, while the FDTD in the transmitted fields. This is reminiscent of wave ampli- solvertakesintoaccountthedispersionbaseddistortionofthe fications in diverse class of mechanical and electromagnetic input pulse, which is naturally not taken into account in the parametric systems, where specific parameters of the system Floquet analysis in (16). are modulated at twice the excitation frequency leading to wave instabilities [11][12]. ∆p=0.05,fm=0.05f0 0.7 1 FDTD FDTD χee=χmm,∆p=0.1 0.6 Floquet 0.8 Floquet 0 20log|Et(z=0+,ω)| 10 -10 0 10 0.5 E(t)ttF{}00..34 E(t)trF{}00..46 ω0 000...246 --02100 0.5 0.2 0.2 yω/p 0.8 -30 0.10 0.8 0.9 1 1.1 1.2 0 0.8 0.9 1 1.1 1.2 frequenc 1 --5400 1 0.7 ∆p=0.1,fm=01.1f0 Modulation 11..24 --7600 1 FDTD FDTD 1.6 .5 0.6 Floquet Floquet -80 0.8 0.5 1.8 -90 E(t)ttF{}00..34 E(t)trF{}00..46 20 0.5 1ffrreeqquueenn1cc.yy5,,ωω//ωω00 2 2.5 3 -100 20log|Et(ω0)| 0.2 0.2 Fig. 3. Transmission spectrum of a time-only modulated metasurface 0.1 for a varying modulation frequency and fixed modulation depth, ∆p. The 0 0 metasurfaceisassumedtobematchedwithχee=χmm. 0.5 1 1.5 0.5 1 1.5 ∆p=0.25,fm=0.25f0 Next,anexampleofaspace-timemodulatedmetasurfaceis 0.7 0.7 considered. Due to the spatial dependence of the susceptibili- FDTD FDTD 0.6 0.6 Floquet Floquet ties following (6), the metasurface is spatially discretized and 0.5 0.5 the Floquet analysis is performed at each x location to de- E(t)ttF{} 00..34 E(t)trF{}00..34 taenrdmEi˜nre(xth,eztr=an0s+m,itωte)d, raensdpercetfliveecltye.dAfineledxsa,mE˜pt(lex,iszs=ho0w+n,ωin) Fig. 4, for an input Gaussian beam, showing the spectrum 0.2 0.2 of the transmitted fields at the output of the metasurface at 0.1 0.1 each location, with the mismatched parameters of Fig. 2. The 0 0 0 0.5 1 1.5 2 2.5 0 0.5 1 1.5 2 2.5 Gaussianwaveformnatureisclearlyobservedforeachspectral frequency,ω/ω0 frequency,ω/ω0 component. More interestingly, the phase profile across the Fig.2. Spectrumoftransmissionandreflectionfieldsfromaninfinitetime- metasurface as a function of frequency is also shown in only (βp = 0) modulated metasurface incident with a normally incident Fig. 4. A linear phase gradient with different spatial slopes, plane-wave, for the case of a) ∆p = 0.05, ωp = 0.05ω0, b) ∆p = 0.1, d∠E (x)/dx, across the metasurface is clearly evident for ωp =0.1ω0 andc)∆p =0.25,ωp =0.25ω0.Totalnumberofharmonics t 2N + 1 = 61 and T0 = 100 fs in (16). The metasurface parameters each harmonic component, except the fundamental frequency. are: ωe0 = 2π(224.63 THz), ωep = 0.36 Trad/s, αe = 500×109, Therefore, each of these harmonics are expected to refract ωm0 = 2π(224.40 THz), ωmp = 0.29 Trad/s, αm = 100×109, and at different angles. Similar observations are also made for excitationfrequencyω0=2π(230THz).Thespectrumisnormalizedtothe non-modulatedcasewith∆p=0. reflected fields (not shown here). To confirm the refraction angles of the transmitted and reflected fields from the metasurface, the output fields at a While Fig. 2 showed results for few discrete points for specific harmonic frequency ω , are then Fourier propagated ∆ and ∆ , Fig. 3 shows a spectral map for a continuously n p m in free-space as [13] variedmodulationfrequencyω ,foragivenmodulationdepth p 7only for simplicity, and clearer compact description of the associated matrices. E˜t(x,z,ωn)=Fx−1[Fx{Et(x,z+,ωn)}exp{−jkz,nz}] MANUSCRIPTDRAFT 6 20log|E(ω0−ωp)| 0 20log|E(ω0)| 0 20log|E(ω0+ωp)| 0 30 30 30 15 15 15 -20 -20 -20 0 0 0 m m m µ µ µ x -15 x -15 x -15 -40 -40 -40 -30 -30 -30 -200 0 200 -200 0 200 -200 0 200 zµm zµm zµm 1 1 1 θ−1=−8.67◦ θ0=0◦ θ+1=+7.08◦ 2,ω)0}±00..68 kxp,e−ak1 2,ω)0}±00..68 20,ω)0}± 00..68 kxp,e+ak1 0 0 = E(x,z=xF{00..24 kx,in E(x,z=xF{00..24 E(x,zxF{ 00..24 kx,in kx,in 0 0 0 -2 -1 0 1 2 -2 -1 0 1 2 -2 -1 0 1 2 spatialfrequencykx×10−6 spatialfrequencykx×10−6 spatialfrequencykx×10−6 Fig.5. Steady-statescatteredfieldsinthetransmissionandreflectionregioncorrespondingtothespace-timemodulatedmetasurfaceofFig.4,obtainedusing analyticalFourierpropagationmethod,forthefirstup-convertedanddown-convertedspectralcomponentalongwiththefundamental,inagreementwiththe predictionsof(4).AlsoshownarethespatialFouriertransformsoftheindividualharmonicswithsolidanddashedcurvescorrespondingtotransmittedand reflectedfields,respectively.Allfieldsarenormalizedtotheirownmaximum. ∆p=0.025,fm=0.1f0 12.5 20log|Et(x,ω)| 0 12.5 6 Et(x,ω)rad 3 θ(ωn =ω0+nωp)=sin−1(cid:34)kxpe,nak(cid:35). (18) k 7.5 -10 7.5 2 n 1 m 2.5 -20 mm 2.5 where kxpe,nak is the peak location of the spatial spectrum of the xµ -2.5 -30 xµxµ -2.5 0 output fields, kn is the wavenumber of the nth harmonic, and -1 the angle θ is measured from the normal of the metasurface. -7.5 -40 -7.5 -2 For the n = 1 harmonics, the refraction angles are found ± -12.5 -50 -12.5 -3 to be θ = 8.67◦ and θ = 7.08◦, respectively, for both 0.5 0.75 1 1.25 1.5 0 0.5 1 1.5 2 transmission and reflection, with a very good agreement with frequency,ω/ω0 frequency,ω/ω0 the theoretical values of 8.33 and 6.81 obtained from (4). ◦ ◦ Fig. 4. Transmission spectrum (amplitude and phase) of a space-time modulated metasurface, when incident with a normally incident Gaussian beam.ThebeamisgivenbyE0(x)=exp{−(x/2wx)2},withwx=5µm. IV. CONCLUSIONS Themetasurfacesize(cid:96)=25µmandtheexcitationfrequencyf0=230THz. Thepumpingspatialfrequencyβp=5π/(cid:96). A rigorous Floquet mode analysis has been proposed for a zero thickness space-time modulated Huygens’ metasur- face to model and determine the strengths of the harmonics E˜ (x,z,ω )= 1[ E (x,z ,ω ) exp +jk z ] of the scattered fields. The proposed method is based on r n Fx− Fx{ r + n } { z,n } GSTCs treating the metasurface as a spatial discontinuity. (cid:112) where k = k2 k2, is the free-space wave number in The metasurface has been modelled using space-time varying z,n n− x the z direction. Fig. 5 shows the resulting field in the x z resonant frequencies of the associated electric and magnetic − plane for the first down-converted (ω ω ) and first up- surface susceptibilities, χ and χ , respectively. The un- 0 p ee mm − converted harmonic (ω +ω ), in addition to the fundamental knownscatteredfieldshavebeenexpressedintermsofFloquet 0 p frequency ω . As expected from Fig. 4 and consistent with modes, which when used with the GSTCs, led to a system of 0 (4), the fundamental frequency ω propagates without any field matrix equations. The resulting set of linear equations 0 refraction. On the other hand, the down-converted and up- were then solved to determine the transmitted and reflected converted harmonics are refracted in the lower and upper half scattered fields. Using an FDTD solver, the proposed method of the x z plane, respectively. Their reflection angles using isvalidatedandconfirmedforseveralexamplesofmodulation − transverse wavenumbers are given by depths (∆ ) and frequencies (ω ). Finally, the computed p p MANUSCRIPTDRAFT 7 steady-state scattered fields are Fourier propagated analyti- cally,forconvenientvisualizationofrefractedharmonics.The [C] proposed method is fast, simple and versatile, and is expected (cid:122) (cid:125)(cid:124) (cid:123) 0 0 0 to be a useful tool in determining the steady-state scattered  0 0 0  fieldsfromaspace-timemodulatedHuygen’smetasurface,for  d[V] arbitrary modulation frequencies and depths.  W1 0 0  dt  0 T1 0  0 0 T 1 V. APPENDIX [G(t)] (cid:122) (cid:125)(cid:124) (cid:123)   A A A Summary of Finite-Difference Time Domain (FDTD) Formu- 1 2 3 lation of Space-Time Modulated Metasurface [16][17]  B1 B2 B3    + W2(t) 0 0 [V]=[E(t)], (21) Consider a space-time modulated metasurface of Fig. 1  0 T2(t) 0  with the plane-wave input and output fields described in (7). 0 0 T (t) 2 The field equations governing the transmitted and reflected where [V] is the solution vector containing the transmission fieldsarethetime-domainLorentzrelations(8)andtheGSTC and reflection field E and E , and [G(t)] contains the exact equations (9). t r description of the metasurface. For a space-time modulated Each Lorentz resonator equation of (8) corresponding to metasurface considered here, the resonant frequencies of the incident, transmitted and reflected fields, is a second order Lorentziansusceptibilitiesareassumedtobeafunctionofboth differential equation. By introducing an auxiliary variable, spaceandtime,i.e.ω (x,t)andω (x,t).Theabovematrix each second-order differential equation can be decomposed e0 m0 equationscanbere-writteninacomplexmatrixequationform into two first-order differential equations. For instance, the as electric and magnetic polarizabilities corresponding to the incident fields are given by d[V] [C] +[G(t)][V]=[E(t)], (22) dt dQ ω Q¯ = 0 +α Q , (19a) which now can be easily solved using standard finite- e0 0 e 0 dt difference technique based on trapezoidal integration as dQ¯ 1 dω ω2 dt0 +Q¯0ωe0 dte0 +ωe0Q0 = ωeep0E0 (19b) [V] =(cid:18)[C]+ ∆t[G] (cid:19)−1(cid:20)∆t[E]i+[E]i−1 i 2 i 2 (cid:18) (cid:19) (cid:21) dM ∆t ω M¯ = 0 +α M , (19c) + [C] [G] [V] . m0 0 dt m 0 − 2 i i 1 − dM¯ 1 dω ω2 0 +M¯0 m0 +ωm0M0 = mp E0, (19d) REFERENCES dt ω dt −η ω m0 0 m0 [1] A.Arbabi,Y.Horie,M.Bagheri,andA.Faraon,“Completecontrolof where Q¯ and M¯ are the two unknown auxiliary variables polarizationandphaseoflightwithhighefficiencyandsub-wavelength 0 0 spatialresolution,”arXiv:1411.1494,pp.4308–4315,Nov2014. in addition to Q0 and M0, for a specified input field E0. [2] N.Yu,P.Genevet,M.A.Kats,F.Aieta,J.-P.Tetienne,F.Capasso,and Similar set of equations, can be developed for the transmitted Z.Gaburro,“Lightpropagationwithphasediscontinuities:Generalized lawsofreflectionandrefraction,”Science,vol.334,no.6054,pp.333– and reflected fields in terms of unknown fields E and E . t r 337,2011. They all represent a set of 12 equations, with 12 auxiliary [3] N.YuandF.Capasso,“Flatopticswithdesignermetasurfaces,”Nature unknowns and 2 primary unknowns. All the scattered fields Materials,vol.13,April2014. [4] M. Kerker, The Scattering of Light and Other Electromagnetic Radia- and their respective polarizations are coupled to each other tion. AcademicPress,NewYork,1969. through the GSTC equations, which can be written in terms [5] M.Decker,I.Staude,M.Falkner,J.Dominguez,D.N.Neshev,I.Brener, of new auxiliary variables as T. Pertsch, and Y. S. Kivshar, “High-efficiency dielectric huygens’ surfaces,”Adv.Opt.Mater.,vol.3,no.6,pp.813–820,2015. [6] S.JahaniandZ.Jacob,“All-dielectricmetamaterials,”NatureNanotech- ω (M¯ +M¯ +M¯ ) α (M +M +M ) nology,vol.2,no.11,pp.23–36,Jan2016. m0 0 t r m 0 t r [7] C. Pfeiffer and A. Grbic, “Metamaterial huygens’ surfaces: Tailoring − 2 wave fronts with reflectionless sheets,” Phys. Rev. Lett., vol. 110, p. = (Et Er E0), (20a) 197401,May2013. µ − − 0 [8] G.V.EleftheriadesandS.Michael,ControllingElectromagneticWave- frontsUsingHuygens’Metasurfaces. JohnWiley&Sons,Inc.,1999. ω (Q¯ +Q¯ +Q¯ ) α (Q +Q +Q ) [9] Y.Hadad,D.L.Sounas,andA.Alu,“Space-timegradientmetasurfaces,” e0 0 t r − e 0 t r Phys.Rev.B,vol.92,p.100304,Sep2015. 2 [10] A.Shaltout,A.Kildishev,andV.Shalaev,“Time-varyingmetasurfaces = (E E E ). (20b) (cid:15) η 0− t− r and lorentz non-reciprocity,” Opt. Mater. Express, vol. 5, no. 11, pp. 0 0 2459–2467,Nov2015. [11] R.S.ChuandT.Tamir,“Wavepropagationanddispersioninspace-time Equations (19) and (20), finally represents 14 equations with periodicmedia,”ElectricalEngineers,ProceedingsoftheInstitutionof, 14 unknowns, and can be expressed in a matrix form as vol.119,no.7,pp.797–806,July1972. MANUSCRIPTDRAFT 8 [12] E. S. Cassedy and A. A. Oliner, “Dispersion relations in time-space periodic media: Part 1 stable interactions,” Proceedings of the IEEE, vol.51,no.10,pp.1342–1359,Oct1963. [13] J. Goodman, Introduction to Fourier Optics. Roberts and Company Publishers;3rdEditionedition,2004. [14] B. E. A. Saleh and M. C. Teich, Fundamentals of Photonics, 2nd ed. Wiley-Interscience,2007. [15] R.-S. Chu and T. Tamir, “Guided-wave theory of light diffraction by acousticmicrowaves,”IEEETrans.Microw.TheoryTech.,vol.18,no.8, pp.486–504,Aug1970. [16] T. Smy and S. Gupta, “Exact finite-difference time-domain mod- elling of broadband huygens’ metasurfaces with lorentz dispersions,” ArXiv:1609.05575,Sept.2016. [17] S. A. Stewart, T. Smy, and S. Gupta, “Finite-difference time- domain (fdtd) modelling of space-time modulated metasurfaces,” ArXiv:1612.02087,Dec.2016. [18] E. F. Kuester, M. A. Mohamed, M. Piket-May, and C. L. Holloway, “Averagedtransitionconditionsforelectromagneticfieldsatametafilm,” IEEE Transactions on Antennas and Propagation, vol. 51, no. 10, pp. 2641–2651,Oct2003. [19] Y.Vahabzadeh,K.Achouri,andC.Caloz,“Simulationofmetasurfaces in finite difference techniques,” IEEE Transactions on Antennas and Propagation,vol.64,no.11,pp.4753–4759,Nov2016. [20] S. Sandeep, J. M. Jin, and C. Caloz, “Finite element model- ing of metasurfaces with generalized sheet transition conditions,” ArXiv:1611.06253,Dec.2016. [21] M. M. Idemen, Discontinuities in the Electromagnetic Field. John Wiley&Sons,2011. [22] K. Achouri, M. Salem, and C. Caloz, “General metasurface synthesis basedonsusceptibilitytensors,”IEEETrans.AntennasPropag.,vol.63, no.7,pp.2977–2991,July2015. [23] S.GuptaandC.Caloz,“Perfectdispersivemediumforreal-timesignal processing,”IEEETransactionsonAntennasandPropagation,vol.64, no.12,pp.5299–5308,Dec2016. [24] K.Achouri,A.Yahyaoui,S.Gupta,H.Rmili,andC.Caloz,“Dielectric resonator metasurface for dispersion engineering,” IEEE Transactions onAntennasandPropagation,vol.PP,no.99,pp.1–1,2016.