Rigorous space-time coupled-wave analysis for patterned surfaces with temporal permittivity modulation [Invited]

Sandeep Inampudi; Mohammad Mahdi Salary; Samad Jafar-Zanjani; Hossein Mosallaei

doi:10.1364/OME.9.000162

1. Introduction

The quest for extreme light manipulation techniques utilizing metamaterials and metasurfaces over the past decade has branched into the discovery of new classes of materials with active-tunable optical properties. The real-time control of the optical response of these active materials has been achieved by exploiting different physical effects such as thermal phase transitions [1–3] and carrier-induced changes [4–8]. Even though very complex, starting from base materials such as silicon, indium tin oxide (ITO), and graphene, the state-of-the-art fabrication techniques have made it possible to realize nano-scale to micro-scale patterns of active tunable materials in all forms including dielectric [9–11], plasmonic [4, 5, 12], and two-dimensional materials [7, 13, 14].

Besides fabrication challenges, the subwavelength patterned active tunable materials also confront the ability of the computational tools to study and optimize the light propagation properties through multi-layered devices. Although standard time-domain methods, such as the finite-difference time-domain (FDTD), can be employed for this purpose, the simulation of realistic time-modulated patterned surfaces can become prohibitively challenging, due to the multiscale features in both spatial and temporal domains. In particular, the great difference in time-scales of modulation and optical frequencies requires a significantly large number of optical cycles in the simulation to accurately capture the steady-state response, resulting in expensive computational cost. Approximate frequency-domain methods based on coupled-mode theory [15–19] and equivalent circuit models [20] have been previously adopted, which are not exact, due to their underlying assumptions, and have limited scope of application. Even though accurate multi-frequency couplings based frequency-domain techniques, such as the finite-difference frequency-domain (FDFD) [21, 22] and transition matrix [23], have been demonstrated recently, the former requires computationally expensive discretization of whole spatial domain and the latter is limited to the cylindrical scatterers with time-modulated boundaries.

Space-time-modulated patterned surfaces, metasurfaces, whose parameters vary in both spatial and temporal domains–rendering a four-dimensional design space–are currently being considered as the next frontiers for extreme light manipulation. Recent studies have demonstrated that introducing spatio-temporal modulation in a metasurface can lift the reciprocity constraint and break the time-reversal symmetry [24–26], for device applications such as isolators [11, 15, 27–31]. In addition, space-time-modulated metasurfaces exhibit frequency-mixing property and generate higher-order frequency harmonics which can extend the degree of light manipulation through space-time photonic transitions [23, 32–35] and offer several new opportunities including realization of temporal photonic crystals and transmission lines with tunable bandgaps [36–40], pure frequency mixing [41], Doppler cloaking of moving targets [42], orbital angular momentum radiation [43], dynamic beam modulation [23], stopping [16,44], localization, and time-reversal of light [18,45,46].

In order to circumvent the above discussed computational challenges in design and optimization of space-time-modulated metasurfaces for tailored performance, here we present a fast and accurate semi-analytical technique called rigorous space-time coupled-wave analysis (RSTWA) developed on the grounds of the rigorous coupled wave analysis (RCWA) [47–49] method. RCWA is a well-established plane-wave-expansion-based Maxwell’s equations solver exclusively developed to compute light propagation properties through patterns with spatial periodicity. Since by nature periodic structures diffract an incident plane wave only to a discrete set of plane wave components (‘diffraction orders’) characterized by the grating equation, the RCWA formulation relies on computing the propagation constants, field profiles, and amplitude relationships of these limited set of plane wave components through the structures. Complete convergence of solutions even in complex multi-patterned layered structures is typically observed with the inclusion of few tens of diffraction orders. The absence of requirement for any spatial or temporal meshing in combination with fast and accurate eigen value solvers, such as LAPACK [50], enable RCWA to outperform other numerical solvers such as FDTD [51], finite element method (FEM) [52], and integral equation (IE) [53] based methods in case of periodic structures with subwavelength scale patterns. In addition, RCWA offers flexibility in treating layered structures of arbitrarily-shaped gratings and has been successfully applied for simulating the performance of static metasurfaces with one-dimensional and two-dimensional periodicities [54–56].

To captivate the above advantages, we extend the extremely functional RCWA technique to RSTWA for fast and accurate modeling of light propagation through material patterns with both spatial and temporal modulation of optical properties. More specifically, we begin with the layout of RCWA formulation for two-dimensional periodic structures (say in the xy− plane), and exchange one of the spatial dimensions with time (say xt− plane), as shown in Fig.1. The otherwise two-dimensional scattering of light into diffraction orders will therefore reduce to one-dimension, but with both spatial and temporal harmonics. For example, in place of Fourier components of a two dimensional periodic pattern in the xy−plane, characterized by q_x_;_m = m2π/Λ_x and q_y_;_n = n2π/Λ_y, we utilize the equivalent space-time Fourier harmonics characterized by q_x_;_m = m2π/Λ_x and q_t_;_n = n2π/Λ_t and given by,

ε^{(m, n)} = \frac{1}{Λ_{x}} \frac{1}{Λ_{t}} \int_{0}^{Λ_{x}} \int_{0}^{Λ_{t}} d x d t ϵ (x, t) \exp (i q_{x; m} x - i q_{t; n} t)

In the following section, utilizing the above equation and similar exchange of space to time dimension in electric and magnetic fields, we present detailed formulation for RSTWA technique to multi-patterned layered structures. We validate the presented formulation by comparing the outcomes against FDTD simulations. We further demonstrate the capabilities of the developed RSTWA algorithm using few representative application examples for parametric characterization and device simulations of spatio-temporally modulated metasurfaces made of silicon nanobars and graphene micro-ribbons. The former case leads to anomalous bending of higher-order frequency harmonics, while the latter leads to the generation of dynamic beams.

Fig. 1 Space-time diffraction from a spatio-temporally modulated grating.

Download Full Size | PDF

2. Matrix formulation for rigorous space-time coupled-wave analysis (RSTWA)

The aim of RSTWA is to provide a scattering matrix type relationship between the amplitudes of incident, reflected, and transmitted diffraction orders (spatio-temporal harmonics) through the structure. A general schematic of the problem under consideration with a simple diffraction grating made of active tunable material is shown in Fig. 1. Any arbitrary shaped unit cell being periodic in terms of both x and t while having a layer-by-layer variation along z can be considered.

A linearly polarized plane wave incidence is assumed from one side of the system with an oblique incident angle (θ_in) and frequency f₀ (also known as the central frequency). The spatial period length of the unit cell along the x-axis is considered as Λ_x, and the temporal period length is taken as Λ_t (= 2π/Ω), where Ω is the modulation frequency. As discussed above, the periodic nature of the system in spatial and temporal coordinates, ‘diffracts’ an incident plane wave into a discrete set of spatial and temporal harmonics (also plane waves) characterized by frequency and wavevector components ω_n and k_x_;_m respectively. The quantities ω_n and k_x_;_m are given by the grating equation,

ω_{n} = ω_{0} + n Ω, a n d k_{x, m} = k_{x, 0} + m 2 π / Λ_{x},

where, ω₀ represents the angular frequency of the incident (central) wave given by (ω₀ = 2π f₀)and k_x_;0 represents the corresponding incident angle as k_x_;0 = (ω₀/c) sin θ_in. m, n are integers spanning from −∞ to ∞ that represent the index of the diffraction orders. Each ordered pair (m, n) corresponds to a coupled spatio-temporal harmonic wave characterized by a frequency (ω_n) and a wavevector component (k_x_;_m). Therefore, the corresponding total electric and magnetic field (x and y) components in a given layer l can be compactly written as,

{E, H}_{{x, y}; l} (\vec{r}) = \sum_{m = - \infty}^{\infty} \sum_{n = - \infty}^{\infty} e^{i k_{x, m} x} e^{- i ω_{n} t} {ℰ, ℋ}_{{x, y}; l}^{(m, n)} (z)

where

{ℰ, ℋ}_{{x, y}; l}^{(m, n)}

represent the amplitude of modal profiles of the corresponding field components in wavevector domain. The summand of the above equation in an xt−plane at a given z-position represents the field distribution of a spatio-temporal harmonic (characterized by k_x_;_m and ω_n) with Floquet-periodicity boundary condition at the unit cell boundaries (Λ_x and Λ_t). The modal amplitudes of electric

(ℰ_{{x, y}, l}^{(m, n)})

and magnetic

(ℋ_{{x, y}, l}^{(m, n)})

components can be further expressed as a sum of forward and backward propagating, transverse electric (TE or s) and transverse magnetic (TM or p) components, in matrix form as,

[\begin{matrix} ℰ_{x; l}^{(α)} (z) \\ ℰ_{y; l}^{(α)} (z) \\ ℋ_{x; l}^{(α)} (z) \\ ℋ_{y; l}^{(α)} (z) \end{matrix}] = [\begin{matrix} W_{l} & W_{l} \\ V_{l} & - V_{l} \end{matrix}] [\begin{matrix} a_{l; s}^{(α) +} e^{i k_{z, l}^{(α)} (z - z_{l - 1})} \\ a_{l; p}^{(α) +} e^{i k_{z, l}^{(α)} (z - z_{l - 1})} \\ a_{l; s}^{(α) +} e^{i k_{z, l}^{(α)} (z_{l - 1})} \\ a_{l; p}^{(α) +} e^{i k_{z, l}^{(α)} (z_{l - 1})} \end{matrix}]

where the column vectors

a_{l; s, p}^{α \pm}

represent the amplitudes of forward (+) and backward (−) traveling (s and p polarized) waves with respect to propagation direction. Note that the double index {(m, n)} is replaced with single index (α) that spans over all possible combinations of (m, n) for numerical implementation purpose.

The motive of the RSTWA technique is to establish a scattering matrix type relationship between the amplitude vectors a_l of the incident $(a_{1}^{+})$ , reflected $(a_{1}^{-})$ , and transmitted $(a_{N}^{+})$ spatio-temporal harmonics of the system.

In case of homogeneous layers with constant permittivity (ϵ_l), by solving Maxwell’s Equations for electric and magnetic fields in unbounded medium with relative permittivity (ϵ_l) and magnetic permeability (µ_l = 1), the matrices W_l and V_l can be expressed as block diagonal matrices defined as,

W_{l} = [\begin{matrix} O & \frac{1}{ϵ_{0} ϵ_{l}} W^{- 1} K_{z; l} \\ I & O \end{matrix}]; V_{l} = [\begin{matrix} - \frac{1}{μ_{0} μ_{l}} W^{- 1} K_{z; l} & O \\ O & I \end{matrix}]

where O and I represent the null matrix and identity matrix of respective dimensions,

W

and K_z_;_l represent diagonal matrices whose elements are given by the angular frequency and z–component of wavevector respectively, defined as:

W^{(α, α)} = ω_{α}; K_{z; l}^{(α, α)} = \sqrt{{(ω_{α} / c)}^{2} ϵ_{l} - {(k_{x, α})}^{2}} .

In case of patterned layers, by taking Fourier transform of Maxwell’s equations and utilizing Fourier convolution theorem on the product $ϵ (x, t) \vec{E}$ (as detailed in appendix I), an eigenvalue equation of tangential electric and magnetic field profiles of the modes can be derived as,

[\begin{matrix} {\dot{ℰ}}_{x; l}^{(α)} (z) \\ {\dot{ℰ}}_{y; l}^{(α)} (z) \\ {\dot{ℋ}}_{x; l}^{(α)} (z) \\ {\dot{ℋ}}_{y; l}^{(α)} (z) \end{matrix}] = [\begin{matrix} O & A_{l}^{(1)} \\ A_{l}^{(2)} & O \end{matrix}] [\begin{matrix} ℰ_{x; l}^{(α)} (z) \\ ℰ_{y; l}^{(α)} (z) \\ ℋ_{x; l}^{(α)} (z) \\ ℋ_{y; l}^{(α)} (z) \end{matrix}]

where

\begin{array}{l} A_{l}^{(1)} = [\begin{matrix} O & i (μ_{0} μ_{r} W - \frac{1}{ϵ_{0}} K_{x} Ξ_{l}^{- 1} W^{- 1} K_{x}) \\ - i μ_{0} μ_{r} W & O \end{matrix}]; \\ A_{l}^{(2)} = [\begin{matrix} O & i (\frac{1}{μ_{0} μ_{r}} K_{x} W^{- 1} K_{x} - ϵ_{0} W Ξ_{l}) \\ i ϵ_{0} W Ξ_{l} & O \end{matrix}] \end{array}

in which K_x is a diagonal matrix whose elements are defined as

K_{x}^{(α, α)} = k_{x, α}

. The matrix Ξ_l is the matrix of Fourier coefficients of the permittivity distribution function of the unit cell ϵ_l(x, t) defined as (further details of derivation are given in appendix I),

Ξ^{(α, β)} = \frac{1}{Λ_{x} Λ_{t}} \int_{0}^{Λ_{x}} \int_{0}^{Λ_{t}} d x d t ϵ_{l} (x, t, ω_{α - β}) \exp {i 2 π [(m_{α} - m_{β}) x / Λ_{x} - (n_{α} - n_{β}) t / Λ_{t}]},

where α and β are index terms that span over all spatio-temporal harmonics, and integers (m_α, n_α)and (m_β, n_β) represent corresponding diffraction order indices. Note that spatio-temporal permittivity function in the integrand also contains a frequency dependence (ω_α₋_β) to include any possible material dispersion. Alternately, the matrix elements of Ξ can be obtained utilizing the pre-computed Fourier coefficients of the space-time unit cell in Eq. (1) as

Ξ^{(α, β)} = ε^{(m_{α} - m_{β}, n_{α} - n_{β})}

The accent dot on top of the field components represents a spatial derivative along the z-direction (∂/∂z). To reduce the computational cost, Eq. (7) can further be reduced to,

[\begin{matrix} {\ddot{ℰ}}_{x; l}^{(α)} (z) \\ {\ddot{ℰ}}_{y; l}^{(α)} (z) \end{matrix}] = A_{l}^{(1)} A_{l}^{(2)} [\begin{matrix} ℰ_{x; l}^{(α)} (z) \\ ℰ_{y; l}^{(α)} (z) \end{matrix}] .

Since the only z dependence of the fields is in the form of the phase factor $e^{\pm i k_{z, l} z}$ , the second derivatives in the above equation will result into the formation of an eigenvalue equation of electric field components as,

(A_{l}^{(1)} A_{l}^{(2)} + K_{z; l}^{2} I) [\begin{matrix} ℰ_{x; l}^{(α)} (z) \\ ℰ_{y; l}^{(α)} (z) \end{matrix}] = O .

Equation (11) shows that the wavevector components (k_z_;_l_,_α) of modes in the patterned layer can be computed as $- \sqrt{eigenvalues}$ of the matrix product $A_{l}^{(1)} A_{l}^{(2)}$ and the electric field profile matrix (similar to W_l in Eq. (4)) can be computed as matrix of eigen vectors of the products $A_{l}^{(1)} A_{l}^{(2)}$ . Consequently, using Eq. (7) the magnetic field profile matrix (V_l) of the patterned layer can be computed as,

V_{l} = A_{l}^{(2)} W_{l} K_{z; l}^{- 1}

with modal field profiles matrices (W_l) and (V_l) of all homogeneous and patterned layers of the system computed and stored in memory, the boundary conditions at the l^th interface of the system shared by the l^th and l + 1^th layer can be expressed as,

[\begin{matrix} ℰ_{x; l + 1}^{(α)} (z_{l}) \\ ℰ_{y; l + 1}^{(α)} (z_{l}) \\ ℋ_{x; l + 1}^{(α)} (z_{l}) \\ ℋ_{y; l + 1}^{(α)} (z_{l}) \end{matrix}] = [\begin{matrix} ℰ_{x; l}^{(α)} (z_{l}) \\ ℰ_{y; l}^{(α)} (z_{l}) \\ ℋ_{x; l}^{(α)} (z_{l}) \\ ℋ_{y; l}^{(α)} (z_{l}) \end{matrix}] .

By substituting Eq. (4) in Eq. (13) and after few algebraic steps, one would arrive at the scattering matrix equation,

[\begin{matrix} a_{l}^{(α) -} \\ a_{l + 1}^{(α) +} \end{matrix}] = [\begin{matrix} r_{l}^{+} & r_{l}^{-} \\ t_{l}^{+} & t_{l}^{-} \end{matrix}] [\begin{matrix} a_{l}^{(α) +} e^{i k_{z, l}^{(α)} (z_{l} - z_{l - 1})} \\ a_{l + 1}^{(α) -} e^{i k_{z, l}^{(α)} (z_{l + 1} - z_{l})} \end{matrix}]

of the interface relating the amplitudes of incoming waves towards the interface to the out going waves away from the interface. The definitions of

r_{l}^{\pm}

and

t_{l}^{\pm}

can be derived as,

\begin{array}{l} r_{l}^{+} = - {(a + b)}^{- 1} (a - b), & t_{l}^{-} = 2 {(a + b)}^{- 1} \\ t_{l}^{+} = 0.5 (a + b) + 0.5 (a - b) r_{l}^{+}, & r_{l}^{-} = (a - b) {(a + b)}^{- 1} \end{array}

where

a = W_{l + 1}^{- 1} W_{l}

and

b = V_{l + 1}^{- 1} V_{l}

The matrices

r_{l}^{\pm}

and

t_{l}^{\pm}

epresent the reflection and ransmission coefficient matrices of waves from both sides of the interface when the interface is isolated from other layers of the system. In the presence of all the interfaces at their respective positions (z_l), the net reflection and transmission matrices (R_l and T_l) of each interface, considering the multiple reflections between layers, can be derived as an iterative set of relations as defined below,

\begin{array}{l} Y_{l} = {(I - r_{l}^{-} e^{i K_{z, l + 1} (z_{l + 1} - z_{l})} R_{l + 1})}^{- 1} (t_{l}^{+} e^{i K_{z, l + 1} (z_{l} - z_{l - 1})}) \\ R_{l} = r_{l}^{+} e^{i K_{z, l} (z_{l} - z_{l - 1})} + t_{l}^{-} e^{i K_{z, l + 1} (z_{l + 1} - z_{l})} R_{l + 1} T_{l} \end{array}

Once the net reflection and transmission matrices of all interfaces are computed, the amplitudes f each plane wave component in each layer can be obtained by using the relations $a_{l + 1}^{(α) +} = T_{l} a_{l}^{(α) +}$ and $a_{l}^{(α) -} = R_{l} a_{l}^{(α) +}$ . The derivation of the above equations has the assumption that the light is ncident only from the first layer whose amplitude spectrum $(a_{1}^{(α) +})$ is completely known. The computed amplitudes along with the W_l and V_l matrices will aid to determine the total electric and magnetic field components at each position in space using Eq. (3).

The proposed RSTWA technique has multiple advantages in comparison to other existing techniques for active-tunable materials. For example, the implementation of the numerical algorithm mainly utilizes Eq. (2), (5), (8), (11), (12), (15), and (16). The algorithm up to Eq. (15) computes the reflection and transmission matrices of interfaces of the structure in an isolated environment (independent of it position z_l). Hence, to minimize computational time and resources, one can compute the ( $r_{l}^{\pm}$ and $t_{l}^{\pm}$ ) matrices ‘once’ during the implementation and perform parametric studies of total reflection and transmission coefficients of diffraction orders from the full structure as a function of layers’ thickness using Eq. (16). In addition, the absence of space or time meshing makes the performance of the algorithm almost independent of the wavelength/time-period of the wave vs pattern’s physical/temporal dimensions. One can verify the convergence of the solution with respect to number of diffraction orders by simply checking the ability of exactly reproducing the shape of the space-time unit cell.

The primary limitation of the proposed RSTWA inherited from RCWA is that it is only applicable to planar layered systems with periodic modulation in both spatial and temporal dimensions. To compute light propagation through complex structured gratings such as sawtooth shape, cylindrical, or arbitrary shape, one may need to divide the cross section of the structure into a large number of rectangular layers which limit the accuracy and speed of computation. A full wave solver such as FDTD (described in Appendix II) is recommended in such cases. Furthermore, RSTWA solution yields the steady-state response of the system in time-domain and cannot capture transient phenomena.

Note that the above presented formulation is limited to 1D+time periodic layered systems. However, the method can be readily extended to 2D+time periodic layered systems by considering two-dimensional spatial diffraction orders and their corresponding frequency harmonics. For instance, the effect of a second spatial dimension can be incorporated by adding the wavevector set k_y_, _j = k_y_;0 + j2π/Λ_y to Eq. (2), where j is an integer number representing the order of spatial harmonic in the y-direction. Since s and p polarized waves cross-couple at oblique illumination of two-dimensional gratings, the null matrix terms in Eq. (5) will then be replaced by non-zero cross-coupling coefficients of 2D systems as described in [57]. Size of the matrices in this case grow to the order of $n_{\max}^{3}$ , with n_max being the maximum order of spatial and temporal harmonics, consuming enormous time and computational resources.

3. Numerical results and discussion

In this section, we present and discuss some numerical results obtained by our proposed RSTWA method. Our goal is to demonstrate the validity of the method, and exemplify its advantages in terms of superior performance and versatility, through two practical, yet complicated, problems of high current interest to the nanophotonics and optical materials community.

3.1. Validation

Here, we report the results of our investigations to verify the accuracy of the RSTWA method. The configuration that we have chosen for this purpose is a periodic arrangement of time-modulated dielectric bars with sinusoidal permittivity variation: ϵ(t) = ϵ_s [1 + δ sin(2π f_mt)], where ϵ_s is the static electric permittivity, δ is the modulation depth, and f_m is the modulation frequency (Ω = 2π f_m is defined in the formulation developed in Sec. 2). The structure is illuminated by the normal incidence of a plane wave with a frequency of f₀ = 3 THz. Both transverse-electric (TE, also known as s-polarization) and transverse-magnetic polarizations (TM, also known as p-polarization) have been considered. The structural parameters are given in the caption of Fig. 2.

Fig. 2 Comparison between the results obtained by the RSTWA, and a time-variant FDTD method (please see Appendix I), corresponding to a periodic arrangement of time-modulated dielectric bars with sinusoidal permittivity variation: ϵ(t) = ϵ_s[1+δ sin(2π f_mt)]. Parameters: f₀ = 3 THz (excitation frequency), D = 50 µm (width and height of the bars), Λ = 95 µm (periodicity), ϵ_s = 3.9, δ = 0.3, f_m = 0.125 f₀. (a) Magnitude of the transmission coefficient versus temporal harmonic index corresponding to the TM(p) polarization. The left pane represents the fundamental, and the right pane represents the ±1 spatial diffraction orders, respectively. (b), (c) Near-field plots of the magnitude of the x-component of the electric field, corresponding to the TM polarizations obtained by RSTWA and FDTD, respectively. Near-fields are plotted for −1, 0, +1 temporal harmonics (left, center, and right panes, respectively). (d)–(f) Same plot as (a)–(c) obtained for the TE(s) polarization.

Download Full Size | PDF

We have employed a full-wave in-house developed time-variant FDTD solver (please see appendix I for more information) to validate the accuracy of the RSTWA method. The simulation time of time-variant FDTD should be considered proportional to the period of the modulating signal, to capture the steady-state response of the time-modulated system. We have therefore considered a relatively high modulation frequency here (f_m = 0.125 f0) to shorten the FDTD simulations. As it will become evident shortly, this consideration, not only does not sacrifice the quality of validation, but actually elevates it, as it is usually more challenging for frequency domain methods, including RSTWA, to handle high modulation frequencies, and also richer physical phenomena start to emerge as the modulation frequency becomes comparable to the frequency of excitation.

We have obtained both far-field and near-field responses of the system by RSTWA and FDTD, and plotted the results in Fig. 2. Figures 2(a)–(c) show the results corresponding to the TM excitation and Figs. 2(d)–(f) correspond to the TE polarization. For more information regarding the plots in Fig. 2 please refer to the figure caption. Careful examination of Fig. 2 reveals excellent agreement between the RSTWA and FDTD results.

An interesting physical phenomenon, which occurs due to the high modulation frequency, is the exclusive presence of the temporal harmonics with harmonic indices larger than +1 (frequencies larger than f₀ + f_m) in the ±1 diffraction orders (spatial harmonics), and is evident in both far-field (right panels of Figs. 2 (a) and (d)) and also near-field plots (the interference pattern in the right pane of Fig. 2 (b)).

It is worth mentioning that, the simulation time required to obtain the results reported in Fig. 2 by the RSTWA method using |m| ⩽ 20 spatial and |n| ⩽ 11 temporal modes, on a workstation computer with Intel Xeon E5-2670 v3 processor, is about 6 s, while the FDTD simulation time is about 43 min. The improvement in simulation time for RSTWA, compared to FDTD, becomes even more drastic when smaller–more practical–modulation frequencies are considered, as in the applications that will be presented in the following, to the extent that simulating the time-modulated grating by FDTD becomes impractical.

3.2. Anomalous beam-steering with spatially phase-variant time-modulated dielectric metasurfaces

Addition of time-modulation to engineered material patterns (metasurfaces) increases dimensionality of the design space which can extend the degree of light manipulation. Our group has recently proposed a new design paradigm for wavefront engineering of light in time-modulated metasurfaces, based on modulation-induced phase shift of higher-order harmonics when a phase delay is introduced in the temporal modulation of elements [23, 32]. However, identification of structural optical design of the pattern is found to be a tedious task utilizing conventional techniques. Thanks to the great efficiency afforded by RSTWA method, we are able to conduct a comprehensive parametric study which provides a deep understanding of photonic transitions in time-modulated metasurfaces and their dependency to the geometrical parameters of the unit cell and associated optical modes.

Here, we design and model a time-modulated dielectric metasurface within the framework of RSTWA method at the near-infrared wavelength of λ₀ =1.55 µm and utilize modulation-induced phase shift to demonstrate anomalous bending of higher-order harmonics. For this purpose, we consider an array of silicon nanobars whose refractive index is sinusoidally modulated in time with a weak modulation depth as n(t) = n₀ [1 + δ sin(ω_mt)] in which n₀ = 3.475, δ = 8.6 × 10⁻³ and f_m = 10 GHz. Such a weak perturbation of silicon’s refractive index can be achieved over large length scales (up to mm) via current-driven double-carrier injection into the intrinsic region of a p-i-n configuration [58,59]. The carrier-induced changes of refractive index via plasma-Drude dispersion has been previously exploited for realization of silicon modulators in the near-infrared regime with modulation frequencies up to several GHz [10]. The schematic of time-modulated metasurface is illustrated in Fig. 3(a). An array of silicon nanobars with a periodicity of Λ = 650 nm is located on top of a SiO₂ spacer layer backed by a silver mirror and the structure is illuminated with a normally incident TE-polarized light. The fabrication of the proposed structure can be carried out in a complementary metal-oxide semiconductor (CMOS)-compatible process by using e-beam lithography to develop the photoresist masks for implantation of the dopants and etching the silicon nanobars [11]. A radio-frequency (RF) biasing signal can be applied between highly p-doped and n-doped regions to modulate the refractive index of the intrinsic region via double-carrier injection. Incorporating RF phase shifters in the biasing network of the elements, each element can be addressed independently to realize a spatially phase-variant temporal profile for the refractive index as shown in Fig. 3(b). The periodic structure is analyzed within the framework of RSTWA method by considering a 2D model (Fig. 3(c)) while using |m| ⩽ 20 spatial and |n| ⩽ 12 temporal harmonics. The refractive index of silica is chosen as 1.444 and the complex permittivity of silver is given by the Drude model. The thickness of silver back-mirror is set as 150 nm and the thickness of dielectric spacer is chosen as 268 nm corresponding to a quarter effective wavelength which leads to formation of an asymmetric Fabry-Pérot cavity and satisfaction of critical coupling condition.

Fig. 3 Schematic depiction of a spatially phase-variant time-modulated dielectric metasurface. (a) An array of silicon nanobars with p-i-n configuration located on top of a layer of SiO₂ backed by a silver mirror are modulated using an RF biasing network incorporating phase shifters. (b) The spatio-temporal profile of the refractive index with spatially variant phase shift, resembling a space-time diffraction grating. The dashed lines in (a) and (b) denote the line along which spatio-temporal refractive index profile is plotted. (c) Cross-sectional view of the unit cell.

Download Full Size | PDF

The frequency conversion efficiency of the time-modulated metasurface is defined as $| \frac{r (f_{0} + f_{m})}{r (f_{0})} |$ and the calculated result using RSTWA is plotted in Fig. 4 as a function of width and height of silicon nanobar. Different branches exhibiting resonant enhancement in the frequency conversion efficiency are identified which correspond to the resonant modes supported by the metasurface. We characterize the resonant optical modes of different branches by plotting the distribution of electric field for different heights and widths of the nanobar. The nearfield results corresponding to the points (b-e) denoted on the colormap of Fig. 4(a) are brought in Figs. 4(b,c-1,d,e-1). The nearfield distribution in Fig. 4(b) exhibits two antinodes of opposite signs at the top and bottom of nanobar which implies that the net electric field inside the silicon nanobar is zero in x-y plane and the displacement current is rotating in y-z plane. This is characteristic of magnetic dipole (MD) resonance which is supported by high-index nanoresonators [60]. The electric field in Fig. 4(c-1) is mainly confined at the center of nanobar while two additional anti-nodes can be observed at the top and bottom of the nanobar, forming a standing-wave pattern along the height. This distribution arises from hybridization of an electric dipole (ED) mode excited in the transverse direction with a Fabry-Pérot resonance excited in the longitudinal direction [61]. Both electric field distributions in Figs. 4(d-1) and (e) exhibit a transverse wavefront corresponding to guided modes inside the grating layer [62], while the distribution in Fig. 4(e) shows additional antinodes in the longitudinal direction which is occurred due to hybridization of the guided mode with Fabry-Pérot resonance at longer heights of the silicon nanobar. As such, it can be concluded that the horizontally-oriented branches in Fig. 4 correspond to hybrid Mie-Fabry-Pérot resonances and the vertically-oriented branches represent hybrid guided-Fabry-Pérot resonances. Unlike Mie-Fabry-Pérot resonant modes, the guided mode resonances result into frequency conversion efficiencies much greater than unity in the account of infinite interaction length of guided light with the modulation. Figures 4(c-2) and (d-2) compare the spectrum of generated frequency harmonics corresponding to the hybrid ED-Fabry-Pérot resonance (point c) and guided mode resonance (point d), respectively. In the former case, a weak photonic transition is observed from the fundamental frequency harmonic to the first-order frequency harmonics and all other higher order harmonics are weakly excited. This is while in the guided mode resonance, a strong photonic transition occurs from the fundamental frequency harmonic into all higher-order frequency harmonics and the guided light undergoes an almost complete frequency conversion. Such strong photonic transitions can be used for achieving stronger non-reciprocal responses from time-modulated metasurfaces with smaller modulation frequencies [63]. The amplitude of generated frequency harmonics shows an exponentially decaying profile vs the harmonic index. It should be noted that the amplitude of up-and down-modulated frequency harmonics are almost equal due to the small frequency of modulation compared to the optical frequency which results into negligible dispersion effects. It should be mentioned that the sub-nm spectral resolution of current near-infrared spectrometers (e.g. FHR spectrometer from Horiba [64]) allows for accurate resolving of generated frequency harmonics separated at < 10 GHz intervals in scattering experiments. Moreover, such narrow spectral spacings have been commonly used in dense wavelength division multiplexing communication systems [65].

Fig. 4 (a) Frequency conversion efficiency of the dielectric time-modulated metasurface of Fig.3 in the logarithmic scale as a function of height and width of silicon nanobar. The dashed lines denote the branches of resonant modes. (b), (c-1), (d-1) and (e) depict the electric field distribution corresponding to the widths and heights of silicon nanobar at the resonant points denoted in the colormap of (a). The distributions exhibit characteristics of MD, ED-Fabry-Pérot, guided and guided-Fabry-Pérot resonant modes, respectively. (c-2) and (d-2) compare the spectrum of generated frequency harmonics at the resonant points of (c) and (d) corresponding to ED-Fabry-Pérot and guided mode resonances.

Download Full Size | PDF

In order to engineer the emerging wavefronts of higher-order harmonics, we introduce a spatially varying phase shift to the temporal modulation of elements as n(x, t) = n₀(1 + δ sin(ω_mt + α(x))) which can be implemented using RF phase shifters in the biasing network as shown in Fig. 3(a). A progressive phase shift of Δα between the adjacent elements translate into $α (x) = \frac{Δ α}{Λ} x$ and creates a linear phase gradient profile of $ϕ_{n} (x) = n \frac{Δ α}{Λ} x$ across the metasurface for the scattered fields of n-th frequency harmonic through modulation-induced phase shift [32]. Such a linear phase gradient leads to anomalous bending of higher-order frequency harmonics into the angles of $θ_{n} = \sin^{- 1} (- \frac{n Δ α}{k_{n} Λ})$ . Due to the non-resonant nature of modulation-induced phase shift, this anomalous bending effect can be achieved for any choice of height and width for subwavelength Si nanobars, while the deflection efficiency is enhanced at the resonant branches shown in Fig. 4(a). Despite the strong photonic transitions due to guided modes enabling total frequency conversion, their efficiency in beam-steering is typically not as high compared to the photonic transitions associated with Mie-Fabry-Pérot resonances. This stems from the fact that in case of guided modes, the radiation is provided by the space-time modal transitions from guided modes to leaky modes possessing complex-valued spatial frequencies which are not perfectly matched with the wavevectors of propagating plane waves. This imperfect phase matching depresses the frequency conversion efficiency in beam-steering. Here, we set the height and width of the nanobars as h = 511 nm and w = 332 nm, corresponding to the hybrid ED-Fabry-Pérot resonance. The progressive modulation phase delay is set as Δα = −k₀Λ sin(20^°) in order to route the first order frequency harmonics into the angles of θ_±1 ≈±20^° and a supercell of 7 elements is chosen to perfectly discretize the required phase distribution for beam-steering. It should be remarked that in the account of modulation frequency being much smaller than the optical frequency, k₊_n ≈ k₋_n ≈ k₀ which leads to a mirror symmetry between up-and down-modulated frequency harmonics (θ₊_n ≈−θ₋_n). Figure 5 demonstrates the wavefronts of electric field corresponding to the fundamental and first-order frequency harmonics. While the fundamental frequency harmonic is governed by the specular reflection due to zero modulation-induced phase shift, first-order frequency harmonics are anomalously bent into opposite directions as predicted by the theory.

Fig. 5 The wavefronts of the electric field corresponding to (a) down-modulated, (b) fundamental and (c) up-modulated frequency harmonics generated by the spatially phase-variant dielectric time-modulated metasurface. The fundamental frequency harmonic is specularly reflected while the first-order frequency harmonics are anomalously bent toward angles of θ_±1 ≈±20^°.

Download Full Size | PDF

3.3. Dynamic beam generation with spatially frequency-variant time-modulated plasmonic metasurfaces

Another exciting opportunity offered by RSTWA is the ability to model and engineer steady-state spatio-temporal dynamics of light through manipulating the interference pattern of space-time harmonics. Generation of dynamic beams by a spatio-temporally modulated metamaterial lens based on a transformation optics approach has been demonstrated in Ref [23]. Here, in order to achieve ultrafast spatio-temporal beam-scanning via ultrathin time-modulated metasurfaces, we exploit the concept of frequency arrays which has been recently proposed by Shaltout et al. [66, 67]. In contrast to phased arrays which operate based on the interference of monochromatic scattered fields from the constituent elements with a progressive phase delay, a frequency array relies on the coherent interference of scattered fields with a progressive shift in the temporal frequency which can be translated to a spatio-temporal phase delay [66]. A frequency array can be realized through frequency comb excitation of a passive metasurface with the ability to focus different wavelengths into an array of diffraction-limited focal spots at a certain focal plane such that the resulting focal spots will act as secondary sources for the frequency array [67]. However, such an implementation requires very large numerical apertures for the metasurface which can increase the footprint of the device and also has a performance limit in the beam directivity and scanning speed due to the limitations in focusing resolution. Here, we demonstrate dynamic beam generation using ultrathin time-modulated based on the notion of frequency arrays which can enable continuous beam-scanning with 180^° angle-of-view, tunable speed and ultracompact footprint. To this end, we design and model a time-modulated plasmonic metasurface based on graphene micro-ribbons at the far-infrared wavelength of λ₀ = 25 µm within the framework of RSTWA method.

Graphene is modeled with a volumetric thickness-dependent permittivity given as $ϵ (ω) = \frac{i σ_{s} (ω)}{ω ϵ_{0} t_{g}}$ , where t_g is the thickness and σ_s(ω) is the surface conductivity of graphene sheet. The thickness of graphene is considered to be t_g = 0.5 nm in RSTWA simulations which is found to be sufficient for convergence of the results in the far-infrared wavelength regime. It should be remarked that RSTWA offers a great advantage in simulation of such ultrathin layers, since no volume discretization is necessary. The conductivity response of graphene sheet is dominated by the intra-band transitions in this wavelength range which can be written as $σ_{s} (ω) = \frac{e^{2} E_{f}}{π ℏ^{2}} \frac{i}{ω + i τ^{- 1}}$ for moderate levels of Fermi energy level at the room temperature [68], where e is electron charge, E_f is Fermi energy level, ħ is reduced Planck’s constant and τ is the energy relaxation time. Here, we have considered τ = 0.5 ps corresponding to a high-quality graphene monolayer, and a sinusoidally-modulated temporal profile for Fermi energy level as E_f = E_f0 [1 + δ sin(ω_mt)]. The modulation of Fermi energy level of graphene has been experimentally demonstrated in the range of 0 – 0.6 eV using electrical gate-biasing [7], and modulation frequencies up to 30 GHz have been achieved in graphene electro-optical modulators [14]. Accordingly, we consider E_f₀ = 0.3 eV and δ = 0.8 such that maximal modulation depth is achieved with moderate Fermi energy levels.

The time-modulated metasurface is consisted of an array of graphene micro-ribbons with a deeply subwavelength periodicity of Λ = λ₀/8 = 3.125 µm located on top of a SiO2 layer backed by a silver mirror with a thickness of 300 nm, under illumination of a normally incident TM-polarized plane wave as schematically depicted in Fig. 6(a). Such micro-ribbon structures have been fabricated by growing graphene using chemical vapor deposition and adopting standard optical lithography followed by oxygen plasma etching [69]. The surface conductivity of each micro-ribbon is modulated in a parallel capacitor configuration by applying an external RF biasing signal between micro-ribbons and the back-mirror. Using different modulation frequencies for the elements, a spatially frequency-variant temporal profile for the surface conductivity can be achieved as shown in Fig. 6(b). The periodic structure is considered as a 2D model (Fig. 6) within the framework of RSTWA method by truncating the spatial harmonics to ±85 and temporal harmonics to ±10. It should be noted that similar to the conventional RCWA method, RSTWA requires larger number of spatial harmonics for convergence to the exact results in case of TM polarization and gratings with large effective depths [47].

Fig. 6 Schematic depiction of a spatially frequency-variant time-modulated plasmonic metasurface. (a) An array of graphene micro-ribbons located on top of a layer of SiO2 backed by a silver mirror are modulated with RF sources of different frequencies. (b) A spatially frequency-variant profile for surface conductivity. The dashed lines in (a) and (b) denote the line along which the surface conductivity profile is plotted. (c) Cross-sectional view of the unit cell.

Download Full Size | PDF

Figure 7(a) depicts the calculated frequency conversion efficiency of the metasurface as a function of width of the micro-ribbons and height of dielectric spacer layer for a modulation frequency of f_m = 1 GHz. The frequency conversion efficiency shows a periodic oscillation by changing the height of dielectric spacer associated with the Fabry-Pérot interference patterns while a resonant enhancement in the efficiency can be observed when the width of graphene ribbons is around w = 1 µm which corresponds to the plasmonic resonance of micro-ribbons [6]. When the thickness of dielectric spacer is equal to an odd integer number of quarter effective wavelength inside the dielectric layer, the frequency conversion efficiency becomes maximal in the plasmonic resonance regime due to the critical coupling condition. In the following, we choose the width of the micro-ribbons as w = 1µm and the hight of dielectric spacer as h = 4.34 µm, corresponding to the first resonant peak in Fig. 7(a). The distribution of electric field is shown at Fig. 7(b) which exhibits an extreme confinement to the micro-ribbons. Figure 7(c) depicts the spectrum of generated frequency harmonics of the time-modulated metasurface showing an exponential decay profile vs harmonic index.

Fig. 7 (a) Frequency conversion efficiency of the time-modulated graphene metasurface as a function of width of the graphene micro-ribbons and height of the dielectric spacer layer. The dashed lines denote the plasmonic resonant regime. (b) Nearfield distribution of the electric field at the plasmonic resonance of the micro-ribbon, when w = 1 µm and h = 4.34 µm. (c) Generated frequency harmonic spectrum by the time-modulated metasurface at the plasmonic resonant regime.

Download Full Size | PDF

Due to the significant confinement of the field in the vicinity of micro-ribbons, the interaction of neighbor elements is negligible and each micro-ribbon acts as a resonant frequency harmonic generator. The frequency of higher-order harmonics can be tuned by changing the modulation frequency without affecting their phase and amplitude, as long as the modulation frequency is much smaller than the optical frequency (f_m ≪ f₀) such that both the structural and material dispersion effects are negligible. Accordingly by applying a progressive shift in the modulation frequency of the micro-ribbons (Δf_m), a multi-frequency array is achieved in which the coherent interference of the frequency harmonics leads to generation of dynamic beams in the space-time domain. In particular, three dominant beams are expected resulting from the interference of fundamental, up-and down-modulated frequency harmonics of the first order. The generated beam from the interference of fundamental frequency harmonics is stationary due to invariance of the temporal frequency and its direction is determined by the specular reflection angle. The other two beams are dynamic whose steering angles are time-dependent and given by [66]:

\sin (θ_{\pm 1} (t)) = \pm (\frac{2 π Δ f_{m}}{k_{0} Λ_{\pm 1}}) t

for |t| < T_±1/2 with

T_{\pm 1} = \frac{1}{| Δ f |}

being the temporal periodicity of the beams. These two beams periodically scan the space-time domain in the opposite directions with a 180^° angle-of-view. For subwavelength periodicities of Λ_±1 <λ₀/2, each beam can experience a dead time-zone region (corresponding to | sin(θ_±(t))| > 1) where the beam is not radiated into free-space. On the other hand, if the periodicity exceeds λ₀/2, each beam can be spitted into multiple directions (corresponding to the case where sin(θ_±(t)) is not unique). In order to avoid occurrence of dead time-zones and beam-splitting, we consider a supercell of N = 44 elements in which each 4 neighbor elements are modulated with a progressive shift of Δf_m = 1 GHz in the modulation frequency such that the periodicity of resulting frequency arrays is Λ_±1 = λ₀/2 and modulation frequency spans from f_m_,1 = 1 GHz to f_m_, _N = 11 GHz. We excite the metasurface with a Gaussian beam having a waist of W₀ = 2.25λ₀ to clearly resolve the generated beams in space. For this purpose, we use a plane-wave expansion of the Gaussian beam into 121 plane waves having different transverse wavevectors [56] and truncate spatial harmonics to ±100 and temporal harmonic to ±11 in RSTWA calculations. Figure 8 demonstrates the spatio-temporal beam-scanning performance of the array by plotting the wavefront of time-domain electric field at different instants within a time frame of T = 1 ns. Three dominant beams are clearly observed with two dynamic beams rotating in opposite directions as the time progresses while covering 180^° angle-of-view. Additional dynamic beams can also be seen in the nearfield plots which are generated due to coherent interference of harmonic orders n > 1. These beams have a temporal periodicity of

T_{n} = \frac{1}{| n Δ f_{m} |}

and are typically much weaker compared to the first-order dynamic beams due to exponential decay in the amplitude with increasing the harmonic index. It should be mentioned that the current state of THz detectors allows for spectral resolutions of tens of MHz and measurement of ultrafast dynamic spatiotemporal beams by femtosecond real-time delays in time-domain spectroscopy [70].

Fig. 8 Wavefront of the time-domain electric field at different time instants within a time frame of T = 1 ns corresponding to the spatially frequency-variant time-modulated graphene metasurface. Three dominant dynamic beams can be seen as a result of the coherent interference of fundamental, up-and down-modulated frequency harmonics. The beam resulting from the interference of fundamental harmonic is stationary while the other two beams are dynamic and scan the space in time with a 180^° angle-of-view in opposite directions.

Download Full Size | PDF

4. Conclusion

In conclusion, we have presented a space-time coupled method, RSTWA, to compute light propagation through patterned layers of time-modulated materials. The validity of the method is verified by comparing the simulation results for a time-modulated dielectric grating with those of a time-variant FDTD technique. Excellent agreement has been achieved, while affording a great efficiency by the RSTWA method. The method can efficiently characterize periodic gratings with spatio-temporal modulation, having deeply subwavelength features in the spatial domain. The developed algorithm typically computes the light propagation characteristics through rectangular structure in few seconds and hence parametric studies of large number of designs can be performed in few minutes. Here, we employed the technique to design and model dielectric and plasmonic time-modulated metasurfaces implemented by electro-optical modulation of silicon and graphene. Spatial variation of the modulation phase and frequency are utilized to anomalously steer the higher-order frequency harmonics, and generate dynamic beams in space-time domain, respectively. The presented formulation can easily be extended to two-dimensional time-modulated gratings, and therefore be employed as a valuable tool for design and characterization of dynamic metasurfaces. The compact, fast, and accurate nature of the algorithm makes it readily available to be integrated with numerical optimization techniques for tailored pattern designs.

Appendix I Derivation of eigenvalue equation for the field profiles

In this appendix, we discuss derivation of Eq. (7) corresponding to the eigenvalue equation in a spatiotemporally modulated grating layer in further details. For this purpose, we start with the Maxwell’s curl equations for the electromagnetic fields $\vec{E} = {E_{x}, E_{y}, E_{z}}$ and $\vec{H} = {H_{x}, H_{y}, H_{z}}$ :

\begin{array}{l} \nabla \times \vec{E} = - μ_{0} μ_{r} \frac{\partial \vec{H}}{\partial t}, \\ \nabla \times \vec{H} = ϵ_{0} \frac{\partial}{\partial t} (ϵ (x, t) \vec{E}) \end{array}

which result into the following set of first order partial differential equations for the field components:

\frac{\partial E_{y}}{\partial z} = μ_{0} μ_{r} \frac{\partial H_{x}}{\partial t}

\frac{\partial H_{y}}{\partial z} = - ϵ_{0} \frac{\partial}{\partial t} (ϵ (x, t) E_{x})

\frac{\partial E_{y}}{\partial x} = - μ_{0} μ_{r} \frac{\partial H_{z}}{\partial t}

\frac{\partial H_{y}}{\partial x} = ϵ_{0} \frac{\partial}{\partial t} (ϵ (x, t) E_{z})

\frac{\partial H_{x}}{\partial z} - \frac{\partial H_{z}}{\partial x} = ϵ_{0} \frac{\partial}{\partial t} (ϵ (x, t) E_{y})

\frac{\partial E_{z}}{\partial x} - \frac{\partial E_{x}}{\partial z} = μ_{0} μ_{r} \frac{\partial H_{y}}{\partial t}

As shown in Eq. (3), we express the field components as summations over space-time diffraction orders as ${E, H}_{{x, y, z}} (\vec{r}) = \sum_{α} e^{i k_{x, α} x} e^{- i ω_{α} t} {ℰ, ℋ}_{{x, y, z}}^{(α)} (z)$ where the summation over α represents a double indices (m, n). In the above set of equations, the spatial derivative along x– axis of all field components can be evaluated as:

\frac{\partial}{\partial x} {E, H}_{{x, y, z}} (\vec{r}) = i \sum_{α} k_{x, α} e^{i k_{x, α} x} e^{- i ω_{α} t} {ℰ, ℋ}_{{x, y, z}}^{(α)} (z)

and the time derivative all field components can be evaluated as:

\frac{\partial}{\partial t} {E, H}_{{x, y, z}} (\vec{r}) = - i \sum_{α} ω_{α} e^{i k_{x, α} x} e^{- i ω_{α} t} {ℰ, ℋ}_{{x, y, z}}^{(α)} (z)

Utilizing Fourier convolution theorem the products ϵ(x, t){E_x, E_y, E_z } can be transformed into summation as follows:

\begin{matrix} ϵ (x, t) E_{{x, y, z}} (\vec{r}) = \sum_{α} ϵ_{α} e^{(- i q_{x, α} x + i q_{t, α} t)} \sum_{α} e^{i k_{x, α} x} e^{- i ω_{α} t} ℰ_{{x, y, z}}^{(α)} (z) \\ = \sum_{α, β} ϵ^{α, β} e^{i k_{x, β} x} e^{- i ω_{β} t} ℰ_{{x, y, z}}^{(β)} (z) \end{matrix}

whose time derivatives are given by:

\frac{\partial}{\partial t} ϵ (x, t) E_{{x, y, z}} (\vec{r}) = - i \sum_{α, β} ω_{β} ϵ^{α, β} e^{i k_{x, β} x} e^{- i ω_{β} t} ℰ_{{x, y, z}}^{(β)} (z)

Now the aim is to obtain coupled differential equations of tangential field components (E_y and H_x) for s-polarization. Plugging Eq. (28) into (20) will directly result into:

\frac{\partial}{\partial z} ℋ_{y}^{α} (z) = i ϵ_{0} \sum_{α, β} ω_{β} ϵ^{α, β} ℰ_{x}^{(β)} (z)

which can be rewritten in the matrix form as:

{\dot{ℋ}}_{y} = i ϵ_{0} W Ξ ℰ_{x} .

Eq. (24) leads to:

\frac{\partial}{\partial z} ℰ_{x}^{α} (z) = i μ_{0} μ_{r} ω_{α} ℋ_{y}^{(α)} (z) + i k_{x, α} ℰ_{z}^{α}

Using Eq. (22), we can obtain $k_{x, α} ℋ_{y}^{α} = - ϵ_{0} \sum_{α, β} ω_{β} ϵ^{α, β} ℰ_{z}^{β}$ which can be used in (31) to obtain the following matrix equation:

{\dot{ℰ}}_{x} = i {μ_{0} μ_{r} W - \frac{1}{ϵ_{0}} K_{x} Ξ^{- 1} W^{- 1} K_{x}} ℋ_{y}

A similar process can be carried for p–polarization field components (E_x and H_y). By substituting Eq. (25) into (19) yields:

\frac{\partial}{\partial z} ℰ_{y}^{α} (z) = - i μ_{0} μ_{r} ω_{α} ℋ_{x}^{(α)} (z)

leading to the following matrix equation:

{\dot{ℰ}}_{x} = - i μ_{0} μ_{r} W ℋ_{x}

From Eq. (21) it can be concluded that $ℋ_{z, α} = \frac{1}{μ_{0} μ_{r}} ω_{α}^{- 1} k_{x, α} ℰ_{y}^{α}$ whose substitution into (23) leads to:

\frac{\partial}{\partial z} ℋ_{x}^{α} (z) = i \sum_{α, β} {\frac{1}{μ_{0} μ_{r}} k_{x, α} ω_{α}^{- 1} k_{x, α} - ϵ_{0} ω_{β} ϵ^{α, β}} ℰ_{y}^{(β)} (z)

which can be expressed in the matrix form as:

{\dot{ℋ}}_{x} = i {\frac{1}{μ_{0} μ_{r}} K_{x} W^{- 1} K_{x} - ϵ_{0} W Ξ} ℰ_{y}

Eqs. (30), (32), (34) and (36) can be recast as the eigenvalue equation given in (7).

Appendix II Time-variant FDTD

Application of the finite-difference time-domain (FDTD) method to the analysis of time-varying materials has been previously reported in the literature [71]. Recently, time-variant FDTD method has been employed to analyze the space-time-modulated metasurface [72, 73]. The permittivity variation of a space-time-modulated Huygens’ metasurface, with finite lateral dimensions, was first mapped to the equivalent Lorentzian susceptibilities of a zero-thickness sheet, using the generalized sheet transition conditions (GSTCs). FDTD method was then employed to solve the electromagnetic scattering from the equivalent transition sheet.

In this paper, we employ FDTD to validate the accuracy of our proposed RSTWA method. Our approach is, however, different from the one proposed in Ref [73], in that it takes the temporal and spatial variations of the material parameters directly into the FDTD formulation. It is therefore simpler and needs less preprocessing, without any need for a transition condition, and can be easily extended to handle more complex material models of interest to the Nanophotonics community.

In this section, we present the formulation of our developed time-variant FDTD method. Starting from the Ampere’s law,

\frac{\partial \vec{D}}{\partial t} = \nabla \times \vec{H} - \vec{J}

where

\vec{D}

,

\vec{E}

, and

\vec{J}

are electric flux density, electric field, and electric current density, respectively, constitutive relation for the electric flux density of linear, time-variant materials can be written as,

\vec{D} = ϵ (\vec{r}, t) \vec{E} = ϵ_{0} ϵ_{r} (\vec{r}, t) \vec{E}

where

ϵ (\vec{r}, t)

represents the space-time-varying electric permittivity, and ϵ₀ is the electric permittivity of vacuum. Also, ϵ_r represents the relative electric permittivity. Electric current density in Eq. 37 can also be expressed as,

\vec{J} = σ_{e} (\vec{r}, t) \vec{E} + {\vec{J}}_{imp}

where

σ_{e} (\vec{r}, t)

is the space-time-varying electric conductivity, and

{\vec{J}}_{imp}

is the imposed electric current density. Substituting Eqs. 38, 39 in 37, we obtain

\frac{\partial}{\partial t} (ϵ_{r} (\vec{r}, t) \vec{E}) = \frac{1}{ϵ_{0}} (\nabla \times \vec{H} - σ_{e} (\vec{r}, t) \vec{E} - {\vec{J}}_{imp}) .

Equation 40 can be discretized, using the Yee algorithm, along with the temporal and spatial central differences [51]. The resulting update-equations for components of the the electric field can be obtained as,

{\vec{E}}^{n + 1} = C_{1} {\vec{E}}^{n} + C_{2} (\nabla \times {\vec{H}}^{n + 1 / 2} - {\vec{J}}_{imp}^{n + 1 / 2})

where

C_{1} = \frac{2 ϵ_{0} ϵ_{r}^{n} - σ_{e}^{n + 1 / 2} Δ t}{2 ϵ_{0} ϵ_{r}^{n + 1} + σ_{e}^{n + 1 / 2} Δ t}; C_{2} = \frac{2 Δ t}{2 ϵ_{0} ϵ_{r}^{n + 1} + σ_{e}^{n + 1 / 2} Δ t} .

where superscripts represent the timing of evaluating the material parameters, and calculating the field components [51].

The time-variant FDTD algorithm, based on the formulation outlined above, is similar to the conventional FDTD algorithm, except that the electric field update-coefficients are not constant throughout the simulation. Instead, material parameters and the corresponding update coefficients need to be updated inside the FDTD time-marching loop, prior to updating the electric field components. It is noteworthy that, although the simulation time for the conventional time-invariant FDTD is proportional to the time period of the excitation signal, it must be considered proportional to the period of material parameter modulation for time-variant FDTD. After the steady-state has been reached, sampled field data in the time domain should be gathered at least in one period of modulation for time to frequency domain and subsequently near-field to far-field transformation.

Funding

U.S. Air Force Office of Scientific Research (AFOSR) (FA9550-18-1-0354).

References

1. A. Tittl, A.-K. U. Michel, M. Schäferling, X. Yin, B. Gholipour, L. Cui, M. Wuttig, T. Taubner, F. Neubrech, and H. Giessen, “A switchable mid-infrared plasmonic perfect absorber with multispectral thermal imaging capability,” Adv. Mater. 27, 4597–4603 (2015). [CrossRef] [PubMed]

2. Z. Zhu, P. G. Evans, R. F. Haglund Jr, and J. G. Valentine, “Dynamically reconfigurable metadevice employing nanostructured phase-change materials,” Nano letters 17, 4881–4885 (2017). [CrossRef] [PubMed]

3. A. Forouzmand and H. Mosallaei, “Dynamic beam control via mie-resonance based phase-change metasurface: a theoretical investigation,” Opt. express 26, 17948–17963 (2018). [CrossRef] [PubMed]

4. Y.-W. Huang, H. W. H. Lee, R. Sokhoyan, R. A. Pala, K. Thyagarajan, S. Han, D. P. Tsai, and H. A. Atwater, “Gate-tunable conducting oxide metasurfaces,” Nano letters 16, 5319–5325 (2016). [CrossRef] [PubMed]

5. J. Park, J.-H. Kang, S. J. Kim, X. Liu, and M. L. Brongersma, “Dynamic reflection phase and polarization control in metasurfaces,” Nano letters 17, 407–413 (2016). [CrossRef] [PubMed]

6. Z. Li, K. Yao, F. Xia, S. Shen, J. Tian, and Y. Liu, “Graphene plasmonic metasurfaces to steer infrared light,” Sci. reports 5, 12423 (2015). [CrossRef]

7. M. C. Sherrott, P. W. Hon, K. T. Fountaine, J. C. Garcia, S. M. Ponti, V. W. Brar, L. A. Sweatlock, and H. A. Atwater, “Experimental demonstration of> 230 phase modulation in gate-tunable graphene–gold reconfigurable mid-infrared metasurfaces,” Nano letters 17, 3027–3034 (2017). [CrossRef]

8. A. Forouzmand, M. M. Salary, S. Inampudi, and H. Mosallaei, “A tunable multigate indium-tin-oxide-assisted all-dielectric metasurface,” Adv. Opt. Mater. 6, 1701275 (2018). [CrossRef]

9. Q. Xu, B. Schmidt, S. Pradhan, and M. Lipson, “Micrometre-scale silicon electro-optic modulator,” nature 435, 325 (2005). [CrossRef] [PubMed]

10. G. T. Reed, G. Mashanovich, F. Y. Gardes, and D. Thomson, “Silicon optical modulators,” Nat. photonics 4, 518 (2010). [CrossRef]

11. H. Lira, Z. Yu, S. Fan, and M. Lipson, “Electrically driven nonreciprocity induced by interband photonic transition on a silicon chip,” Phys. review letters 109, 033901 (2012). [CrossRef]

12. J. A. Dionne, K. Diest, L. A. Sweatlock, and H. A. Atwater, “Plasmostor: a metal-oxide-si field effect plasmonic modulator,” Nano Lett. 9, 897–902 (2009). [CrossRef] [PubMed]

13. M. Liu, X. Yin, and X. Zhang, “Double-layer graphene optical modulator,” Nano letters 12, 1482–1485 (2012). [CrossRef] [PubMed]

14. C. T. Phare, Y.-H. D. Lee, J. Cardenas, and M. Lipson, “Graphene electro-optic modulator with 30 ghz bandwidth,” Nat. Photonics 9, 511 (2015). [CrossRef]

15. Z. Yu and S. Fan, “Complete optical isolation created by indirect interband photonic transitions,” Nat. photonics 3, 91 (2009). [CrossRef]

16. S. Longhi, “Stopping and time reversal of light in dynamic photonic structures via bloch oscillations,” Phys. Rev. E 75, 026606 (2007). [CrossRef]

17. M. S. Kang, A. Brenn, and P. S. J. Russell, “All-optical control of gigahertz acoustic resonances by forward stimulated interpolarization scattering in a photonic crystal fiber,” Phys. review letters 105, 153901 (2010). [CrossRef]

18. Y. Sivan and J. B. Pendry, “Time reversal in dynamically tuned zero-gap periodic systems,” Phys. review letters 106, 193902 (2011). [CrossRef]

19. B. Dana, L. Lobachinsky, and A. Bahabad, “Spatiotemporal coupled-mode theory in dispersive media under a dynamic modulation,” Opt. Commun. 324, 165–167 (2014). [CrossRef]

20. Y. Hadad, J. C. Soric, and A. Alu, “Breaking temporal symmetries for emission and absorption,”Proc. Natl. Acad. Sci. 13201517363 (2016).

21. Y. Shi, W. Shin, and S. Fan, “Multi-frequency finite-difference frequency-domain algorithm for active nanophotonic device simulations,” Optica 3, 1256–1259 (2016). [CrossRef]

22. Y. Shi, A. Cerjan, and S. Fan, “Invited article: Acousto-optic finite-difference frequency-domain algorithm for first-principles simulations of on-chip acousto-optic devices,” APL Photonics 2, 020801 (2017). [CrossRef]

23. M. M. Salary, S. Jafar-Zanjani, and H. Mosallaei, “Time-varying metamaterials based on graphene-wrapped microwires: Modeling and potential applications,” Phys. Rev. B 97, 115421 (2018). [CrossRef]

24. A. Shaltout, A. Kildishev, and V. Shalaev, “Time-varying metasurfaces and lorentz non-reciprocity,” Opt. Mater. Express 5, 2459–2467 (2015). [CrossRef]

25. Y. Hadad, D. Sounas, and A. Alu, “Space-time gradient metasurfaces,” Phys. Rev. B 92, 100304 (2015). [CrossRef]

26. Y. Shi, S. Han, and S. Fan, “Optical circulation and isolation based on indirect photonic transitions of guided resonance modes,” ACS Photonics 4, 1639–1645 (2017). [CrossRef]

27. D. L. Sounas, C. Caloz, and A. Alu, “Giant non-reciprocity at the subwavelength scale using angular momentum-biased metamaterials,” Nat. communications 4, 2407 (2013). [CrossRef]

28. D. Correas-Serrano, J. Gomez-Diaz, D. Sounas, Y. Hadad, A. Alvarez-Melcon, and A. Alù, “Nonreciprocal graphene devices and antennas based on spatiotemporal modulation,” IEEE Antennas Wirel. Propag. Lett. 15, 1529–1532 (2016). [CrossRef]

29. S. Taravati and C. Caloz, “Mixer-duplexer-antenna leaky-wave system based on periodic space-time modulation,” IEEE Transactions on Antennas Propag. 65, 442–452 (2017). [CrossRef]

30. S. Taravati, N. Chamanara, and C. Caloz, “Nonreciprocal electromagnetic scattering from a periodically space-time modulated slab and application to a quasisonic isolator,” Phys. Rev. B 96, 165144 (2017). [CrossRef]

31. D. L. Sounas and A. Alù, “Non-reciprocal photonics based on time modulation,” Nat. Photonics 11, 774 (2017). [CrossRef]

32. M. M. Salary, S. Jafar-Zanjani, and H. Mosallaei, “Electrically tunable harmonics in time-modulated metasurfaces for wavefront engineering,” New J. Phys. :1801.10575 (2018).

33. M. Liu, D. A. Powell, Y. Zarate, and I. V. Shadrivov, “HuygensâĂŹ metadevices for parametric waves,” Phys. Rev. X 8, 031077 (2018).

34. N. Chamanara, Y. Vahabzadeh, and C. Caloz, “Simultaneous control of the spatial and temporal spectra of light with space-time varying metasurfaces,” arXiv preprint arXiv :1808.03385 (2018).

35. S. Taravati and A. A. Kishk, “Dynamic modulation yields one-way beam splitting,” arXiv preprint arXiv :1809.00347 (2018).

36. J. R. Zurita-Sánchez, P. Halevi, and J. C. Cervantes-Gonzalez, “Reflection and transmission of a wave incident on aslab with a time-periodic dielectric function Ït (t),” Phys. Rev. A 79, 053821 (2009). [CrossRef]

37. J. R. Zurita-Sánchez and P. Halevi, “Resonances in the optical response of a slab with time-periodic dielectric function ε (t),” Phys. Rev. A 81, 053834 (2010). [CrossRef]

38. J. Reyes-Ayona and P. Halevi, “Observation of genuine wave vector (k or β) gap in a dynamic transmission line and temporal photonic crystals,” Appl. Phys. Lett. 107, 074101 (2015). [CrossRef]

39. J. S. Martínez-Romero, O. Becerra-Fuentes, and P. Halevi, “Temporal photonic crystals with modulations of both permittivity and permeability,” Phys. Rev. A 93, 063813 (2016). [CrossRef]

40. J. R. Reyes-Ayona and P. Halevi, “Electromagnetic wave propagation in an externally modulated low-pass transmission line,” IEEE Transactions on Microw. Theory Tech. 64, 3449–3459 (2016). [CrossRef]

41. S. Taravati, “Aperiodic space-time modulation for pure frequency mixing,” Phys. Rev. B 97, 115131 (2018). [CrossRef]

42. D. Ramaccia, D. L. Sounas, A. Alù, A. Toscano, and F. Bilotti, “Doppler cloak restores invisibility to objects in relativistic motion,” Phys. Rev. B 95, 075113 (2017). [CrossRef]

43. A. Mock, D. Sounas, and A. Alù, “Tunable orbital angular momentum radiation from angular-momentum-biased microcavities,” Phys. Rev. Lett. 121, 103901 (2018). [CrossRef] [PubMed]

44. M. F. Yanik, W. Suh, Z. Wang, and S. Fan, “Stopping light in a waveguide with an all-optical analog of electromagnetically induced transparency,” Phys. review letters 93, 233903 (2004). [CrossRef]

45. M. F. Yanik and S. Fan, “Time reversal of light with linear optics and modulators,” Phys. review letters 93, 173903 (2004). [CrossRef]

46. M. Minkov and S. Fan, “Localization and time-reversal of light through dynamic modulation,” Phys. Rev. B 97, 060301 (2018). [CrossRef]

47. M. Moharam, E. B. Grann, D. A. Pommet, and T. Gaylord, “Formulation for stable and efficient implementation of the rigorous coupled-wave analysis of binary gratings,” JOSA a 12, 1068–1076 (1995). [CrossRef]

48. V. Liu and S. Fan, “S⁴ : A free electromagnetic solver for layered periodic structures,” Comput. Phys. Commun. 183, 2233 – 2244 (2012). [CrossRef]

49. P. Lalanne and G. M. Morris, “Highly improved convergence of the coupled-wave method for tm polarization,” JOSA A 13, 779–784 (1996). [CrossRef]

50. “LAPACK–Linear Algebra PACKage,” http://www.netlib.org/lapack/. Accessed: October 1, 2018.

51. A. Taflove and S. Hagness, Computational Electrodynamics: The Finite-difference Time-domain Method, Artech House Antennas and Prop (Artech House, 2005).

52. B. Szabó, B. A. Szabo, and I. Babuška, Finite Element Analysis(John Wiley & Sons, 1991).

53. J. Volakis, Integral equation methods for electromagnetics (The Institution of Engineering and Technology, 2012).

54. X. Ni, Z. Liu, A. Boltasseva, and A. V. Kildishev, “The validation of the parallel three-dimensional solver for analysis of optical plasmonic bi-periodic multilayer nanostructures,” Appl. Phys. A 100, 365–374 (2010). [CrossRef]

55. A. O. Korotkevich, X. Ni, and A. V. Kildishev, “Fast eigensolver for plasmonic metasurfaces,” Opt. Mater. Express 4, 288–299 (2014). [CrossRef]

56. J. Cheng, S. Inampudi, and H. Mosallaei, “Optimization-based dielectric metasurfaces for angle-selective multifunctional beam deflection,” Sci. reports 7, 12228 (2017). [CrossRef]

57. S. Inampudi, V. Toutam, and S. Tadigadapa, “Robust visibility of graphene monolayer on patterned plasmonic substrates,” Nanotechnology 30, 015202 (2019). [CrossRef]

58. J. Mayer, R. Baron, and O. Marsh, “Observation of double injection in long silicon p-i-n structures,” Phys. Rev. 137, A286 (1965). [CrossRef]

59. S. Colburn, A. Zhan, and A. Majumdar, “Tunable metasurfaces via subwavelength phase shifters with uniform amplitude,” Sci. reports 7, 40174 (2017). [CrossRef]

60. A. Forouzmand and H. Mosallaei, “All-dielectric c-shaped nanoantennas for light manipulation: Tailoring both magnetic and electric resonances to the desire,” Adv. Opt. Mater. 5, 1700147 (2017). [CrossRef]

61. Y. Yang, A. E. Miroshnichenko, S. V. Kostinski, M. Odit, P. Kapitanova, M. Qiu, and Y. S. Kivshar, “Multimode directionality in all-dielectric metasurfaces,” Phys. Rev. B 95, 165426 (2017). [CrossRef]

62. R. Magnusson, “Wideband reflectors with zero-contrast gratings,” Opt. letters 39, 4337–4340 (2014). [CrossRef]

63. S. Taravati, “Giant linear nonreciprocity, zero reflection, and zero band gap in equilibrated space-time-varying media,” Phys. Rev. Appl. 9, 064012 (2018). [CrossRef]

64. “HORIBA FHR Spectrometer Brochure,” http://www.horiba.com/fileadmin/uploads/Scientific/Documents/Mono/FHR6401000.pdf . Accessed: October 1, 2018.

65. C. A. Brackett, “Dense wavelength division multiplexing networks: Principles and applications,” IEEE J. on Sel. Areas Commun. 8, 948–964 (1990). [CrossRef]

66. A. M. E. A. Shaltout, “Photonic metasurfaces for spatiotemporal and ultrafast light control,” Ph.D. thesis, Purdue University (2015).

67. A. Shaltout, K. Lagoudakis, J. van de Groep, S. Kim, J. Vučković, V. Shalaev, and M. Brongersma, “Continuous angle beamsteeringusingspatiotemporal frequency-combcontrolindielectric metasurfaces,” in CLEO: QELS_Fundamental Science, (Optical Society of America, 2018), pp. FW3H–2.

68. L. Falkovsky and A. Varlamov, “Space-time dispersion of graphene conductivity,” The Eur. Phys. J. B 56, 281–284 (2007). [CrossRef]

69. L. Ju, B. Geng, J. Horng, C. Girit, M. Martin, Z. Hao, H. A. Bechtel, X. Liang, A. Zettl, Y. R. Shen, et al., “Graphene plasmonics for tunable terahertz metamaterials,” Nat. nanotechnology 6, 630 (2011). [CrossRef]

70. T. Yasui, E. Saneyoshi, and T. Araki, “Asynchronous optical sampling terahertz time-domain spectroscopy for ultrahigh spectral resolution and rapid data acquisition,” Appl. Phys. Lett. 87, 061101 (2005). [CrossRef]

71. X. Liu and D. A. McNamara, “The use of the fdtd method for electromagnetic analysis in the presence of indepedently time-varying media,” Int. J. Infrared Millim. Waves 28, 759–778 (2007). [CrossRef]

72. T. J. Smy and S. Gupta, “Exact finite-difference time-domain modelling of broadband huygens’ metasurfaces with lorentz dispersions,” arXiv preprint arXiv :1609.05575 (2016).

73. S. A. Stewart, T. J. Smy, and S. Gupta, “Finite-difference time-domain (fdtd) modelling of space-time modulated metasurfaces,” IEEE Transactions on Antennas Propag. (2017).

Rigorous space-time coupled-wave analysis for patterned surfaces with temporal permittivity modulation [Invited]

Abstract

1. Introduction

2. Matrix formulation for rigorous space-time coupled-wave analysis (RSTWA)

3. Numerical results and discussion

3.1. Validation

3.2. Anomalous beam-steering with spatially phase-variant time-modulated dielectric metasurfaces

3.3. Dynamic beam generation with spatially frequency-variant time-modulated plasmonic metasurfaces

4. Conclusion

Appendix I Derivation of eigenvalue equation for the field profiles

Appendix II Time-variant FDTD

Funding

References

Cited By

Figures (8)

Equations (42)

Optical Materials Express

Mohammad Mahdi Salary	https://orcid.org/0000-0003-2056-3348
Samad Jafar-Zanjani	https://orcid.org/0000-0002-8818-4021