Time-varying Huygens' meta-devices for parametric waves

Huygens' metasurfaces have demonstrated almost arbitrary control over the shape of a scattered beam, however, its spatial profile is typically fixed at fabrication time. Dynamic reconfiguration of this beam profile with tunable elements remains challenging, due to the need to maintain the Huygens' condition across the tuning range. In this work, we experimentally demonstrate that a time-varying metadevice which performs frequency conversion can steer transmitted or reflected beams in an almost arbitrary manner, with fully dynamic control. Our time-varying Huygens' metadevice is made of both electric and magnetic meta-atoms with independently controlled modulation, and the phase of this modulation is imprinted on the scattered parametric waves, controlling their shapes and directions. We develop a theory which shows how the scattering directionality, phase and conversion efficiency of sidebands can be manipulated almost arbitrarily. We demonstrate novel effects including all-angle beam steering and frequency-multiplexed functionalities at microwave frequencies around 4 GHz, using varactor diodes as tunable elements. We believe that the concept can be extended to other frequency bands, enabling metasurfaces with arbitrary phase pattern that can be dynamically tuned over the complete 2\pi range.


I. INTRODUCTION
Recent advances in Huygens' metasurfaces and metadevices provide new insight into highly efficient wavefront shaping and scattering manipulation by incorporating electric and magnetic responses [1][2][3][4][5][6][7][8][9][10][11][12][13], extending earlier studies of blazed gratings [14]. While a large amount of work has been done on static Huygens' metasurfaces and meta-devices, advanced applications require that the performance of meta-devices can be tuned to adapt to varying operating conditions or requirements. Ideally, this would mean that the amplitude and phase response of a Huygens' meta-device can be tuned dynamically in an arbitrary fashion. Here, "arbitrary" means the amplitude and phase response can be independently tuned to all the possible states; for example, to tune the phase response over a complete 2π range without changing the amplitude response, or to change the amplitude from zero to maximum while keeping the phase unchanged. While there have been some recent attempts to achieve independent tuning of Huygens' metasurfaces [15], fully arbitrary control remains very challenging to realize in practice, since truly arbitrary tuning requires independent control of not only the electric and magnetic responses, but also of the gain and loss.
To solve this problem, we propose parametric metadevices, where some parameters of the structure can be dynamically modulated. When an electromagnetic wave with central frequency ω 0 interacts with a system having properties modulated with frequency Ω ω 0 , new sideband frequencies ω n = ω 0 + nΩ (n ∈ Z) are generated. The manifestation of this process can be found in * mingkai.liu@anu.edu.au systems of different scales, ranging from the energy level splitting of trapped cold atoms [16,17], Raman scattering of vibrational molecules [18][19][20], Brillouin scattering due to opto-acoustic coupling [21][22][23], to the amplitude and frequency modulation of radio signals [24], and the micro-Doppler effect widely employed in radar sensing [25,26]. The manipulation and detection of these sidebands are of great importance from both fundamental and application points of view. A variety of novel effects explored recently in optical systems with dynamic modulation also rely on sideband control, such as sideband cooling [16,[27][28][29], magnet-free optical isolation [30][31][32][33], and optomechanical interaction in the resolved sideband regime [34,35].
Understanding and harnessing the effects of dynamic modulation in artificial subwavelegnth electromagnetic systems such as resonant particles and metasurfaces, could lead to various ultra-compact tunable devices and novel functionalities via the introduction of an additional degree of freedom -time varying properties [36][37][38][39][40][41][42][43][44]. One important feature of sideband generation in subwavelength systems is that the scattered phases of the sidebands have a gauge freedom and can be controlled by the dynamic modulation process. However, due to the symmetry of dipole scattering, a thin-layer of modulated elements scatters sidebands in both forward and backward directions. In addition, the limited modulation strength available in practice further suppresses the efficiency of energy conversion from the carrier wave to sidebands at the subwavelength scale.
Here, we develop and experimentally verify a theory of time-varying Huygens' meta-devices for parametric waves, which consist of electric and magnetic meta-atoms that can be modulated independently via external stimuli. Compared to conventional static Huygens' metadevices, these time-varying Huygens' meta-devices do not aim to manipulate the carrier wave at frequency ω 0 , but instead efficiently convert the carrier wave into sideband frequencies ω n via dynamic modulation (see Fig. 1). Unlike static Huygens' meta-devices, where the amplitude and phase response to the carrier wave are determined by the intrinsic resonances of the meta-atoms, the sidebands scattered from a time-varying meta-device are fully controlled by the amplitude and phase of modulation.
Since the parities of electric and magnetic dipole radiation are preserved under modulation, the directivity of sideband scattering can be manipulated by the relative modulation phase of electric and magnetic meta-atoms, and high conversion efficiency from the carrier frequency to sidebands can be achieved within a layer of subwavelength thickness. By introducing a gradient of modulation phase along the unit cells of a metasurface, the sidebands can be steered in the transverse direction. We note that some recent studies already provided vivid numerical demonstrations of sideband generation and steering in time-varying Huygens' metasurfaces [45][46][47]; in this study we reveal the physical mechanism and highlight the key components for achieving high conversion efficiency and full spatial control of parametric waves. We demonstrate this concept experimentally in the microwave frequency range, where we measure the sideband scattering of independently modulated electric and magnetic meta-atoms. Using optimized modulation signals, we achieve a high conversion efficiency of over 75% from the carrier wave to the target sidebands and successfully demonstrate various sideband scattering regimes with controlled directionality, including single-sideband unidirectional scattering [ Fig. 1 (top)], double-sideband unidirectional scattering [ Fig. 1 (middle right)], and doublesideband bi-directional scattering [ Fig. 1 (middle left)]. Using a finite time-varying Huygens' meta-device, we demonstrate novel effects including all-angle (360-degree) beam steering and frequency-multiplexed functionalities [ Fig. 1 (bottom)]. These results set the foundation of more advanced Huygens' meta-devices for parametric waves.

II. ANALYTICAL MODEL
A. General description based on the effective impedance model The physics of time-varying Huygens' meta-devices can be described by the boundary conditions of a thin sheet with space-time-dependent electric and magnetic polarizations P(r, t) and M(r, t), dubbed the generalized sheet transition conditions (GSTCs) [48]: n where n is the normal vector of the surface; the superscripts 'i', 'f' and 'b' represent the incident, forward scattered and backward scattered fields, respectively; ' ' denotes the components parallel to the surface. The more general form of GSTCs also includes the normal components of polarizations P ⊥ and M ⊥ (see Appendix A), but in this work, we limit our discussion to a metadevice with polarizations only in the lateral directions (P ⊥ = M ⊥ = 0). When a slow periodic modulation with a fundamental frequency of Ω is introduced in the polarizations, i.e. P = P (r, Ωt) and M = M (r, Ωt), all the timedependent components of Eqs. (1) and (2) can be decomposed in the frequency domain in the form of F (r, t) = +∞ n=−∞ F n (r)e −iωnt , where ω n = ω 0 + nΩ are the sideband frequencies [24], with ω 0 the carrier frequency of the incident field. The electric and magnetic polarizations can be achieved using time-varying meta-atoms with effective Fourier components of electric and magnetic dipole moments p n = P n A, and m n = M n A/µ, respectively, where A is the area of the unit-cell, and µ is the permeability of the medium surrounding the metadevice. The boundary equations (1) and (2) for sideband ω n are then given by δ n0 is the Kronecker delta, which is 1 for n = 0, and 0 otherwise. It is clear that the space-time dependency of the scattered fields on the two sides of a meta-device is determined by the space-time variation of the effective dipole moments; only in the case of static meta-devices (Ω = 0) is the frequency of the scattered fields identical to the incident one. However, the Fourier components of the effective dipole moments are not independent but interact with each other via parametric modulation. In order to describe this parametric process in the meta-device, we extend our previous models for static meta-devices [49,50]. Generally, the interaction between the electric (magnetic) meta-atom at position r and the incident wave can be expressed in a compact form: is the effective impedance of the meta-atom in the array with all the mutual interaction taken into account (see Appendix C for details); I E(M) is the mode amplitude (current amplitude); V E(M) is the effective electromotive force evaluated by the overlap of the incident field and the mode profile [49,50]. When a periodic modulation f (Ωt − ϕ) is introduced to the meta-atoms, Z, I and V become time-dependent, and they can be expanded in the frequency domain by a Fourier series. Representing each meta-atom as a series RLC circuit, the input electromotive force is equal to the total voltage at any point in time, which requires that the coefficients of the Fourier decomposition of Z(t), I(t) and V (t) satisfy the following equation (see Appendix C for details) Here, l = n − m. The subscripts denote the order of the Fourier coefficients of the time-varying Z(t), I(t) and V (t), while the superscript in Z (m) l denotes that the Fourier coefficient is calculated at frequency ω m . ϕ is the relative phase introduced from the modulation, which is effectively a time-delayed replica of the modulation waveform. The dynamics of the modulated meta-atom at position r can then be expressed in a compact matrix form: where I = · · · , I m e imϕ , · · · , I n e inϕ , · · · T , The zeroth order impedance Z (m) 0 describes the linear response of the meta-atom at frequency ω m , where the effect of mutual interaction with other meta-atoms has been taken into account; the (m − n) th order component Z (n) m−n characterizes the frequency conversion from ω m to ω n , which is originated from the dynamic modulation of the self-impedance. For more details regarding the construction of Eq. (7), see Appendix C.
Under normally incident plane wave excitation E i , the Fourier component of the effective electromotive force acting on the electric and magnetic meta-atoms can be expressed as u E,n = −iω n p n /I E,n and u M,n = m n /I M,n are the Fourier components of the normalized effective electric and magnetic dipole moments, which can be defined by the Fourier component of the normalized current distribution j n [49,50] (see Appendix C for more details). Solving Eq. (6) independently for both electric and magnetic meta-atoms, we can find the current amplitudes I E(M),n , then the dipole moments at frequency ω n can be readily found as m n (r) = I M,n (r)u M e inϕM(r) .
Here, ϕ E(M) is the relative phase introduced from the modulation of the electric (magnetic) meta-atom. Once we substitute Eqs. (12) and (13) into Eqs. (3) and (4), it becomes clear that an additional phase of nϕ is introduced to the sideband fields (E n , H n ) via modulation. This feature of time-varying Huygens' meta-devices offers a unique opportunity in controlling wave scattering -since the modulation phase ϕ can be chosen arbitrarily for different unit cells, one can generate any desired spatial phase pattern along the meta-device simply by modifying the local modulation phase ϕ(r). As a result, the phase of each unit cell can be dynamically tuned over a complete 2π range, giving our structure advantage over existing static Huygens' meta-devices that rely on the deliberate overlap and balance of electric and magnetic resonances. Moreover, unlike conventional sideband scattering processes, the functionality of time-varying Huygens' meta-devices can be actively tuned by changing the modulation phase difference between the electric and magnetic meta-atoms.

B. Modulation with uniform amplitude and a linear phase gradient
To give an illustrative example, we consider the simple case of a homogeneous metasurface consisting of identical units, where the response of the meta-atoms is polarization-independent. The incident carrier wave propagates in the normal direction, and the dynamic modulation of the meta-atoms (of the same type) has identical amplitude but a periodic linear phase gradient: ϕ(y) = ϕ 0 + Gy, with G being the spatial frequency of modulation and ϕ 0 being the modulation phase at y = 0. As a result, the solution of the system is a series of Floquet modes characterized by (ω n , β n ) I(y) = · · · , I n (y)e i(nϕ0+βny) , · · · T , where each sideband frequency ω n is scattered with a transverse wave vector β n = β 0 + nG. When the modulation frequency Ω is much smaller than the linewidth of the resonance, the difference between the impedance evaluated at ω l and ω m is negligible, and the effective impedance matrix in Eq. (7) can be simplified as (for details, see Appendix C): Z β,n is the zeroth order effective impedance of Floquet mode (ω n , β n ). For a normally incident carrier wave (β 0 = 0), the sidebands ω n and ω −n have opposite transverse wave vectors: β n = −β −n , i.e. they are deflected towards opposite transverse directions, and their impedances can be approximated as identical: Z β,n = Z β,−n . In the special situation of uniform modulation without any phase gradient (G = 0), all the sidebands are scattered in the normal direction (β n = 0), and the effective impedance can be considered the same for all sidebands: Z β,n = Z β,m , in this special situation, Eq. (15) becomes a Toeplitz matrix [Eq. (C33) in Appendix C]. We further assume that the change of the mode profile is negligible during modulation (i.e. j n = 0 for n = 0), and only the zeroth order terms remain in the effective electromotive force V, as well as in the normalized dipole moments u E and u M of Eq. (10) and (11), then with V n = V δ n0 . These approximations work well in the adiabatic limit of modulation, and the sideband spectra calculated based on the impedance model converge and agree well with the full-wave simulation when the order of truncation increases. See Appendix E for details. In order to give a simple explicit expression, we truncate both the sideband and the modulation to the first order; for purely reactive modulation, we define Z 1 = Z −1 = −iξ, where ξ ∈ R (see Appendix E for more discussion). The effective impedance and the electromotive force at position y are simplified as Z 0 = Z β,0 and Z β = Z β,±1 are the effective impedances for the carrier wave and the first order sidebands. The mode amplitudes can be expressed in a compact form: where σ = Z 0 Z β + 2ξ 2 . Equation (19) shows that while the zeroth order response I 0 is only affected by the modulation amplitude ξ, both the amplitude ξ and the modulation phase ϕ play a key role at the sideband frequencies.
For a particular Floquet mode, the effective impedance of the meta-atom can be characterized with Z = −iX + R (rad) + R (ohm) , where X, R (ohm) and R (rad) correspond to frequency-dependent reactance, ohmic loss and radiative loss, respectively. From the passivity condition, the following relations for electric and magnetic meta-atoms hold for the TE polarized Floquet mode (ω n , β n ) (see Appendix D for details): η is the wave impedance of the surrounding space; k n = ω n /c is the total wave vector and κ n = k 2 n − β 2 n is the longitudinal component. Applying Eq. (19) for both electric and magnetic meta-atoms via Eqs. (12) and (13), and using the relations in Eqs. (10), (11) and (20), we can derive the dipole moments p n and m n , as well as the forward and backward fields E f n and E b n via Eqs. (3) and (4) [for details see Eqs. (A5) and (A6)]. The generalized scattering parameters for TE polarization can be defined as r n = κ n /k n E b n /E i and t n = κ n /k n E f n /E i , where the factor κ n /k n originates from the local power conservation [51] (see Appendix D for more details). The final expressions are given by To have frequency conversion with maximal efficiency (limited only by the fraction of energy dissipated in the meta-atoms), no energy may be transmitted or reflected at the carrier frequency, i.e. t 0 = r 0 = 0. Applying this condition to Eqs. (21) and (22), we get the requirement for the modulation term: , which should be satisfied independently for electric and magnetic meta-atoms. It is clear from Eqs. (21) and (22) that due to the symmetry of the dipole scattered field, perfect conversion is impossible if only electric or magnetic dipole response is employed. For a lossless metasurface, R (ohm) = 0, the condition for complete frequency conversion can be simplified as |ξ| 2 = 1 √ 2 Z * 0 Z β , which can be strictly satisfied when Z 0 and Z β have the same phase angle. From the relations in Eq. (20), this implies that Z E,β = kn κn Z E,0 , Z M,β = κn kn Z M,0 ; as a result, the required impedance modulation |ξ E | 2 = kn 2κn |Z E,0 | 2 , |ξ M | 2 = κn 2kn |Z M,0 | 2 . The corresponding sideband scattering parameters in Eqs. (23) and (24) are then simplified as ψ E(M) = Arg(|Z E(M),0 |/Z E(M),0 ) is the intrinsic linear phase response of the electric (magnetic) meta-atoms and ∆ψ = ψ M − ψ E is the intrinsic phase difference; ∆ϕ = ϕ M −ϕ E is the phase difference between the modulation signals applied to the electric and magnetic metaatoms.
Equations (25) and (26) highlight the capability of time-varying Huygens' meta-devices to control the directionality of sidebands by changing the modulation phase difference between electric and magnetic meta-atoms. For example, for overlapping electric and magnetic resonances (∆ψ = 0), |r ±1 | = 0 and |t ±1 | = 1/ √ 2 when the relative modulation phase ∆ϕ = 0; while |r ±1 | = 1/ √ 2 and |t ±1 | = 0 when ∆ϕ = 180 • . This indicates that a single time-varying Huygens' meta-surface can function as either a transmissive (see the middle-right panel of Fig. 1) or reflective device simply by changing the relative phase ∆ϕ. Another example is when ∆ψ = ∆ϕ = ±90 • , the +1 st and −1 st orders can be well separated in opposite directions: Fig. 1). These effects will be demonstrated experimentally in Sec. IV and highlighted in Figs. 6 and 7.
In contrast to static metasurfaces, overlapped and balanced electric and magnetic resonances at the carrier frequency is not always required in time-varying Huygens' meta-devices, for two reasons: 1) different types of directive sideband scattering require different detuning between the electric and magnetic resonances; 2) the unbalanced linear response of electric and magnetic metaatoms can be compensated in the parametric process by adjusting the amplitudes, phases and waveforms of the dynamic modulation of electric and magnetic metaatoms.
We emphasize that while the above discussion provides some useful insight, the condition for perfect conversion shown above is based on the approximation where both the impedance modulation and the sidebands can be truncated to the first order. In practice, however, this approximation is difficult to achieve since sideband generation is a cascaded process. Moreover, impedance modulation of resonators generally changes both the amplitude and the phase of the carrier wave, and the phase changes nonlinearly with the modulation signal, which inevitably introduces high order sidebands. This effect becomes more pronounced when strong modulation is applied. In order to maximize the energy conversion to the first order sidebands, we need to introduce higher order correction terms in the modulation waveform to suppress the undesirable high order sidebands, as will be shown below.

III. DESIGN OF HUYGENS' UNITS FOR PARAMETRIC WAVES
To validate the concept, we design electric and magnetic meta-atoms working in the microwave regime using full-wave simulation (CST Microwave studio). As a first step, we study a Huygens' unit (a pair of electric and magnetic meta-atoms) in a metallic rectangular waveguide, which is easier to characterize experimentally. We emphasize that while the theoretical discussion above focused on a meta-device excited by a normally-incident plane wave, the impedance model is very general and can be easily extended to other situations. Here, the pair of closely spaced resonators are positioned at the center of the waveguide and are excited by the fundamental waveguiding mode. Although the effective impedance Z of the meta-atom in a rectangular waveguide is different from the one in a periodic array, the whole system is still closely related to the situation discussed in Sec. II B, since in both cases the scattered fields from the electric and magnetic meta-atoms have opposite parities, and the waveguide system can be considered as a special situation discussed in Sec. II B, in which the carrier wave and the sidebands share the same scattering channels (forward and backward in the normal direction). Therefore, we expect that the features discussed in Sec. II B can also be observed in this system. Indeed, we find that even a basic Huygens' unit made of a pair of electric and magnetic meta-atoms is sufficient to demonstrate the effect of directional sideband scattering in conjunction with a high conversion efficiency.
The Huygens' unit consists of an electric and a mag- netic split-ring resonators, which are printed individually on a Rogers RO4003 substrate ( r = 3.5, loss tangent 0.0027, substrate thickness 0.4 mm), and the whole unit is positioned at the center of a WR229 rectangular waveguide that supports a TE 10 mode [see geometry and parameters of meta-atoms in the caption of Fig. 2 (a)]. We use varactor diodes (SMV1405) as the voltagetunable capacitors. Under appropriate DC bias voltages, the electric and magnetic resonances overlap. Figure 2 (b) shows the simulated linear response of the unit at ω 0 /(2π) =3.7 GHz as a function of the DC bias voltages (U E , U M ). There exists a parameter regime around U E /U M ≈ 2, where the reflection approaches zero and the transmission phase changes by around 2π, which is a typical manifestation of the Huygens' condition.
We are interested in the behavior of the unit around the overlapped resonance when dynamic modulation is introduced. We choose the carrier frequency to be ω 0 /(2π) =3.7 GHz, and assume that the modulation frequency Ω is sufficiently low such that the sidebands pro-duced by the dynamic modulation can be calculated in the adiabatic limit (see Appendix F for more discussion), where we approximate the time-varying transmission and reflection signals with voltage-dependent stationary linear responses (see Appendix F for more discussion). In this paper, we only study the behavior in the adiabatic limit of modulation, while a detailed theoretical discussion of modulation beyond the adiabatic limit can be found in Ref. [52].
,amp , and U E(M),offs being the modulation amplitude and the DC offset; is the normalized waveform constructed using a Fourier series. As has been noted in Sec. II B, we introduced high order terms in the modulation waveform to suppress the undesirable higher order sidebands, and we truncate the highest order to N = 8 since higher order terms (N > 8) only bring in negligible improvement. For each set of modulation signals [Ũ E (t),Ũ M (t)], we can obtain the modulated scattering parameters at each time step via a one-to-one mapping of the voltage-dependent stationary linear response from Fig. 2 (b), and the corresponding sidebands can be calculated via Fourier transformation (see Appendix F for more details).
We employ a genetic algorithm to optimize the waveform in order to maximize both the directivity and the power of the chosen sidebands (see Appendix Sec. H for details). As an example, Figs. 3 (a) to (c) depict the normalized spectra for three different types of directive sideband scattering. The spectra are normalized to the total scattered power of all sidebands including the carrier frequency, and more than 80% of the scattered power is converted to the targeted sidebands. The corresponding modulation waveforms and the resulting modulation in transmission and reflection are plotted in Figure 3 (d) to (f). To simplify the optimization procedure, we use the same normalized waveform forŨ E andŨ M . This restriction is not essential, however our empirical tests indicated no improvement in efficiency when allowing these waveforms to differ. The optimized coefficients of the modulating signals are detailed in Table I of Appendix H.
The results of the optimization process reveal that different types of sideband generation require very different modulation waveforms. To achieve a high conversion efficiency for single-sideband forward scattering [ Fig. 3 (a)], the transmission coefficient is ideally modulated linearly over 360 • (see Appendix B for more discussion); the optimized result indeed shows that the phase of transmission changes almost linearly over 360 • , and the required modulation waveform has a highly asymmetric sawtooth-like shape [ Fig. 3 (d)]. In contrast, for double-sideband unidirectional forward scattering [ Fig. 3 (b)], the transmission is ideally amplitude modulated by a purely sinusoidal waveform (see Appendix B), and thus the optimized U E and U M have a sinusoidal-like wave-form and are modulated in phase (∆ϕ ≈ 0) [ Fig. 3 (e)]. Note that the reflection remains low during the modulation cycle and thus the transmission amplitude modulation is achieved via the absorption modulation of the meta-device. The most demanding situation is doublesideband bidirectional scattering [ Fig. 3 (c)], which ideally requires that both the transmission and reflection phases are modulated linearly over 360 • , but with opposite signs of the slopes (see Appendix B). This challenging requirement can be maximally satisfied using electric and magnetic meta-atoms simultaneously since we can find a parameter regime from Fig. 2 (b) where the transmission and reflection phases change dramatically. The final optimized result from Fig. 3 (f) confirms that the transmission phase is indeed modulated over 360 • with a positive slope while the reflection phase is modulated with a negative slope; U E and U M are modulated with a phase-lag ∆ϕ ≈ −90 • , indicating there is an intrinsic phase difference of ∆ψ ≈ −90 • between the electric and magnetic response.
These results are consistent with our previous discussion based on Eqs. (25) and (26), which shows the ∆ϕ required for different types of sideband control. In fact, the optimized sideband spectra can also be reproduced using the impedance model introduced in Sec. II B by taking into account the higher order terms in the impedance matrix shown in Eq. (C33). To evaluate the accuracy of the impedance model, we calculate the Fourier coefficients Z E(M),n using the time-varying scattering coefficients shown in Figs. 3 (d) to (f), and apply these Fourier coefficients in the impedance model to calculate the mode amplitudes I E(M),n as well as the sideband spectra (see Appendix E for details). The results clearly show that as the order of truncation increases, the spectra calculated with the impedance model converge and agree well with the full-wave simulation. In fact, the maximum relative error of the dominant sidebands is already below 5% when we truncate the harmonics to the first order sidebands (see Fig. 10 in Appendix E). Note that different from the simplified case discussed in Eq. (17), the second order terms Z ±2 are also preserved in Eq. (C33), which allows us to calculate the case where the modulation waveform is asymmetric, such as single-sideband conversion.

A. A Huygens' unit in a rectangular waveguide
To demonstrate the dynamic control of sideband scattering, we fabricated the electric and magnetic metaatoms separately on a Rogers RO4003 printed circuit board; varactor diodes were soldered into the additional gaps in the meta-atoms as voltage-tunable elements. The schematic of the microwave setup is shown in Fig. 4. The pair of meta-atoms were fixed on a polystyrene foam holder (with a relative permittivity close to 1) forming a Huygens' unit, and positioned at the center of a rectangular waveguide (see Fig. 4). We first employed the network analyzer to measure the resonance spectrum of the sample and set the DC bias of the arbitrary function generator to tune the two meta-atoms such that their resonances overlap. The resonant frequencies of the fabricated meta-atoms red-shift compared to the simulation, due to the slight difference in substrate permittivity and the connection to the bias network; nevertheless, the resonant frequencies are still in the designed range and can be easily tuned to overlap. When the electric and magnetic resonances are separated, two reflection peaks were observed; as the DC bias voltages (U E,offs , U M,offs ) ≈ (−2.0V, −2.5V ), the resonances of the pair overlap, with a reflection dip observed at around 3.64 GHz and the measured transmission phase varies over 2π, which is a clear evidence for the Huygens condition. We set this point as the initial point for optimization and introduced a microwave CW carrier signal with ω 0 /(2π) = 3.64 GHz from the RF signal generator. Periodic modulation signals with a base frequency Ω/(2π) = 2 MHz were generated by the arbitrary function generator to modulate the meta-atoms; the interaction between the carrier wave and the modulated metaatoms generate sidebands on both sides of the carrier frequency, which were measured by the spectrum analyzer.
The modulation voltage signals were constructed using the same definition of Fourier series shown in Sec. III. The offset voltage U E(M),offs , the amplitude U E(M),amp , the phase ϕ E(M) as well as the coefficients of the waveforms were optimized using a genetic algorithm to max- imize the power conversion to the desirable sidebands (see Appendix H for details). Figure 5 (a), (b) and (c) show the normalized scattering spectra for three different types of directive sideband scattering; the normalized power of the dominant sidebands is around 75% for single-sideband unidirectional scattering and doublesideband bidirectional scattering, and around 84% for double-sideband unidirectional scattering. The corresponding modulation waveforms and the time-varying transmission and reflection signals (measured with I/Q mode of the spectrum analyzer) are shown in Fig. 5 (d), (e) and (f). All the important features of time-varying amplitude and phase match well the theoretical predictions shown in the Fig. 3 (d) to (f). The optimized coefficients of the modulation signals are listed from Table II to Table IV in Appendix H. These results confirm that, with appropriate modulation waveforms, a Huygens' design can achieve various directional sideband scattering with high conversion efficiency.
Further, to demonstrate the capability to tune the directionality dynamically, we kept the optimized waveforms unchanged and only tuned the relative modulation phase of electric and magnetic meta-atoms. Remarkably, for all three types of directive sideband scattering, the directionality of the sideband scattering can be tuned dynamically in between the forward and backward mode, without sacrificing the total conversion efficiency [ Fig. 6 (a), (b) and (c)]. This property is particularly attractive since a time-varying Huygens' meta-device can function as a lens (transmissive mode), a mirror (reflective mode), or something in between, without requiring additional optimization of the waveform but a simple change of the relative phase of modulation. The sinusoidal-like change of the power with respect to ∆ϕ = ϕ M −ϕ E is also consistent with the prediction from Eqs. (25) and (26): |r ±1 | 2 ∝ 1−cos(∆ψ±∆ϕ) and |t ±1 | 2 ∝ 1+cos(∆ψ±∆ϕ).

B. A finite Huygens' meta-device in a parallel-plate waveguide
One advantage of time-varying Huygens' meta-devices is that they allow dynamic complex wavefront control in ways impossible using conventional linear tunable metasurfaces, such as all-angle beam steering and frequencymultiplexed functionalities. To demonstrate these effects, we designed and fabricated an array of Huygens' units working inside a two-dimensional field scanner based on a parallel-plate waveguide, as shown in Fig. 7 (a) (for the detailed designed geometries, see Fig. 11 of Appendix F).
The electric and magnetic meta-atoms are printed on the opposite sides of a single substrate to enhance the mechanical stability.
The array is composed of 16 identical units, each consisting of an electric and a magnetic meta-atom that can be modulated independently. To mimic normally incident plane wave excitation, we generated a collimated beam using a parabolic mirror to direct the signals from a monopole source antenna. We employed four dual channel arbitrary function generators to introduce 8 channels of independent modulation, where four of them were applied to the electric meta-atoms, and the other four used to modulate the magnetic meta-atoms. For the details of the experimental setup, see Fig. 12 and 13 and the corresponding discussion in Appendix G.
To demonstrate dynamic metasurface functionality, we first tuned the bias voltages such that the two resonances overlap around 4 GHz. Then we introduced a CW carrier signal at ω 0 /(2π) =4 GHz and modulated the array with 4 synchronized arbitrary function generators. The modulation waveforms were first optimized to achieve doublesideband forward scattering for sidebands ω ±1 under a uniform modulation phase pattern G1 [see Fig. 7 (h)]; the optimized coefficients of the modulation signals are listed in Table V in Appendix H. For simplicity, we only introduced different modulation phases at different units while keeping the modulation waveforms of all the electric or magnetic meta-atoms identical (see Fig. 7 (h) for the phase patterns). Due to the limited number of independent modulation channels, we only generated phase patterns with 4 discrete phase levels.
Since different sidebands experience different imprinted phases, the time-varying meta-device allows frequency-multiplexed functionalities. When the first order term is dominant in the modulation waveforms, the phase imprinted on sidebands ω ±1 can be considered as conjugated. We confirmed this unique feature by mapping the two-dimensional distribution of the field intensity under the modulation phase patterns G3 and G5, as shown in Fig. 7 (b) to (e). The phase pattern G3 allows large angle beam steering (>60 degrees) of sidebands ω −1 and ω 1 towards different directions, which was further confirmed by the scattering pattern in Fig. 7 (f). A more interesting effect occurred when the phase pattern G5 was applied -the time-varying meta-device functioned as a convex lens for ω 1 and a concave lens for ω −1 at the same time. By properly engineering the modulation waveforms, additional functionalities can be multiplexed into higher order sidebands.
We measured the angular distribution of sideband power via a circular scan around the center of the sample. By introducing different modulation phase patterns, we achieved directional beam steering of the sidebands over the whole 360 • . Figure 7 (f) depicts the measured angle-dependent power for sideband ω −1 . For clarity, we normalize each measurement and plot them in Fig. 7 (g). As has been shown in Fig. 6, adding an additional π difference between the modulation phases ϕ E and ϕ M al- lows the beam to be steered from the forward (marked by circles) to the backward (marked by squares) direction. Such an all-angle directional beam steering highlights the unique capability of time-varying meta-devices as it is a highly nontrivial task for conventional tunable metasurfaces or beam deflectors.
The asymmetry of the measured scattered power in Fig. 7 (f) under conjugated phases originates from the inhomogeneous response of the Huygens' units. This was mainly introduced during the hand-soldering of the varactor diodes and resistors, and it can be overcame by employing a more sophisticated fabrication process. We also note that the peak power reduces at large beamsteering angle, which can be attributed to two main reasons. The first one is due to the fact that we kept the modulation waveforms unchanged during beam steering. In fact, due to the mutual interaction among the units, the effective impedance of the Floquet mode will have a noticeable change at large steering angles compared to the scattering in the normal direction. We expect that additional improvement of the conversion efficiency can be obtained if the modulation waveforms are optimized for each steering angle. The second reason is the reduced directivity due to the decreased effective aperture at large steering angle. In fact, this is a common problem in all array antenna systems [53], and the directivity can be further improved by increasing the sample size (currently ∼ 2λ) and adapting more sophisticated designs of the meta-atoms.

V. DISCUSSION AND CONCLUSION
The concept of time-varying Huygens' meta-devices studied here could provide new tools and possibilities for research and applications that generate and manipulate sidebands. The ability to achieve a high conversion efficiency of sidebands is particularly attractive since it could enable various ultra-compact devices. One potential example is compact isolators. A naive design could consist of a bandpass filter that only transmits ω 0 at the input side, followed by a time-varying Huygens' meta-device that achieves single sideband up-conversion (ω 0 → ω 1 ), and a bandpass filter that only transmits ω 1 at the output side. Since the process of sideband generation is nonreciprocal, a wave propagating in the reverse direction with frequency ω 1 would be further upconverted to ω 2 after passing the time-varying Huygens' meta-device, and be blocked by the filter at the input side.
Another potential application based on the effects shown in Fig. 7 is a compact beam deflector that could steer one or more sidebands at different directions and scan them over the full 4π solid angle. This functionality could benefit the development of new compact radar systems for full-angle and multi-target detection. For example, existing frequency-modulated continuous wave (FMCW) radars widely employed in vehicles for obstacle detection can only see things within a certain angle range, and thus multiple radars are required to cover 360 degrees of view. This solution might become too costly and bulky when it is applied to light-weight platforms such as drones and small robots. A naive vision is to utilise a time-varying metasurface as the transmitter of the radar, from which different sidebands can be steered to different directions. The angular directions of the obstacles all around can be detected simultaneously with a single monopole receiver by identifying the frequencies of the signals bounced back.
The proof-of-concept design presented in this paper shows a promising start, yet further optimization and even new designs are required to facilitate implementation in more practical areas and in a more scalable fashion. We note that while the relatively high quality factor of the meta-atoms employed in the current design allows large modulation with a small voltage (< 10 V), the peak absorption of the system can reach 65% when the electric and magnetic resonances overlap. From simulation we notice that the percentages of the energy lost in the varactor diodes, metal tracks and substrates are around 66%, 20% and 14%, respectively. To reduce the absorption loss, one can design meta-atoms with a lower Q factor by increasing the scattering loss and avoiding the self-resonance of the diode. In addition, one can employ tuning mechanisms that have lower loss, such as MEMS/NEMS [54][55][56][57], tunable capacitors based on ferroelectric thin films [58,59], or transistors.
While the demonstrated dynamic sideband control is limited to the simple case of a one-dimensional array with only 8 independent channels, this is not a fundamental limitation. Using an FPGA or micro-controller to achieve independent modulation over a large two-dimensional array is technically possible, as has been shown recently in the static tuning of a Huygens' meta-device [60]. The extension of the idea to terahertz and even optical frequencies is feasible by using other high-speed modulation mechanisms, such as electro-optical [61][62][63], optomechanical/acousto-optical [64][65][66] and nonlinear optical effects [67][68][69][70].
To conclude, we introduced the concept of timevarying Huygens' meta-devices and studied the sideband generation and manipulation in this type of system. We showed both theoretically and numerically that dynamic modulation of Huygens' meta-devices provides unique opportunities in manipulating parametric wave scattering. Importantly, we designed and fabricated prototype meta-devices working at microwave frequencies, and successfully demonstrated controlled generation and directive scattering in the experiment, with a high conversion efficiency (> 75%) from the carrier wave to the target sidebands. We also demonstrated novel effects that are difficult to achieve with conventional tunable linear metadevices, including all-angle beam steering and frequencymultiplexed functionalities. Our study provides new insights for realizing highly efficient and ultra-compact devices based on dynamic modulation, allowing dynamic tuning of electromagnetic waves in an almost arbitrary way, which should find potential applications in many areas including radar and compressive sensing. D. A. Powell, Y. Zarate and I. V. Shadrivov in the microwave experiment.

Appendix A: Boundary conditions
The physics behind the concept of time-varying Huygens' metasurfaces can be described by the generalized sheet transition conditions (GSTCs) [48]: where n is the normal vector of the surface; '+' and '−' represent field components on the two sides of the surface; '⊥' and ' ' denote the components being normal and parallel to the surface. In the paper, we limit our discussion to a flat meta-device with polarizabilities only in the transverse direction (P ⊥ = M ⊥ = 0), which is described by Eqs. (1) and (2). Without loss of generality, we assume that the meta-device is on the plane z = 0, with P = P x · x + P y · y and M = M x · x + M y · y, and the wave is propagating in the y − z plane. When a slow periodic modulation is introduced, the time-dependent components of Eqs. (1) and (2) can be decomposed in the frequency domain in the form of F (r, t) = +∞ n=−∞ F n (r)e −iωnt , with ω n = ω 0 + nΩ, with ω 0 and Ω being the carrier frequency of the incident field and the modulation base frequency of the polarizations, respectively. The boundary conditions Eqs. (1) and (2) for TE (transverse electric) polarization at frequency ω n are then given by The subscripts 'i', 'b' and 'f' denote the incident, backscattered and forward-scattered field components, respectively; δ n0 is the Kronecker delta function. The conditions for TM (transverse magnetic) polarization can be obtained from duality In the simple case discussed in Sec. II B, the response of the meta-device is isotropic (i.e. polarizationindependent in x and y directions); the carrier wave propagates in the normal direction while the scattered sidebands are Floquet mode with transverse wave vectors β n .
For TE polarized propagating waves, the ratio of the transverse electric and magnetic fields is given by E n /H n = ηk n /κ n . η = µ/ is the wave impedance of the surroundings, k n = ω n /c is the wave vector and κ n = k 2 n − β 2 n is the longitudinal component. By substituting the magnetic field H n with the electric field E n in the boundary equations, we can get the transverse field components in the forward and backward directions as p n = P n A, m n = M n A/µ are the Fourier components of the effective electric and magnetic dipole moments, respectively, where A is the area of the unit-cell and µ is the permeability of the surroundings. Equations (A5) and (A6) are employed to calculate the generalized scattering parameters: r n = κ n /k n E b n /E i and t n = κ n /k n E f n /E i , as shown in Sec. II B.
Appendix B: The ideal condition for different types of sideband control Equations (A5) to (A6) directly link the field scattered from a Huygens' meta-device to the Fourier components of the modulated electric and magnetic dipoles. By assigning the desired scattered fields E n in the forward and backward directions, we can retrieve the required Fourier components of the modulated dipole moments p n and m n based on Eq. (A5) and (A6).
As an example, we examine the required dipole moments for the three types of perfect directive sideband scattering that are discussed in Sec. III, where all the sidebands are scattering in the normal direction (κ = k), as indicated in Fig. 8. In the ideal situation where all the energy from the carrier wave is converted to the desired sidebands: We can get the time varying electric fields n e −inΩt , as well as the required dipole moments p n and m n for the three different types of directive sideband scattering: (a) single-sideband forward scattering which requires where which requires where which requires where The resulting normalized time-varying signals t(t) = E f (t)/E i and r(t) = E b (t)/E i are given in Fig. 8. In the case where only one sideband exists on one side (front or back side) of the meta-device, it requires a pure linear phase modulation of t or/and r with a phase variation of 360 • [see Figs. 8 (a) and (c)]. When both sidebands exist on the same side of the meta-device [see Fig. 8 (b)], a pure sinusoidal amplitude modulation is required. It should be noted that there exist instances where the amplitude of the normalized signal goes above one for double-sideband forward scattering. Such an oscillation is due to the temporal interference of the two sideband waves, but the total energy remains conserved after time averaging over one modulation cycle.
However, we do notice that in the strict adiabatic limit of modulation (Ω/ω 0 → 0), the signals t and r at each time-instance should correspond to the stationary transmission and reflection coefficients of the carrier wave, which should not exceed one at any instance for a purely passive meta-device. Therefore, the required signal in Fig. 8 (b) does imply that at least for adiabatic modulation, balanced gain and loss are required during each modulation cycle in order to achieve a 100% lossless conversion to double sidebands. This is related to the recent findings in static Huygens meta-devices, which showed that balanced spatial-dependent gain and loss are required for perfect beam-steering [13] in the absence of bianisotropy or spatial dispersion [10,11].

Appendix C: The interaction impedance matrix
In general, a mutual interaction exists when the metaatom is positioned in an array or in certain environments such as a waveguide or close to a ground-plane. The linear response of the meta-atom at position r j can be described by the coupled equation: (C1) Z self (r j ) is the self-impedance of the meta-atom at position r j , and Z mut (r j , r i ) is the mutual-impedance that describe the interaction between meta-atoms at r m and r n , and V i (r j ) is the effective electromotive force (input voltage). When dynamic modulation is introduced, Z self (r m ) can be expanded into a matrix to describe the nonlinear parametric process that occurs locally within the metaatom, while Z mut (r m , r n ) can be expanded into a matrix to describe the linear mutual interaction at each individual sideband frequency (as will be shown below).
The effect of mutual interaction can be incorporated with the self-impedance as an effective impedance Z eff , and Eq. (C1) can be simplified as For simplicity, we have written the effective impedance of the electric (magnetic) meta-atom as Z E(M) in the main text, and its interaction with the incident wave can be expressed as Z E(M) I E(M) = V E(M) . Note that if we retrieve the effective impedance of the meta-atom in a periodic array or waveguide via the scattering parameters from full-wave simulations, then the interaction effects are automatically taken into account. The form of V i depends on the incident polarization and the mode profile of the meta-atom. For simplicity, we assume that the incident wave propagates in the y − z plane, the response of the meta-atoms is isotropic in the x − y plane, and the meta-atoms are electrically thin in the z direction. The effective electromotive force acting on the electric and magnetic meta-atoms under a plane wave excitation E i e i(κz+βy) can be expressed explicitly. For TE wave excitation, For TM wave, we have where κ, β and k are the longitudinal, transverse and total wave vectors, respectively. j E (r) and j M (r) are the normalized electric current distributions of the modes on the electric and magnetic meta-atoms. At normal incidence, the expressions are simplified as Eqs. (10) and (11), respectively. u E and u M are the normalized effective electric and magnetic dipole moments of the meta-atoms, which can be defined from the normalized electric current distribution j(r) where the integral is performed over the volume of the meta-atom. Note that since the normalized current j(r) has a unit of m −2 , u E has a unit of m, while u M has a unit of m 2 . The effective impedance can also be defined rigorously based on the mode of the meta-atoms, where the details can be found in our previous studies [49,50].
Below, we derive Eq. (5) in the main text. When dynamic modulation is introduced to the meta-atom at position r j , the effective impedance Z eff , mode amplitude I and effective electromotive force V i become timedependent. Based on Kirchoff's voltage law, the input voltage should equal the total voltage across the series RLC equivalent circuit (see Fig. 9): Here S = 1/C is the elastance, and we use charge Q since it yields purely differential equations rather than integrodifferential equations. V i is the input voltage generated by a particular set of polarization and mode, represented by one of Eqs. (C3) to (C6). For illustration, we use When the temporal modulation of the impedance is periodic in time with a frequency Ω ω, each time varying effective circuit parameter in the above equations can be expanded in a Fourier series of the form: F (t) = +∞ l=−∞ F l e −il(Ωt−ϕ) , and the charge can also be expanded as Q(t) = e −iωt +∞ l=−∞ Q l e −il(Ωt−ϕ) = +∞ l=−∞ Q l e −iω l t e ilϕ . ϕ is the modulation phase corresponding to a constant time-offset for the modulation waveform. Then we have For electric meta-atoms excited by a normally incident field Satisfying Eq. (C9) at every instant in time requires that the Fourier coefficients satisfy the following equation Orthogonality of the exponential functions implies that equivalent terms on each side of Eq. (C17) must be equal. This yields the following equation for each order n: where l = n − m. We further define which is Eq. (5) in the main text. Note that for brevity, we have omitted the notation r j in the above equations. Now we bring the position notation r j back and write the coefficients of Eq. (C19) in a matrix form: I(r j ) = · · · , I m (r j )e imϕj , · · · , I n (r j )e inϕj , · · · T , (C22) V(r j ) = · · · , V m (r j )e imϕj , · · · , V n (r j )e inϕj , · · · T . (C23) ϕ j is the modulation phase at r j . It is important to note that while we have incorporated the effect of mutual interaction in the effective impedance ← → Z eff , we have assumed that only the modulation of self-impedance contributes to the higher order terms Z (m) l (l = 0), i.e.
Z (m) self (r j ) is the zeroth order self-impedance of the meta-atom at position r j , evaluated at frequency ω m .
The mutual impedance matrix describes the linear interaction at each sideband frequency, and thus it is a diagonal matrix with only the zeroth order terms Therefore, the value of the zeroth order effective impedance Z m 0 in Eq. (C21) is determined by not only the zeroth order self-impedance but also the mutual impedances and the mode amplitudes of all other interacting meta-atoms at sideband frequency ω m . In general, Z if m = n, except for some special cases.
When the modulation frequency Ω is much smaller than the linewidth of resonance of the meta-atoms, the difference between the impedance elements of the same order evaluated at different sideband frequencies ω l and ω m becomes negligible. We can simplify the notation as mut , and the impedance matrices Eq. (C25) and (C26) can be simplified as Toeplitz matrices: · · · Z self (r j ) · · · Z n (r j )e inϕj · · · . . . . . . . . .
For the special example discussed in Sec. II B, where the dynamic modulation of meta-atoms has an identical modulation amplitude but a periodic linear phase gradient along y direction: ϕ(y) = ϕ 0 +Gy, the solution of the system Eq. (C22) should be a series of Floquet modes: where I(y j ) = · · · , I n (y j )e i(nϕ0+βnyj ) , · · · T , G is the spatial frequency of modulation, β n = β 0 + nG is the transverse wave vector of the Floquet mode, and ϕ 0 is the modulation phase at y = 0. Applying Eqs. (C27) to (C31) in Eq. (C24), it becomes clear that the effective impedance matrix of the series of Floquet modes is given by where the higher order components Z n are generated directly from the modulation of self-impedance, while Z β,m (y j ) = Z self (y j ) + i Z mut (y j , y i )e imG(yj −yi) is the zeroth order effective impedance of Floquet mode (ω m , β m ). With a normally incident carrier wave (β 0 = 0), we have Z β,m = Z β,−m , and the effective impedance matrix is given by Eq. (15). Only in the special case where the sidebands are also scattered in the normal direction (G = 0), the effective impedance matrix becomes a Toeplitz matrix: Generally, the Fourier coefficient Z n and Z −n are complex but not conjugated. In the special situation where only reactive modulation exists, Z 1 = −Z * −1 = −iξ, and the phase of the complex value ξ is determined by the time delay of the modulation waveform. Since we already introduce the phase parameter ϕ to describe the time delay of the modulation signal, it is convenient to set the phase of ξ to zero in Eq. (17) such that ξ ∈ R.
Appendix D: The relation between radiative loss and normalized dipole moments Below we derive the relation between the radiative loss term R (rad) and the normalized dipole moments for the Floquet mode (ω n , β n ) under TE polarization, as shown in Eq. (20) of the main text: Note that the above relations are only valid for Floquet mode (ω n , β n ) under TE polarization. These intrinsic relations are determined by the current distributions of the modes, and thus should hold regardless of the modulation condition. Therefore, it is convenient to derive the relations in an unmodulated meta-device consisting of identical meta-atoms in each unit, where both the excitation and the scattering are propagating plane waves with a transverse wave vector β n . For a lossless metadevice (R (ohm) = 0), the fields in the forward and backward directions are given by where E p and E m are the scattered plane waves generated by the array of electric and magnetic meta-atoms, respectively. From the passivity condition, the power should be conserved: From Eqs. (D2), (D3) and (D4), we have the following relation: This relation should hold regardless of the ratio of E p and E m , i.e. it should also be satisfied even for an array of pure electric or pure magnetic meta-atoms, and thus requires the electric fields generated by electric and magnetic meta-atoms satisfy the following relations individually: Substitute with the following relations: and use the identity for lossless meta-atoms Re(1/Z) = Re(Z * )/|Z| 2 = R (rad) /|Z| 2 , we get Eq. (20). error of the impedance model when comparing to the sideband spectra calculated with full-wave simulation: |ζ full − ζ imp |/ζ full ; for clarity, we only show the relative error of the dominant sidebands in each type of modulation. It is clear that as the order of truncation increases in the impedance model, the result converges and approach to the full-wave calculation. Surprisingly, the maximum relative error in the three cases studied is already below 5% even when we truncate to the first order sideband, as shown in Fig. 10 (b) for the case of double-sideband forward scattering. Note that different from the simplified case discussed in Eq. (17), the second order terms Z ±2 are also preserved in the impedance matrix Eq. (C33) when we truncate the sidebands to N = 1. The introduction of higher order impedance terms Z n (|n| > 2) allows us to calculate the case where the modulation waveform is asymmetric, as is required for the case of single-sideband conversion.

Appendix F: Numerical simulation
To validate the practicality of the concept, we design electric and magnetic meta-atoms working in the microwave regime using full wave simulation (CST Microwave studio).
For a single Huygens' unit in a rectangular waveguide, we employ electric boundary condition on the waveguide sidewalls and open boundary condition otherwise. The meta-atoms are positioned in the middle of the waveguide and excited with the fundamental mode of the rectangular waveguide.
For the Huygens' array, we re-optimize the designs in the parallel-plate waveguide. We simulate the structure under periodic boundary conditions along the direction of the array and perfect electric boundaries on the top and bottom waveguide walls. The unit is excited with the fundamental TEM mode. The detailed geometries can be found in Fig. 11. To enhance the mechanical stability, the electric and magnetic meta-atoms are designed on the opposite sides of a 0.8 mm thick substrate (Rogers RO4003), and the center-to-center distance be-tween neighboring units is 10 mm. High-resistance (10 kω) resistors are employed for the bias-lines to avoid the out-coupling of microwave signals to the external modulators.
To calculate the sidebands in the adiabatic limit, we first run a full-wave simulation of the system by replacing the diodes with discrete ports; the full impedance matrix of the system can be extracted and used to perform a circuit simulation, where we define the property of the varactor diode based on its SPICE model. The transmission and reflection coefficients under different bias voltage U E and U M are calculated by a parameter scan in the circuit simulation of CST. The simulated response of the array [ Fig. 11 (b)] has the same feature as the one of a single unit in the waveguide [ Fig. 2 (b)]. The time-dependent transmission and reflection under periodic modulation of U E (t) and U M (t) are found via a one-to-one voltagescattering parameter mapping; finally, a Fourier transformation of the time-dependent signals gives the information of sidebands in both forward and backward directions.
This adiabatic approximation is valid when the modulation frequency Ω, the range of resonant frequency modulation ∆ and the linewidth of the resonance γ satisfy ∆Ω γ 2 [52]. In our studied system, Ω/γ ≈ 0.009, ∆/γ < 1.7, therefore the adiabatic approximation is valid. This approximation could provide more physical insight and higher efficiency during optimization compared to first-principle calculation methods such as FDTD, which become very inefficient when critical time scales vary by several orders of magnitude.
For the scattering of a finite array under different modulation phase patterns, due to the limited numerical simulation capacity, we can only perform a linear simulation to give a reference, as we assume that the mode profiles of the meta-atoms will not have a substantial difference at the sidebands and the carrier frequency when the modulation frequency satisfies the adiabatic approximation. To do so, we perform a linear simulation of a periodic array and calculate the far-field scattering of a single unitcell at the frequency with minimal reflection, which gives directional scattering. Then we assign different phase patterns and simulate the scattering power pattern of an array of 16 units to mimic the situation shown in Fig. 11. The results have a reasonably good agreement with the measured ones in Fig. 7 (f), even though in the simulation we assume that the scattering amplitude from each unit is identical and does not change under different modulation phase patterns [see Fig. 11 (c)].

Appendix G: Experimental measurement
To confirm the phase control of sidebands can also be employed in dynamic beam steering, we fabricated an array composed of 16 Huygens' units and positioned the array inside a two dimensional field scanner based on parallel plate waveguide. Each unit has 2 ports (4 patches) connected to the external modulation, as can be seen from Fig. 11 -two patches connected to the ground, one patch for the electric meta-atom and one patch for the magnetic meta-atom, and all 64 patches of the array were independently connected with thin enamel wires that connect to the arbitrary function generators.
We identified the resonances of the meta-atoms by measuring the scattering in the forward direction with the network analyzer, which shows two dips in the spectrum. We tuned the bias voltages such that the two resonances overlap around 4 GHz when the bias voltages are U E = U M = −2V , which were chosen as the starting point of optimization. The array was then excited with a CW signal at 4 GHz (λ ≈ 75mm), which was a collimated beam generated by a parabolic mirror with a monopole source antenna positioned around 315 mm (∼ 4.2λ) away from the array center (see Figs. 12 and 13). We employed four double-channel arbitrary function generators to provide eight channels of modulation waveforms, four of them were applied on the electric meta-atoms and the other four were used to modulate the magnetic meta-atoms. Due to the limited number of modulation channels available, we only generated modulation phase patterns with 4 discrete phase levels. The patterns were reconfigured by rewiring the connections on a breadboard circuit.
The optimization was performed using the scattered fields measured in the forward (θ = 0 • ) and backward (θ = 180 • ) directions, under the uniform modulation phase pattern of G1. After optimization, a 2D line scan over a circle centered around the Huygens array was performed to get the scattering pattern. Due to the limited scanning range, the radius of the circle is limited to 2λ. Since the sidebands scattered away from the meta- Schematic of the arithmetic-crossover operation. Two selected parental data sets (chromosomes) from generation N are hybridized to generate a new data set for generation N+1 via convex combination. g ∈ (0, 1) is a random number.
surface decay with distance, the effect of reflection from the parabolic mirror and the non-perfect absorbers at the boundaries is relatively small. The measured level of sidelobes for ω ±1 under the phase pattern G1 is around -10 dB. However, for the carrier wave, since the incident beam width is slightly larger than the sample size, the scattered field from the metasurface, the reflections from the parabolic mirror and the non-perfect absorbers at the boundaries interfere with the incident field, which makes an accurate estimation of the spurious scattering quite complicated. Therefore, in the optimization of the array, we only aimed at achieving directive scattering for sidebands ω ±1 and maximizing their power ratio over all other sidebands, except for the carrier wave scattered in the backward direction. From the two-dimensional distribution of the measured power under the modulation phase pattern G1, we estimate the levels of the scattered power over the incident one to be around 24% for sidebands ω ±1 and around 3% for all other higher order sidebands; the rest are the spurious scattering of the carrier wave and the loss due to absorption. is the normalized waveform constructed using a Fourier series, where we take N=8 as the highest order term. For brevity, we will neglect the subscript of "E(M)" below and write the sets of coefficients a (l) and b (l) as {a} and {b}. The problem we need to solve is to find the optimal set of coefficients, i.e. ({a}, {b}, U amp , U offs , ϕ), in order to achieve the desired types of sideband scattering.
The flow of optimization, as depicted in Fig. 14 (a), includes many rounds (generations) of operations. In each round, a total number of N sets of coefficients (genes) are generated, and they are used to construct N different modulation signals, which gives N different values of the objective function (O.F.). Here we define the objective function as O.F. = P T × B P ; P T is the total scattered power of the desired sidebands, and B P = min(ζ n /ζ m , ζ m /ζ n ) is a function that gives maximum value when the power of the two desired sidebands are equal. This objective function is employed to maximize the conversion efficiency while maintaining a balanced power when two sidebands are involved.
The best M genes that give the highest values of O.F. are selected as parents to produce C 2 M = M (M − 1)/2 new genes (offspring) for the next generation via crossovers and mutation. There are many different approaches for crossover operations [73], such as one-point-crossovers, multi-point-crossovers, uniformcrossovers. Here we employ arithmetic-crossovers that produce offspring via a convex combination of the genes from two parents, as depicted in Fig. 14 (b). After the crossover, mutation is performed by introducing small random values into both the new genes and the old genes to ensure the diversity; then, genes with the lowest values of O.F. are replaced by the best parents (without mutation) and the best "ancestors" from previous generations in the new generation. This iterative process continues until the stopping criteria (e.g., the desired value of O.F.) are met. In both the simulation and the experiment, we employed 40 genes in each generation, and the results converge well after around 30 to 40 generations.
The optimized coefficients for the three types of sideband scattering shown in Fig. 3 (d) to (f) are summarized in Table I, where we used the same normalized waveform for electric and magnetic meta-atoms. For doublesideband scattering, we employed symmetric waveforms such that b (l) = 0.
In the experiment, we set all coefficients to free param-eters, and they are summarized in Table II, III , IV and V, respectively. We note that ϕ is not an independent parameter since the effect of time-delay can also be accomplished by changing all a (l) and b (l) ; however, it is much more convenient to use ϕ, in particular when we change the relative time-delay of the two modulation signals while maintaining the waveforms unchanged, as has been shown in Fig. 6 and 7.