Electrically Tunable Harmonics in Time-modulated Metasurfaces for Wavefront Engineering

Modulation of metasurfaces in time gives rise to several exotic space-time scattering phenomena by violating the reciprocity and generation of higher-order frequency harmonics. We introduce a new design paradigm for time-modulated metasurfaces, offering electrically tunable engineering of the generated frequency harmonics and their emerging wavefronts by controlling the phase delay in modulation. It is demonstrated that the light acquires a dispersionless phase shift regardless of incident angle and polarization, upon undergoing frequency conversion in a time-modulated metasurface which is linearly proportional to the modulation phase delay and the order of generated frequency harmonic. The conversion efficiency to the frequency harmonics is independent of modulation phase delay and only depends on the modulation depth and resonant characteristics of the metasurface element, with the highest efficiency occurring in the vicinity of resonance, and decreasing away from the resonant regime. The design approach allows for creating tunable spatially varying phase discontinuties with 2{\pi} span in the wavefronts of generated frequency harmonics for a wide range of frequencies and incident angles. Specifically, we apply this approach to a time-modulated metasurface in the Teraherz regime consisted of graphene-wrapped silicon microwires. For this purpose, we use an accurate and efficient semi-analytical framework based on multipole scattering. We demonstrate the utility of the design for tunable beam steering and focusing of the generated frequency harmonics. Furthermore, we rigorously verify the broadband and wide-angle performance of the metasurface in manipulation of the generated frequency harmonics. The proposed design approach enables a new class of high-efficiency tunable metasurfaces with wide angular and frequency bandwidth, wavefront engineering capabilities and multi-functionality.

Modulation of metasurfaces in time enables several exotic scattering phenomena by breaking the time-reversal symmetry and generation of frequency harmonics. We introduce a novel design paradigm for time-modulated metasurfaces, offering electrically tunable engineering of the generated frequency harmonics and their emerging wavefronts by controlling the phase delay in modulation. It is demonstrated that the light acquires a dispersionless phase shift regardless of incident angle and polarization, upon frequency transitions in a time-modulated metasurface which is determined by the modulation phase delay and the order of generated harmonic. The conversion efficiency to the frequency harmonics is independent of modulation phase delay and only depends on the modulation depth and resonant characteristics of the metasurface element, with the highest efficiency occurring in the vicinity of resonance, and decreases away from the resonance. The design approach allows for creating tunable spatially varying phase discontinuties with 2π span in the generated frequency harmonics for a wide range of frequencies and incident angles. Specifically, we apply this approach to a time-modulated metasurface in the Teraherz regime consisted of graphene-wrapped silicon microwires. For this purpose, we use an accurate and efficient semi-analytical framework based on transition matrix. We demonstrate the utility of the design for tunable beam steering and focusing of the generated harmonics. Furthermore, we rigorously verify the broadband and wideangle performance of the metasurface in manipulation of the generated harmonics. The proposed design approach enables a new class of high-efficiency tunable metasurfaces with wide angular and frequency bandwidth, wavefront engineering capabilities and multi-functionality.
Optical metasurfaces have had a significant impact on development of novel optical devices owing to their capability in the modulation of electromagnetic wavefronts of light and their thin structures which can reduce the size of optical platforms by replacing the bulky components. They have been used for numerous applications such as beam steering [1][2][3][4], focusing [5,6], holography [7,8], etc. These functionalities are realized by engineering the phase discontinuities across the metasurface according to the generalized Snells law [9,10]. In conventional metasurfaces, generally two different approaches have been adopted for this purpose. In the first approach, the structural parameters of a resonant structure are varied across the metasurface which can provide a desired phase profile for linearly polarized light due to the phase agility at the resonance [5,6]. In the second approach, a half-wave plate element is rotated around its axis which leads to a geometric phase shift equal to twice as the rotation angle for a circularly polarized incident wave based on the ParancharatnamBerry (PB) design rule [11][12][13].
However, one major limitation associated with these conventional metasurfaces is their static optical response after fabrication which ties them to a specific application they were designed for. To overcome this limitation, an extensive effort has been put into post-fabrication tuning of geometrically-fixed metasurfaces by exploiting mechanical [14,15], thermal [16,17], optical [18,19] and electrical [20,21] effects which renders more flexibility * hosseinm@ece.neu.edu in the operation, allowing for multifunctionality and ondemand manipulation of light. Among all the tunability mechanisms, electrical tunability of electro-optical materials via field-effect modulation holds the greatest promise due to the continuous tunability, wide tunability range and shorter response time comparing to mechanical reconfiguration and thermal phase transition. Furthermore, unlike all-optical tunability, electrical tunability easily allows for independent addressing and biasing of each element without the need for complex lens systems which enables wavefront engineering. In particular, graphene and indium tin oxide (ITO) have attracted a lot of attention in the account of their compatibility with silicon technology, large scale fabrication capability and low-dimensionality of the active regions. The carrier concentration in these materials can be tuned through electrostatic gating in a parallel capacitor configuration which can be translated into the change in the optical constants of the material through carrier-dependent dispersion models in the infrared (IR) and Terahertz (THz) frequencies. However, the change in the carrier concentration is limited to a thin accumulation layer for ITO and to a two-dimensional surface for graphene. Due to this fact, these materials have been integrated into resonant geometries which can enhance the light-matter interaction in the active regions and enable tunable modulation of phase over a wide range in reflection [20][21][22][23][24] or transmission [25] through spectral shift of the resonance. Despite the fruitful progress that has been made toward electrically tunable modulation of light wavefront, all the proposed designs thus far suffer from a limited efficiency as they cannot cover the 2π span for phase modulation.
Moreover, the afforded phase modulation comes at the cost of scarifying the amplitude efficiency which hinders the efficient performance of these metasurfaces. Furthermore, these metasurfaces exhibit extremely narrow bandwidths in both angular and frequency responses due to the strong resonant dispersion which makes them more susceptible to the fabrication errors and limits their functionality to a specific wavelength and incident angle. As such, realization of an electrically tunable metasurface covering the 2π span for the phase modulation with wide frequency and angular bandwidths is a challenge yet-tobe addressed. In this contribution, we offer a new route for overcoming these limitations by introducing a dispersionless non-resonant phase shift of light upon frequency transitions in time-modulated metasurfaces.
Despite the real-time tunability offered by electrooptical metasurfaces, their operation has remained mostly quasi-static as the variations in time are considered to be negligible. The electro-optical materials such as graphene and ITO offer the opportunity to implement time-varying metasurfaces by changing the external bias in time. Due to the quick response of these materials, they allow for modulation frequencies up to several gigahertz [26]. As such, they can be biased using radiofrequency (RF) circuits to realize time-modulated metasurfaces in IR and THz frequencies. Time modulation in a metasurface enables several exotic space-time scattering phenomena by breaking the time-reversal symmetry and generation of frequency harmonics [27][28][29][30][31][32][33][34][35][36][37]. In addition to the electrical tunability, time-modulation adds a new dimensionality to the design space which is expected to overcome several limitations associated with static or quasi-static metasurfaces.
In this work, we develop a novel design paradigm for time-modulated metasurfaces which is able to manipulate the generated frequency harmonics and shape their wavefronts with electrical tunability. We demonstrate that by changing the phase delay in the electrical modulation of a metasurface using phase shifters, each generated frequency harmonic will pick up a dispersionless phase shift proportional to its corresponding order in forward and backward scattered plane waves covering 2π span while maintaining a constant amplitude. This opens up the possibility for efficient wavefront engineering of frequency harmonics. The frequency conversion efficiency of the metasurface is determined by the geometry of building blocks and their resonant characteristics, exhibiting a peak in the vicinity of the resonance, and decreases moving away from the resonance. A time-modulated metasurface is implemented in the THz regime by using graphene-wrapped silicon microwires as building blocks and modulation of graphene conductivity through an RF biasing network incorporating phase shifters. We adopt an accurate and efficient semi-analytical framework based on transition matrix for modeling the metasurface which eliminates the need for computationally expensive numerical methods and allows for the fast design and characterization of time-modulated metausr-faces. The proposed metasurface can act simultaneously as tunable generator and manipulator of frequency harmonics. The application of the design in wavefront shaping functionalities is investigated and the broadband wide angle performance is rigorously verified. The rest of this manuscript is organized as follows. In Section 2, the theory of the design rule is established for a general timemodulated metasurface. In Section 3, we implement an electrically tunable time-modulated metasurface based on graphene-wrapped silicon microwires. Tunable wavefront engineering functionalities are investigated and the results are discussed in Sec. 4. The broadband and wide angle performance is demonstrated in Section 5. Finally, the conclusions are drawn in Section 6.

I. DESIGN RULE
In order to simplify the analysis and gain insight to the scattering of light from time-varying metasurfaces, we utilize the equivalent surface admittance representation in this section to establish the required design rule through analytical relations for reflection and transmission coefficients. In the subsequent sections, we verify the design rule for a time-modulated metasurface with realistic structure implemented by graphene-wrapped microwires as building blocks, analyzed without any approximations using a robust semi-analytical framework based on transition matrix formulation and multipole scattering theory [38].
Let us consider a periodic array consisted of subwavelength elements varying in time. We assume the size of elements are small enough comparing to the wavelength, so that the array can be modeled with an equivalent timevarying surface admittance. This homogenization model is valid only for deeply subwavelength elements in a periodic arrangement and it has not been used in following sections.
Extending the well-known transfer matrix formalism, the tangential electromagnetic fields in both sides of the surface admittance are expanded into backward and forward propagating plane waves with time-varying complex amplitudes a m (t) and b m (t) as shown in Fig. 1. We assume the temporal variations of the admittance to be small compared to oscillations of excitation frequency such that the dispersion effects induced by timemodulation can be ignored. The validity of this assumption is ensured by the range of accessible modulation frequencies for electro-optical materials in the THz and IR frequency regimes. As such, the transfer matrix equation relating the amplitudes of the fields on opposite sides of the surface admittance can be written in multiplicative form in time-domain for each excitation frequency (ω 0 ) as [39]: where T (t) is the time-varying transfer matrix, Y s (ω 0 , t) is the time-varying surface admittance of the metasurface and Y 0 is the admittance of plane wave in free-space which depends on the incident angle and polarization.
Here, we do not set any limit on Y 0 , so that the results are valid irregardless of incident angle and polarization. It should be remarked while ignoring the modulationinduced dispersion effects, we consider dependency of admittance temporal profile to the excitation frequency, thus accounting for the change in material parameter with optical frequency. The transfer matrix equation can be taken into angular frequency domain by taking the Fourier transform of both sides as: where * denotes the convolution. Equation (2) implies that an input frequency ω 0 will be converted to a spectrum of output frequencies. When the elements are varying with a periodic modulation having a modulation frequency of ω m , the transfer matrix is also periodic and can be expanded in form a Fourier series as: And the Fourier transform takes the following form: Therefore, the generated frequency harmonics will be a discrete spectrum consisted of the input frequency upand down-modulated by the time modulation frequency. Substituting Eq. (4) into (2) while choosing ω = ω q = ω 0 + qω m and ω = ω p = ω 0 + pω m with q, p = · · · , −1, 0, +1, · · · , we will arrive at the following system of equations: The above matrix equation relates the complex amplitudes of the forward and backward propagating frequency harmonics on opposite sides of the metasurface in the frequency domain. By setting B 2 (ω p ) = 0 for all p's and A 1 (ω q ) = 0 for all q's except A 1 (ω 0 ) = 1, we can solve the equation to obtain the reflection and transmission coefficients of all generated harmonics for a monochromatic excitation of ω 0 incident from the top, as: By assuming the time-modulation profile as a simple offset sinusoid, the surface admittance of the metasurface can be expressed as: In the above relation Y 0 s (ω 0 ) is the average surface admittance and ∆Y s is the admittance modulation depth which depends on the resonant characteristics of the building block. Substituting the admittance into the time-varying transfer matrix and comparing with (3), we will obtain the Fourier coefficients of the transfer matrix as: Plugging (9)-(11) into (5) results into a tri-block diagonal matrix equation which can be readily solved to obtain the transmission and reflection coefficients of the generated harmonics. Now, let us consider adding a phase delay of α to the sinusoidal modulation profile. The Fourier coefficients of the transfer matrix will be subsequently changed as where subscript α denotes the quantities corresponding to a time-modulation with a phase delay of α. Writing transfer matrix equation (5) for the phase-delayed modulation and using T (ω 0 ), one can easily obtain: which results into: These equations imply that by changing the modulation phase delay (α), the generated frequency harmonic of order n will acquire a phase shift of nα in forward and backward scattering planes while maintaining a constant amplitude. This provides a span of 2nπ for the phase shift of n-th harmonic by changing α from 0 to 2π which can be used as a design principle for wavefront engineering of generated frequency harmonics. Furthermore, according to the formulation, the phase shift introduced by modulation phase delay upon frequency transitions in the time-modulated metasurface is non-resonant, dispersionless and independent of incident angle and polarization. Due to the conjugate phase of up-and down-modulated harmonics and the difference between the acquired phase of different harmonic orders, the design principle leads to several exotic phenomena such as anomalous bending, spatial decomposition of harmonics and dual-polarity focusing/diffusing. In this context, the proposed design rule resembles the PB design rule in which a circularlypolarized light undergoing polarization conversion in a half-wave plate picks up a dispersionless geometric phase shift simply by rotating elements where the phase shift is conjugate for circular polarizations of opposite handness [11,12].
The established design rule allows for realization of advanced functionalities through electrically phased timemodulated metasurfaces implemented by using electrooptical materials such as graphene and ITO. The optical constants of these materials can be modulated up to several GHz which can be achieved by using a radiofrequency (RF) biasing network. The phase delay in modulation can be realized and electrically tuned using RF phase shifters. It should be remarked that the phase shifting of electrical modulation is sufficient instead of a true time-delay. Multiple techniques can be adopted to realize passive and active RF phase shifter circuits [40] which have been extensively used in phased array microwave antennas. It should be noted that conventional phased-array antennas are time-invariant and utilize RF phase shifters in their feeding network to generate a phased radiation. On the contrary, a phased time-modulated metasurface relies on an optical pump (plane wave incidence) for radiation while the scatterers are modulated in time with an RF electrical network.
The developments made in this section are agnostic to the shape and materials of the building blocks; implying that the design rule can be implemented with arbitrary-shaped unitcells incorporating different types of tunable materials in different frequency regimes. In the following, we will synthesize a time-modulated metasurface in the THz regime using graphene-wrapped silicon microwires and apply this design rule to explore beamshaping and exotic phenomena in more depth. All the results in the following sections are obtained with an exact semi-analytical method based on multipole scattering [38] without using equivalent admittance approximation.

II. IMPLEMENTATION OF TIME-MODULATE METASURFACE
Graphene is an excellent candidate as a tunable electro-optical material in the mid-IR and THz regimes in the account of its low-loss, low-dimensionality, large-scale fabrication feasibility, broadband tunability and short response time. The conductivity of graphene can be modulated over a wide range using electrical gating in parallel capacitor configuration, at speeds up to several GHz [26,[41][42][43][44][45]. It has been used in many active nanophotonics devices for a variety of applications such as wide-band absorption [46] and electrically tunable phase modulation [21]. For these applications, graphene is integrated into resonant geometries to enhance light-graphene interaction and achieve a large tunability. In particular, in ref [21], graphene has been incorporated into a resonant plasmonic metasurface with metal-insulator-metal configuration to achieve a large phase modulation in the reflection. However, the phase shift is limited to 230 degrees and exhibits a narrow bandwidth. We overcome these limitations through applying the proposed design rule to a time-modulated metasurface and exploiting the non-resonant phase shift of light upon frequency transitions.
We choose to operate in the THz regime where graphene yields the largest tunability. In order to implement the metasurface, we have used graphenewrapped micorwires as the building blocks which can enhance light-graphene interaction through resonant characteristics associated with the Mie scattering. Such wire elements have been fabricated via chemical-vapordeposition (CVD) growth of monolayer graphene on wires [47,48] or draping graphene flakes over wires with an adhesive tape [19]. Furthermore, they have been envisioned for a variety of applications such as tunable modulation of absorption [49], cloaking [50,51], wave guiding [52][53][54], nonlinear harmonic generation [55], etc. Moreover, the well-defined geometry of microwires allow for developing semi-analytical techniques based on multipole scattering which enable fast design and accurate characterization of space-time gradient metasurfaces. It should be emphasized that simulation of time-modulated metasurfaces poses a great challenge for a time-domain analysis like finite-difference time-domain (FDTD) method as the required simulation times to capture the steadystate solution are much longer comparing to the timeinvarient counterparts which needs tremendous computational resources and can quickly become prohibitive for large structures with subwavelength geometrical features. As such, a semi-analytical framework is of great benefit for unraveling the novel physics of time-modulated metasurfaces which have remained unexplored.
Here, we use a recently developed technique based on transition matrix (T-matrix) developed by authors which allows for robust characterization of nearfield and farfield performance of periodic and aperiodic arrays of multimaterial microwires coated by time-modulated graphene layers above layered substrates [38]. To this purpose, we obtain a generalized T-matrix for time-varying cylindrical interfaces which relates the scattered field coefficients corresponding to a specific cylindrical mode and frequency harmonic to the incident field counterparts. The aggregate T-matrix approach is used to treat multimaterial configurations. The coupling effects in periodic and aperiodic arrays of nanowires are dealt with, by translating cylindrical harmonics through Graf's addition theorem and the substrate effect is modeled by using scattering matrix approach in periodic arrangements and evaluation of Weyl-type integrals in aperiodic arrangements. All the presented results in the paper are obtained with this exact semi-analytical framework. The validity of the method is rigorously verified by comparing the results with full-wave FDTD simulations which is brought in the supplemental material.
The unitcell of the implemented metasurface is demonstrated schematically in Fig. 2. It is composed of an array of subwavelength microwires with a periodicity of Λ = 130µm on a silica substrate with a thickness of 80µm. The topology of microwires are chosen to meet two requirements: (i) support Mie resonances to increase modulation depth and frequency conversion efficiency, and (ii) allow for biasing graphene in a parallel capacitor configuration. The cores are made of silicon with radius of R core = 40µm which are subsequently coated by thin insulating SiO 2 layers with radial thickness of 20 nm, in accordance to the biasing requirements. The electrical biasing can be done by applying a voltage between the silicon core and the graphene which requires doping of the silicon [56]. Here, we use an n-type silicon with the background concentration of n = 10 15 cm −3 . The silicon cores are connected to a ground plane while a time-varying electrical voltage is applied to the graphene layers through a 2D biasing grid which allows addressing and biasing each element, independently. It should be mentioned such electrical gating configurations have been used for realization of field-effect transistors using graphene-wrapped silicon microwires [57].
The complex permittivity of silicon is obtained via a Drude model [58] as (ω) = inf − is the high-frequency permittivity, and ω p and Γ are plasma frequency and damping constant, respectively. The plasma frequency is related to carrier concentration (n) through ω 2 p = ne 2 0m * in which e is the electron charge, 0 is the vacuum permittivity and m * is the effective mass of electron. For silicon, we use inf = 11.7, Γ = 180THz and m * = 0.27m e .
The frequency-dependent conductivity of monolayer graphene is dominated by intraband transitions in the THz regime which can be expressed as [59]: where e is the electron charge, is the reduced Planck's constant, k B is the Boltzmann constant, T is temperature, τ and E f are the scattering time and Fermi energy level of graphene, respectively. At the room temperature k B T = 25.7 meV and E f >> k B T for a moderately doped graphene. As such, the conductivity can be approximated as: The relationship between the Fermi level, E f , of the graphene sheet and the gate voltage, V g , in a parallel capacitor configuration can be expressed as [59]: where C ox is the geometric capacitance of the gate oxide and V Dirac is the gate voltage at which minimum conductance is observed. According to (17) and (18), the conductivity of graphene follows a temporal profile proportional to the square root of the modulation voltage.
Applying an external bias with a sine squared waveform such that E f = E f 0 (1 + ∆E f sin(ω m t + α)), the temporal profile of graphene conductivity for an excitation frequency of ω 0 can be written as: Throughout this work, the scattering time is considered as τ = 0.5 ps [60]. It has been shown experimentally that the Fermi energy level of graphene can be modulated in the range of 0-0.6 eV using electrical gating [21]. In order to achieve the largest modulation depth and highest efficiency for frequency conversion in the time-modulated metasurface, we choose E f 0 = 0.3 eV and ∆E f = 0.3 eV. The transmission amplitude spectrum of the metasurface in the un-modulated case is shown in Fig. 2(b) with a Fermi energy level of E f 0 = 0.3 eV for a normally incident plane wave with transverse electric (TE) polarization in which the magnetic field is along microwire axis. The transmission dips correspond to the resonant modes which are characterized according to their magnetic fields mode profiles shown in the figure. The sinusoidal modulation of graphene in time will lead to generation of frequency harmonics ω 0 ± nω m for a monochromatic excitation frequency of ω 0 . The modulation frequency is chosen as f m = 1 GHz which is well within the practical range [41,42]. Moreover, the generated frequency harmonics can be resolved with the resolution of THz detectors [61]. The transmission and reflection amplitude spectra of the generated harmonics by the timemodulated metasurface are shown in Figs. 2(c) and (d), respectively. As it can be observed from the results, the higher-order frequency harmonics are generated in both forward and backward directions due to symmetry and absence of incident component at these frequencies. The highest frequency conversion efficiencies are obtained in the vicinity of resonant wavelengths of fundamental harmonic due to larger scattering modulation depth and stronger light-graphene interaction. Higher conversion efficiency can be maintained over the operating frequency range by exploiting higher order modes and increasing diameter of microwires. Here, the diameter is chosen to ensure the coupling between neighboring elements is weak which leads to a better performance in local wavefront control. An extended discussion on the conversion efficiency of the time-modulated metasurface is included in section 1 of the supplemental material.
We do remark that despite the high conversion efficiency afforded by resonant characteristics of microwires, the amplitudes of generated frequency harmonics are still small compared to the fundamental frequency. Increasing the modulation depth (∆E f ) will enhance the frequency conversion efficiency of the metasurface. Here, we have limited ourselves to the range E f = 0 − 0.6 eV which has been shown to be accessible by electrical gating in experiments. It is noteworthy that Fermi energy levels as high as E f = 1−2 eV have been reported [62] which can potentially increase the modulation depth. Alternative geometries with more exotic shapes may be adopted to increase the conversion efficiency which has not been the focus of this contribution. Moreover, cascaded time-modulated layers hold a great promise for up-and down-conversion of frequencies with near-unity efficiencies through engineered asymmetric dispersions [35,36].
In order to apply the proposed design rule, we introduce a phase delay of α in the time-modulation which can be implemented and electrically tuned by incorporating passive and active RF phase shifters in the biasing network of the elements. Figure 3 demonstrates the transmission amplitude and phase shift spectra of generated frequency harmonics versus the modulation phase delay (α), for a normally incident TE-polarized plane wave. The phase shifts of frequency harmonics are shown as the pseudo-colors in the 3D plots of Fig. 3. The phase reference at each wavelength is chosen as the phase corresponding to α = 0 to show the linear variation of transmission phase with respect to modulation phase delay more clearly. The results verify that the by introducing a modulation phase delay of α, the transmission amplitude remains constant for all frequency harmonics while n-th harmonic acquires a phase shift of nα. Moreover, the phase variation with respect to modulation phase delay has the same linear profile for all wavelengths implying that the phase gradient across the metasurface is dispersionless and depends solely on the modulation phase delay which allows for broadband functionality of the metasurface. In the account of small modulation frequency, the amplitudes of up-and down-modulated harmonics are almost equal. Increasing the modulation frequency will induce stronger asymmetries in the amplitudes of up-and down-modulated frequency harmonics which is studied in section 2 of supplemental material. Similar results can be obtained for the reflection of generated harmonics which are not brought here for the sake of brevity. Furthermore, the phase response of metasurface to transverse magnetic (TM) polarization follows the same principle which is demonstrated in section 3 of supplemental material.
The amplitude and phase variations with respect to the modulation phase delay for oblique incidence are demonstrated by plotting transmission amplitude and phase shift of generated frequency harmonics as functions of modulation phase delay (α) and incident angle (θ i ) for a TE-polarized plane wave with an excitation wavelength of 300µm in Fig. 4. The phase shifts are indicated by FIG. 4. The amplitude and phase shift of transmission coefficients corresponding to the frequency harmonics generated by the metasurface, as functions of incident angle and modulation phase delay for an oblique incidence of TE-polarized plane wave with an excitation wavelength of λ = 300µm. The phase shifts are indicated by pseudo-colors with the phase reference at each incident angle chosen as the phase corresponding to α = 0 to show the phase variations with respect to modulation phase delay more clearly.
pseudo-color in 3D plots where the phase at each incident angle is measured with respect to the phase correspond-ing to α = 0 to show the linear phase variation profile with respect to modulation phase delay more clearly. As it can be observed, the frequency conversion efficiency to the harmonics changes with incident angle due to the resonant characteristics of graphene-wrapped microwires while the scattering amplitudes of generated harmonics is independent of modulation phase delay. Furthermore, the results clearly verify that by introducing a modulation phase delay of α, the of n-th harmonic acquires a phase shift of nα which implies that changing the modulation phase delay of each element can be used to realize a phase distribution across the time-modulated metasurface for all incident angles enabling a wide-angle performance.
The design principle is also verified by full-wave FDTD simulations and the comparisons between FDTD and semi-analytical simulations are included in the supplemental material section 4.

III. WAVEFRONT ENGINEERING OF GENERATED HARMONICS
In this section, we demonstrate the utility of the proposed design rule and metasurface for wavefront engineering of the generated harmonics. We consider a large finite array of 101 microwires which is geometrically fixed but the modulation phase delay can be changed by addressing and biasing each microwire, independently. As such, the metasurface is modulated in both space and time (space-time gradient). Figure 5 demonstrates the schematic of such a phased time-modulated array. The required phase shifts in modulation can be provided by an array of passive or active phase shifters in RF modulation circuit [40]. In the following, we consider a normally incident plane wave with TE polarization and wavelength of λ 0 = 300µm and investigate beam steering and focusing of generated frequency harmonics as classical examples of wavefront engineering. Due to the conjugate phase shift of the up-and down-modulated harmonics and the difference between the acquired phase shifts by different harmonic orders, several exotic behaviors can arise which are demonstrated and discussed.

A. Steering of Generated Frequency Harmonics
As the first application, we aim at steering the generated frequency harmonics toward specific directions. To this purpose, the metasurface should provide a linear phase gradient along one direction at the desired frequency harmonic, as [5]: where ∆φ n is the phase gradient provided by the metasurface for n-th frequency harmonic, k n = (ω 0 + nω m )/c is the wavenumber corresponding to the n-th frequency harmonic, θ n is the steering angle of n-th frequency harmonic and Λ is the periodicity. This phase gradient can be achieved by varying the modulation phase delay (α) linearly across the metasurface. According to the theory and the results presented in previous sections, a linear change in modulation phase delay with a gradient of ∆α results into a linear change in the phase of n-th generated frequency harmonic with a gradient of ∆φ n = n∆α. This will lead to spatial decomposition of generated frequency harmonics as they will be bent toward well-separated directions defined by (θ n = − arcsin(n∆α/k n )). In particular, up-and down-modulated frequency harmonics will have opposite phase gradients and will anomalously bend into opposite directions. We set ∆α = −k 1 Λ sin(15 • ) which is chosen to bend the first up-modulated frequency harmonic toward the angle of θ 1 = 15 • . The modulation phase delay profile across the metasurface and the wavefronts of generated frequency harmonics are shown in Figs. 6(a) and (b), respectively. As it can be observed, the metasurface is quasi-transparent for the fundamental frequency harmonic (n = 0) as it does not experience any phase shift. The first up-modulated harmonic is bent toward the desired angle of θ 1 = 15 • in both forward and backward scattering planes. The down-modulated harmonic will bend into the opposite direction by virtue of experiencing opposite phase gradient. Moreover, the second generated harmonics will bend toward θ ±2 = ± arcsin(2k 1 /k ±2 sin(15 • )) ≈ ±31 • . The bending angles of generated frequency harmonics can be tailored at will by electrically tuning the modulation phase delay across the geometrically-fixed metasurface.
We emphasize that the generation of higher-order frequency harmonics at opposite sides of the metasurface is a result of symmetry and absence of incident field component at these frequencies. A similar behavior can be observed in single-layer nonlinear metasurfaces [63] and linear metasurfaces operating at the cross-polarization [3,4,12,64] which are bound to re-radiate and scatter the light symmetrically on both sides. It should be mentioned that placing a back-mirror can be used to create a reflect-array by blocking the transmission. In this case, the design rule persists while the conversion efficiency can also be increased due to the multiple scattering between the time-modulated elements and the back mirror.
The results for such a time-modulated reflect-array design are demonstrated in section S5 of the supplemental materials.

B. Focusing of Generated Frequency Harmonics
An opportunity offered by electrically tunable metasurfaces is achieving multiple functionalities through a single platform. Here, we demonstrate focusing of generated frequency harmonics using the proposed metasurface as another application. To this purpose, the phase gradient across the metasurface should follow a quadratic profile, as [5]: where F n is the focal distance of n-th frequency harmonic. The quadratic phase profile can be realized by changing the modulation phase delay across the metasurface, accordingly. However, in this case the focusing requirements for different frequency harmonics will be different as a linear scaling of the quadratic phase profile does not correspond to a focusing phase profile.
Here, we choose the modulation phase delays in accordance to the focusing requirement of the first upmodulated harmonic at a focal distance of F 1 = 5λ = 1.5mm. The required modulation phase delay profile across the metasurface and the field intensities of the fundamental and first generated frequency harmonics are demonstrated in Figs. 7(a) and (b), respectively. Similar to the previous case, the fundamental frequency does not experience any phase gradient and is partially reflected and transmitted. The intensity of the first up-modulated harmonic shows nearly perfect focusing at the desired focal distance from the metasurface plane in both forward and backward scattering planes. Moreover, due to opposite phase gradient experience by the down-modulated harmonic, it will be strongly diffused and diverged away from the center metasurface creating a shadow region at the center. As such, this time-modulated lens has dualpolarity functionality as it is convex for up-modulated harmonic and concave for down-modulated harmonic.

IV. BROADBAND AND WIDE-ANGLE PERFORMANCE
As established by formulation and verified from the results shown in Figs. 3 and 4, the phase shift of generated harmonics solely depends on modulation phase delay and is independent of wavelength and incident angle. This offers a significant advantage for tunable beam-shaping of generated harmonics as it enables broadband and wideangle performance i.e. the metasurface can function over a wide range of frequencies and incident angles. It should be mentioned that despite the preserved functionality, the efficiency will be different at different wavelengths and incident angles due to the changes in the frequency conversion efficiency. In this section, we rigorously verify the broadband and wide-angle performance of the metasurface for focusing and anomalous beam steering.
First, we consider the time-modulated metasurface with quadratic modulation phase delay in accordance to the focusing requirement of the first up-modulated frequency harmonic at a focal distance of F 1 = 1.5 mm illuminated by a normally incident TE-polarized plane wave and study the focusing performance at six different wavelengths in the range of 160 µm-360 µm . Fixing the focal distance and size of the array, the numerical aperture (NA) of the lens remains constant for all wavelengths. The results of normalized field intensity of the first upmodulated frequency harmonic generated by the metasurface for different excitation wavelengths are shown in Fig. 8. An almost perfect focusing is achieved at the desired location for all the wavelengths which verifies the broadband functionality of the metasurface. The width of focal spot is determined by the diffraction limit and becomes larger by increasing the excitation wavelength.
The highest focusing efficiency is observed at the wavelength of λ = 300µm. It should be noted that despite the larger frequency conversion efficiency at the wavelength of λ = 210µm, the focusing efficiency is not as high which can be attributed to the larger phase variations between neighbor elements. The lowest focusing efficiency is observed at λ = 250µm which is due to the poor frequency conversion efficiency of the metasurface at this wavelength.
FIG. 8. The normalized field-intensity of the first upmodulated frequency harmonic generated by the timemodulated metasurface for different excitation wavelengths. The modulation phase delay is changing quadratically across the metasurface to focus the first up-modulated harmonic at a focal distance of F1 = 1.5 mm.
Next, we investigate the functionality of the metasurface for oblique incidence. It should be noted that when the metasurface is incident by an oblique plane wave, a phase gradient due to the tangential momentum of the incidence is added to the phase gradient achieved by modulation phase delay. As such, the acquired phase shift by the n-th frequency harmonic generated by the timemodulated metasurface upon illumination by an oblique plane wave with incident angle of θ 0 and frequency of ω 0 is equal to nα + k 0 sin(θ 0 )x with k 0 = ω 0 /c and x being the transverse direction of the metasurface.
We consider a linear phase delay modulation with a gradient ∆α = −k 1 Λ sin(30 • ). Figures 9(a) and (b) rep-resent the results for the wavefronts of fundamental and first generated frequency harmonics for a TE-polarized plane wave with incident angles of θ 0 = 15 • and θ 0 = 30 • , respectively. Due to the addition of a phase gradient by oblique incidence, the n-th frequency harmonic will bend toward the angle of θ n = − arcsin((n∆α + k 0 Λ sin(θ 0 ))/k n ). As such, the bending angles of the upand down-modulated harmonics will not be opposite of each other anymore. For θ 0 = 15 • , the bending angles are obtained as θ +1 ≈ 14 • and θ −1 ≈ −49 • . Similarly, for θ 0 = 30 • , we will have the corresponding bending angles as θ +1 ≈ 0 and θ −1 ≈ −90 • . The results clearly verify the bending angles predicted by the theory. The metasurface remains quasi-transparent for the fundamental harmonic while up-and down-modulated are bent toward different angles. It should be mentioned that the phase gradient introduced by the oblique incidence can be canceled out by tuning the phase modulation delay across the metasurface so that the up-and down-modulated harmonics bend toward opposite directions similar to the case of normal incidence. The real-time tunability along with broadband and wide-angle functionality of the metasurface allows for compensation of chromatic aberrations and aberrations resulting from oblique incidence by adaptively tuning the phase modulation delay across the metasurface to achieve aberration-free lenses for a wide range of frequencies and incident angles.
It should be remarked that the proposed design paradigm in this paper is applicable to any timemodulated metasurface with arbitrary-shaped unitcells incorporating different types of electro-optical materials such as graphene, transparent conducting oxides, semiconductors and transition metal nitrides. This allows for achieving electrical manipulation of harmonics in different frequency regimes from visible to mid-infrared. Furthermore, a great potential can be foreseen for exploration of alternative geometries to improve the frequency conversion efficiency.

V. CONCLUSION
In conclusion, we proposed a novel paradigm for designing time-modulated metasurfaces by changing the phase delay in modulation using phase shifters which allows for obtaining a dispersionless phase shift covering 2π span in the forward and backward scattering of generated frequency harmonics, regardless of incident angle and polarization. This principle can be used for designing efficient electrically tunable time-modulated metsurfaces to engineer the wavefronts of generated frequency harmonics. The applicability of the design principle is illustrated for steering and focusing the generated frequency harmonics in a time-modulated metasurface using graphenewrapped silicon microwires which is optically excited in the THz regime and electrically modulated in RF regime. Furthermore, the broadband and wide-angle performance of the metasurface is rigorously verified. This paradigm is a clear departure from the previous electrically tunable metasurface designs as it provides a full control over the phase modulation with a constant amplitude while exhibiting broadband and wide-angle functionality.