Stacked Time-Varying Metasurfaces

TL;DR

This paper presents a novel approach for designing stacked spatiotemporal metasurfaces using independently engineered uniform layers with periodic modulation, enabling advanced wave control with fewer elements and broad practical applications.

Contribution

It introduces a general theoretical framework for designing and optimizing stacked metasurfaces with independent modulation functions, reducing complexity and enhancing functionality.

Findings

Achieved magnet-free isolation at fundamental frequency.

Demonstrated a temporal analogue of circulators.

Implemented a metasurface stack with few modulated elements in a waveguide.

Abstract

Spatiotemporal metasurfaces offer unique opportunities for wave manipulation, however, their practical realization is often constrained by the requirement for in-plane spatial modulation, which necessitates a large number of time-varying elements. In this work, we introduce an alternative architecture based on a cascade of spatially uniform metasurfaces subjected to periodic temporal modulation. Although all metasurfaces share the same modulation frequency, their individual modulation functions are independently engineered to achieve a desired complex electromagnetic response. We develop a general theoretical framework for the design and optimization of such stacked metasurface systems, composed of dense arrays of cylindrical meta-atoms with time-varying plasma and/or collision frequencies. The effectiveness of the approach is demonstrated through the optimization of metasurface designs that enable magnet-free isolation at the fundamental frequency and a temporal analogue of circulators. Furthermore, we show that a metasurface stack can be implemented using only a few time-modulated elements embedded within a parallel-plate waveguide, opening new avenues for extremely compact, versatile, and scalable spatiotemporal platforms for next-generation photonic and microwave systems.

Equations

$$\frac{{\mathrm{d}}^2\mathbf{P}(t)}{\mathrm{d}{t}^2}+\gamma (t)\frac{\mathrm{d}\mathbf{P}(t)}{\mathrm{d}t}={\omega }_{\mathrm{p}}^2(t){\epsilon }_0\mathbf{E}(t),$$

$$\frac{\mathrm{d}{\mathbf{J}}_{\mathrm{m}}(t)}{\mathrm{d}t}+\gamma (t){\mathbf{J}}_{\mathrm{m}}(t)={\omega }_{\mathrm{p}}^2(t){ϵ}_0\mathbf{E}(t).$$

$$\frac{1}{{\omega }_{\mathrm{p}}^2(t){ϵ}_{0}S}\frac{\mathrm{d}{I}_{\mathrm{m}}(t)}{\mathrm{d}t}+\frac{\gamma (t)}{{\omega }_{\mathrm{p}}^2(t){ϵ}_{0}S}{I}_{\mathrm{m}}(t)=E(t).$$

$$\frac{\mathrm{d}}{\mathrm{d}t}\left[L_{\mathrm{m}}(t)I_{\mathrm{m}}(t)\right]+R_{\mathrm{m}}(t)I_{\mathrm{m}}(t)=v(t)/\Delta z,$$

$$L_{\mathrm{m}}(t)=\frac{1}{{\omega }_{\mathrm{p}}^2(t){ϵ}_{0}S},R_{\mathrm{m}}(t)=\frac{\gamma (t)}{{\omega }_{\mathrm{p}}^2(t){ϵ}_{0}S}-\frac{\mathrm{d}}{\mathrm{d}t}\left[\frac{1}{{\omega }_{\mathrm{p}}^2(t){ϵ}_{0}S}\right].$$

$$Z_{\mathrm{g}}=Z_{\mathrm{m}}d+j\frac{{\eta }_0}{2}{\alpha }_{\mathrm{ABC}},$$

$$R_{\mathrm{g}}(t)=R_{\mathrm{m}}(t)d,$$

$$L_{\mathrm{g}}(t)=\left[L_{\mathrm{m}}(t)+\frac{{\eta }_0}{2}\frac{\sqrt{{ϵ}_0{\mu }_0}}{\pi }\ln \left(\frac{d}{2\pi r_0}\right)\right]d.$$

$${\omega }_{\mathrm{p},i}(t)$$

$${\gamma }_i(t)$$

$$R_{\mathrm{g},i}(t)=\sum _{q=-\infty }^{+\infty }{r}_{q,i}{e}^{jq{\omega }_{\mathrm{m}}t}$$

$$L_{\mathrm{g},i}(t)=\sum _{q=-\infty }^{+\infty }{l}_{q,i}{e}^{jq{\omega }_{\mathrm{m}}t}.$$

$$\begin{array}[]{l}\displaystyle I_{\rm L}(t)=I_{\rm R}(t)=\sum_{n=-\infty}^{+\infty}i_{n}^{\rm L}e^{j\omega_{n}t},\vskip 5.69054pt\\ \displaystyle V_{\rm L,R}(t)=\sum_{n=-\infty}^{+\infty}v_{n}^{\rm L,R}e^{j\omega_{n}t}.\end{array}$$

$${\mathbf{v}}^{\mathrm{L},\mathrm{R}}=[\cdots ,v_{-1}^{\mathrm{L},\mathrm{R}},v_0^{\mathrm{L},\mathrm{R}},v_{+1}^{\mathrm{L},\mathrm{R}},\cdots {]}^T\text{and}{\mathbf{i}}^{\mathrm{L}}=[\cdots ,i_{-1}^{\mathrm{L}},i_0^{\mathrm{L}},i_{+1}^{\mathrm{L}},\cdots {]}^T,$$

$$I_{\mathrm{L}}(t)L_{\mathrm{g}}(t)={\int }^t{V}_{\mathrm{L}}(t)\mathrm{d}t.$$

$$\left(\begin{array}[]{cccc}jl_{0}\omega_{-N}&jl_{-1}\omega_{-N}&\cdots&jl_{-2N}\omega_{-N}\\ jl_{1}\omega_{1-N}&jl_{0}\omega_{1-N}&\cdots&\vdots\\ \vdots&\vdots&\ddots&\vdots\\ jl_{2N}\omega_{N}&jl_{2N-1}\omega_{N}&\cdots&jl_{0}\omega_{N}\end{array}\right)\cdot\left(\begin{array}[]{c}i^{\rm L}_{-N}\\ i^{\rm L}_{1-N}\\ \vdots\\ i^{\rm L}_{N}\end{array}\right)=\left(\begin{array}[]{c}v^{\rm L}_{-N}\\ v^{\rm L}_{1-N}\\ \vdots\\ v^{\rm L}_{N}\end{array}\right),$$

$$V_{\mathrm{R}}(t)=R_{\mathrm{g}}(t)I_{\mathrm{L}}(t).$$

$$\left(\begin{array}[]{cccc}r_{0}&r_{-1}&\cdots&r_{-2N}\\ r_{1}&r_{0}&\cdots&\vdots\\ \vdots&\vdots&\ddots&\vdots\\ r_{2N}&r_{2N-1}&\cdots&r_{0}\end{array}\right)\cdot\left(\begin{array}[]{c}i^{\rm L}_{-N}\\ i^{\rm L}_{1-N}\\ \vdots\\ i^{\rm L}_{N}\end{array}\right)=\left(\begin{array}[]{c}v^{\rm R}_{-N}\\ v^{\rm R}_{1-N}\\ \vdots\\ v^{\rm R}_{N}\end{array}\right),$$

$$\left(\begin{array}{c}{\mathbf{v}}_{\text{tot}}^{\text{inp}}\\ {\mathbf{i}}_{\text{tot}}^{\text{inp}}\end{array}\right)=\left(\begin{array}{cc}\overline{\overline{A}}& \overline{\overline{B}}\\ \overline{\overline{C}}& \overline{\overline{D}}\end{array}\right)\cdot \left(\begin{array}{c}{\mathbf{v}}_{\text{tot}}^{\text{out}}\\ {\mathbf{i}}_{\text{tot}}^{\text{out}}\end{array}\right)=\overline{\overline{T}}\cdot \left(\begin{array}{c}{\mathbf{v}}_{\text{tot}}^{\text{out}}\\ {\mathbf{i}}_{\text{tot}}^{\text{out}}\end{array}\right).$$

$$\overline{\overline{T}}_{\mathrm{m},i}=\left(\begin{array}[]{cc}\overline{\overline{I}}&0\\ \overline{\overline{Y}}_{\mathrm{g},i}&\overline{\overline{I}}\end{array}\right),$$

$$\overline{\overline{T}}_{{\rm a},i}=\left(\begin{array}[]{cc}\overline{\overline{I}}{\rm cos}\,k_{n}D_{i}&\overline{\overline{I}}j\eta_{0}\,{\rm sin}\,k_{n}D_{i}\\ \overline{\overline{I}}j\eta_{0}^{-1}\,{\rm sin}\,k_{n}D_{i}&\overline{\overline{I}}{\rm cos}\,k_{n}D_{i}\end{array}\right).$$

$${\overline{\overline{T}}}_{\mathrm{F}}={\overline{\overline{T}}}_{\mathrm{m},1}\cdot {\overline{\overline{T}}}_{\mathrm{a},1}\cdot {\overline{\overline{T}}}_{\mathrm{m},2}\cdot {\overline{\overline{T}}}_{\mathrm{a},2}\cdot \dots \cdot {\overline{\overline{T}}}_{\mathrm{a},M-1}\cdot {\overline{\overline{T}}}_{\mathrm{m},M},$$

$${\overline{\overline{T}}}_{\mathrm{B}}={\overline{\overline{T}}}_{\mathrm{m},M}\cdot {\overline{\overline{T}}}_{\mathrm{a},\mathrm{M}-1}\cdot {\overline{\overline{T}}}_{\mathrm{m},M-1}\cdot {\overline{\overline{T}}}_{\mathrm{a},\mathrm{M}-2}\cdot \dots \cdot {\overline{\overline{T}}}_{\mathrm{a},1}\cdot {\overline{\overline{T}}}_{\mathrm{m},1},$$

$$\left(\begin{array}[]{c}{\bf v}_{\rm i}+{\bf v}_{\rm r}\\ \eta_{0}^{-1}({\bf v}_{\rm i}-{\bf v}_{\rm r})\end{array}\right)=\overline{\overline{T}}_{\rm F}\cdot\left(\begin{array}[]{c}{\bf v}_{\rm t}\\ \eta_{0}^{-1}{\bf v}_{\rm t}\end{array}\right)=\begin{pmatrix}\overline{\overline{A}}_{\rm F}&\overline{\overline{B}}_{\rm F}\\ \overline{\overline{C}}_{\rm F}&\overline{\overline{D}}_{\rm F}\end{pmatrix}\cdot\left(\begin{array}[]{c}{\bf v}_{\rm t}\\ \eta_{0}^{-1}{\bf v}_{\rm t}\end{array}\right).$$

$${\overline{\overline{\tau }}}_{\mathrm{F}}=2{\left({\overline{\overline{A}}}_{\mathrm{F}}+{\eta }_0^{-1}{\overline{\overline{B}}}_{\mathrm{F}}+{\eta }_0{\overline{\overline{C}}}_{\mathrm{F}}+{\overline{\overline{D}}}_{\mathrm{F}}\right)}^{-1}$$

$${\overline{\overline{\Gamma }}}_{\mathrm{F}}=\left({\overline{\overline{A}}}_{\mathrm{F}}+{\eta }_0^{-1}{\overline{\overline{B}}}_{\mathrm{F}}\right){\overline{\overline{\tau }}}_{\mathrm{F}}-\overline{\overline{I}}.$$

$${\omega }_{\mathrm{p},1}(t)$$

$${\gamma }_1(t)$$

$${\omega }_{\mathrm{p},2}(t)$$

$${\gamma }_2(t)$$

$${\Delta }_{\mathrm{is}}=\left|\left|E_{\mathrm{t},0}^{\mathrm{F}}\right|-1\right|+\left|\left|E_{\mathrm{t},0}^{\mathrm{B}}\right|-0\right|.$$

$${\Delta }_1$$

$${\Delta }_2$$

$${\Delta }_3$$

Full text

Mostafa Movahediqomi Sergei Tretyakov Viktar Asadchy* Xuchen Wang

\dedication

{affiliations}

Mostafa Movahediqomi, Sergei Tretyakov, Viktar Asadchy

Department of Electronics and Nanoengineering

Aalto University

15500, FI-00076 Aalto, Finland.

Email Address: viktar.asadchy@aalto.fi

Xuchen Wang

Qingdao Innovation and Development Base

Harbin Engineering University

266400, Qingdao, China.

1 Introduction

In recent years, the pursuit of artificial materials with tailored electromagnetic responses has led to the rapid development of metamaterials and, more recently, metasurfaces [1]. Traditionally, such structures were designed exclusively through spatial engineering, by modifying the geometry, size, or arrangement of the constituent scatterers. Yet, this approach is inherently restricted. For time-invariant and linear scatterers, the scattered frequency and total energy are conserved, preventing access to many desirable functionalities. These limitations have spurred growing interest in exploiting the temporal dimension of material parameters, which provides a powerful route to fundamentally new scattering behaviors [2, 3, 4].

Time-varying platforms have already revealed the potential to unlock unconventional wave phenomena, such as unbounded field accumulation [5] and temporal analogues of Wood’s anomalies [6]. In this context, spatiotemporal metasurfaces have emerged as a particularly promising architecture, enabling frequency conversion [7, 8], isolation [9, 10], circulation [11, 9], power combination [12], and nonreciprocal phase manipulation [9]. The central principle is the induction of a synthetic motion along the metasurface plane, which effectively breaks time-reversal symmetry [13]. However, realizing such metasurfaces remains technically demanding. Approaches based on continuous space-time gradients require carefully engineered impedance profiles [14], while discrete space-time metagratings rely on dense arrays of synchronized meta-atoms [15]—both posing major implementation barriers due to the in-plane spatial modulation requiring non-local design and many time-varying elements.

In this work, we introduce a different strategy that mitigates these challenges by shifting the spatial modulation into the out-of-plane direction. Our system consists of a stack of metasurfaces that are uniform across their planes but undergo periodic modulation in time. This longitudinal arrangement enhances the design freedom compared to conventional in-plane spatiotemporal metasurfaces, while providing easier practical realization. Our theoretical modeling is performed in terms of bulk parameters of the meta-atoms, such as plasma and collision frequencies, with explicit consideration of frequency dispersion—an aspect particularly advantageous for optical realizations. Conceptually, the proposed stacked structures resemble traveling-wave modulation in transmission-line and optical systems [16, 17] or tandem-type modulation [18], yet they are more general and substantially more compact. Furthermore, these metasurface cascades represent a generalization of recently proposed single-layer arrays of time-modulated spherical meta-atoms [19, 20]. We showcase their versatility by designing stacked metasurfaces capable of magnet-free isolation at the fundamental frequency and a temporal analogue of circulators. Remarkably, both designed metasurface stacks have electrical thickness below the operational wavelength. Furthermore, we demonstrate that such architectures can be straightforwardly implemented in parallel-plate waveguides with only a few time-modulated elements.

2 Theoretical framework

Let us consider the configuration of stacked time-varying metasurfaces shown in Fig. 1(a), illuminated by a plane wave along the normal direction. Each metasurface comprises a dense array of long cylinders with material properties varying periodically in time. The cylinders have extent along the $z$-direction much larger than the operational wavelength and, therefore, are modeled as infinitely long. The incident-wave electric field is polarized along the cylinders. The cylinder radii $r_{0}$ and the spacing between the neighboring cylinders $d$ are both assumed to be sub-wavelength to ensure that each metasurface is effectively uniform in the $xz$-plane. The separation distance between the $i$-th and $(i+1)$-th metasurfaces is denoted as $D_{i}$ and can be arbitrary but larger than $d$ to ensure that there is no near-field coupling between the cylinders in neighboring metasurfaces [21]. We model the material of each cylinder in the $i$-th metasurface with Drude dispersion with periodically time-varying plasma frequency $\omega_{{\rm p},i}^{2}(t)=\frac{N_{i}(t)e^{2}}{m\epsilon_{0}}$ ($N_{i}(t)$ is the time-dependent electron density, $m$ and $e$ are the effective electron mass and electric charge) and periodically time-varying electron collision frequency $\gamma_{i}(t)$. It is important to note that since such cylinder arrays can be considered as impedance sheets, other physical realizations are equally possible.

2.1 Modeling a single time-varying cylindrical meta-atom

We begin by analyzing a single cylinder with time-varying material properties. Here, for brevity, we omit the use of the subscript $i$ denoting which metasurface in the stack the considered meta-atom belongs to. The dynamics of its polarization density ${\bf P}$ is governed by

[TABLE]

where ${\bf E}(t)$ is the total electric field in the cylindrical meta-atom. Since the current density inside the meta-atom is related to its bulk electric polarization as ${\bf J}_{\rm m}(t)=\frac{{\rm d}{\bf P}(t)}{{\rm d}t}$, Eq. (1) can be rewritten as

[TABLE]

Assuming that the current density is approximately uniform over the cylinder cross section, the total current in the meta-atom is the product of the current density and the cross-sectional area $S=\pi r_{0}^{2}$, that is, $I_{\rm m}(t)=J_{\rm m}(t)S$. Substituting this into Eq. (2) and taking the projection of the equation to the $z$-axis yields:

[TABLE]

We can notice that Eq. (3) can be interpreted as an equation describing voltage and current in a series connection of time-varying resistance $R_{\rm m}(t)$ and time-varying inductance $L_{\rm m}(t)$:

[TABLE]

where $\Delta z$ is the differential length of the cylinder along the $z$-direction. Therefore, a single time-varying cylindrical meta-atom can be modeled by a series connection of time-varying inductance and resistance given by

[TABLE]

Note that the physical meaning of $L_{\rm m}$ and $R_{\rm m}$ are the resistance and inductance per unit length of the cylindrical meta-atom, with units of [H/m] and [$\Omega$/m], respectively. The equivalent circuit of a differential length $\Delta z$ of the meta-atom is shown in Fig. 1(a).

2.2 Modeling an array of time-varying cylindrical meta-atoms

The above relations were derived for a single meta-atom. We now extend the analysis to a dense array of such meta-atoms with periodicity $d$. If the meta-atoms were time-invariant, the grid impedance of such an array would be given by [21, Eq. (4.39)]

[TABLE]

where $Z_{\rm m}=R_{\rm m}+j\omega L_{\rm m}$ is the complex impedance per unit length of a single meta-atom defined in the frequency domain, $j$ is the imaginary unit (we adopt $e^{j\omega t}$ notation for harmonic oscillations), and $\eta_{0}=\sqrt{\mu_{0}/\varepsilon_{0}}$ is the wave impedance of the background material (here assumed to be free space). Parameter $\alpha_{\rm ABC}=\dfrac{k_{0}d}{\pi}\ln\left(\dfrac{d}{2\pi r_{0}}\right)$ is the grid parameter established by Kontorovich et al. [22] in the framework of averaged boundary conditions (here $k_{0}$ is the wavenumber of the background material). It is related to the frequency and geometry of the array, but independent of the material properties. The additional inductive term in Eq. (6) arises purely from the geometry.

For the array of time-varying meta-atoms, this geometric contribution must also be included. Therefore, the equivalent grid resistance and inductance of the metasurface is given by

[TABLE]

Note that here $R_{\mathrm{g}}(t)$ and $L_{\mathrm{g}}(t)$ have units of [$\Omega$] and [H]. Thus, the stacked metasurfaces can be modeled by a cascade of shunt resistances and inductances in an equivalent transmission-line model shown in Fig. 1(b). The spacing between neighboring metasurfaces is modeled by a finite-length section of a transmission line with characteristic impedance equal to the wave impedance of the background medium. In Supplementary Material, we compare the transmission-line model developed here with the full-wave simulations of the metasurface.

2.3 Mode-matching method

So far, we have considered a generic, not necessarily periodic, temporal variation of material properties within the meta-atoms. We now focus on the special case of periodic temporal modulations of material parameters. Specifically, we aim to analyze the reflection and transmission of different frequency harmonics from a single metasurface subject to a time-periodic modulation of the plasma frequency $\omega_{\rm p}$ and the collision frequency $\gamma$ of its constituent meta-atoms. Given the temporal profiles of $\omega_{\rm p}$ and $\gamma$, we derive the corresponding time-varying grid resistance $R_{\rm g}$ and inductance $L_{\rm g}$. Employing the mode-matching method [23], we subsequently determine the ABCD (transfer) matrix of the time-varying metasurface. This matrix formalism can be readily extended to compute the reflection and transmission characteristics of arbitrarily complex stacks composed of multiple metasurfaces.

The material parameters $\omega_{{\rm p}}$ and $\gamma$ of different metasurfaces and their temporal modulation profiles can be chosen independently. Nevertheless, in our analysis, they all share the same fundamental modulation frequency $\omega_{\rm m}$, resulting in the fact that the entire stack of metasurfaces has a temporal periodicity $T=2\pi/\omega_{\rm m}$. Therefore, we expand the material parameters of the $i$-th metasurface into the Fourier series:

[TABLE]

where $w_{q,i}$ and $g_{q,i}$ determine the harmonic coefficients for the $i$-th metasurface.

Since $R_{{\rm g},i}(t)$ and $L_{{\rm g},i}(t)$ are the functions of $\omega_{{\rm p},i}(t)$ and $\gamma_{i}(t)$, they are also periodic with the same period $T$ and can be written as

[TABLE]

The harmonic coefficients $r_{q,i}$ and $l_{q,i}$ can be expressed in terms of $w_{q,i}$ and $g_{q,i}$ using Eq. (5). Below, for brevity, we omit the metasurface index “$i$”, with all quantities henceforth referring to the same chosen metasurface.

The modulation induces an infinite number of harmonics with the frequencies $\omega_{n}=\omega_{0}+n\omega_{\rm m}$, where $\omega_{0}$ is the frequency of the incident on the metasurface plane wave. Therefore, the current going through and the voltage across the equivalent resistor and inductor can be written as

[TABLE]

In conventional static circuits, the voltage and current are scalar values. Moreover, the impedance of lumped components, which relates the scalar current and voltage, is a scalar. However, in this multi-mode system, the voltage and current can be written in the form of vectors, i.e.,

[TABLE]

where each element in the vector denotes the complex amplitude of each harmonic. The impedance or admittance of the modulated components should be a matrix which relates the current and voltage vectors to one another. In order to calculate the scattering harmonics, we should know the matrix representations of the resistor and inductor under modulation.

The current-voltage relation of the time-varying inductance in the time domain is given by

[TABLE]

By substituting Eqs. (10), and (9b) in Eq. (11), we can obtain the corresponding equation in matrix form as $\overline{\overline{Z}}_{\rm L}\cdot{\bf i}^{\rm L}={\bf v}^{\rm L}$ or explicitly as [23]

[TABLE]

when truncating the Floquet harmonics $n$ beyond the range from $-N$ to $N$.

For the resistance, the relation between the current and voltage of harmonics has a simple form:

[TABLE]

Following the same procedure, we obtain the following matrix equation:

[TABLE]

which can be written in a short form as $\overline{\overline{Z}}_{\rm R}\cdot{\bf i}^{\rm L}={\bf v}^{\rm R}$.

The total impedance matrix of the metasurface is a summation of the inductance and resistance matrices, that is, $\overline{\overline{Z}}_{\mathrm{g}}=\overline{\overline{Z}}_{\rm L}+\overline{\overline{Z}}_{\rm R}$.

2.4 Transfer-matrix method

For the general case when multiple metasurfaces are arranged in a cascade, it is convenient to define the transfer (ABCD) matrix that relates the currents and voltages at the input and output of a given metasurface as

[TABLE]

The transfer matrix for the $i$-th metasurface can be described as [24, Sec. 4.4]:

[TABLE]

where $\overline{\overline{I}}$ is the $(2N+1)\times(2N+1)$ unity matrix and $\overline{\overline{Y}}_{\mathrm{g},i}=\overline{\overline{Z}}_{\mathrm{g},i}^{-1}$.

The contribution of the air gaps between the neighboring metasurfaces should also be taken into account in the calculation. These gaps are equivalent to transmission-line sections with lengths $D_{i}$ and wavenumbers for each Floquet harmonic $k_{n}=\omega_{n}/c$, where $c$ is the speed of light in free space. The transfer matrix of the $i$-th air gap is given by [24, Sec. 4.4]

[TABLE]

Eventually, the transfer matrix of the entire metasurface stack for forward and backward wave illuminations can be calculated by a multiplication of the transfer matrices:

[TABLE]

[TABLE]

where $M$ is the number of stacked metasurfaces.

Finding the transmission and reflection coefficients of the metasurface cascade, we should remember that the transfer matrix relates the total voltages and currents at the input and output. In other words, at the input, they include both incident and reflected waves, while at the output, only the transmitted waves:

[TABLE]

The transmission and reflection coefficients for forward illumination scenario, defined as ${\bf v}_{\rm t}=\overline{\overline{\tau}}_{\rm F}\cdot{\bf v}_{\rm i}$ and ${\bf v}_{\rm r}=\overline{\overline{\Gamma}}_{\rm F}\cdot{\bf v}_{\rm i}$, can be derived as [24, Sec. 4.4]

[TABLE]

[TABLE]

Likewise, one can obtain similar relations for the transmission and reflection coefficients for the backward illumination scenario $\overline{\overline{\tau}}_{\rm B}$ and $\overline{\overline{\Gamma}}_{\rm B}$.

In summary, the developed theoretical formalism allows us to calculate complex amplitudes of all reflected and transmitted Floquet harmonics as described by Eqs. (20) and (21) for a known material modulation given by Eqs. (8a) and (8b). In the next section, we apply this formalism to an optimization algorithm to determine optimal parameters of temporal modulation of each metasurface for obtaining desired reflection and transmission responses of the metasurface cascade.

3 Demonstrators

Although the proposed theoretical framework is general and applicable to metasurface stacks with an arbitrary number of layers $M$ and arbitrary inter-metasurface distances $D_{i}$, for simplicity, we focus in this section on demonstrator designs, each consisting of two metasurfaces. As a form of temporal modulation in both metasurfaces, we assume a sinusoidal variation in the plasma and collision frequencies:

[TABLE]

The coefficients defined in these equations can be uniquely converted to those in Eqs. (8a) and (8b). Without loss of generality, here we assume that the relative modulation strengths $\delta_{1,2}$ and modulation phases of both plasma and collision frequencies are the same for a given metasurface. We also assume a phase difference $\phi$ between the material properties of the two metasurfaces. Coefficients shown with the bars on top denote the time-average (mean) values of the corresponding material parameters.

3.1 Isolation at the fundamental frequency

The quest of designing isolators using time-varying systems has been an active field of recent research [25, 26, 27]. Efficient isolation (nonreciprocal transmission) has been typically achieved in bulk systems where the electrical thickness is much larger than the operational wavelength [17, 28]. The challenge of obtaining isolation in compact time-varying systems is difficult to overcome, despite relying on simultaneous frequency conversion [29, 30] or the use of bianisotropy [23]. The former feature is undesired, as many applications require isolation at the fundamental (the same as the incident) frequency. The latter approach, although known, has not had practical realization to date. Here, we realize isolation at the fundamental frequency in the stacked metasurface geometry depicted in Fig. 1(a). As we show below, only two metasurfaces inside the stack are sufficient to obtain isolation with high efficiency and without frequency conversions. This two-metasurface cascade has an effective time-varying bianisotropic response [3].

For the optimization of material parameters, the fmincon optimizer in MATLAB is selected due to its ability to handle nonlinear constraints and its compatibility with the “MultiStart” setting, which enhances convergence speed and quality. The optimized parameters include those listed in Eqs. (22a)–(22d), that is, $\bar{w}_{1,2}$, $\bar{\gamma}_{1,2}$, $\delta_{1,2}$, $\phi$, and $\omega_{\rm m}$, as well as the distance between the layers $D_{1}\gg d$. Without loss of generality, we fix $d=\lambda_{0}/20$ and the radius of the cylinders $r_{0}=\lambda_{0}/150$, where $\lambda_{0}$ is the wavelength in free space at the fundamental (incident) frequency $\omega_{0}$.

The optimization objective is to reach unit and zero amplitudes of the transmitted waves for the forward and backward illuminations, respectively, at the fundamental frequency $\omega_{0}$. The incident wave is represented by a single frequency $\omega_{0}$ with an electric field amplitude of 1 V/m in both illumination scenarios. As mentioned before, the reflected and transmitted propagating waves in the considered system shown in Fig. 1(a) can be interpreted as the corresponding voltages in the equivalent transmission-line model, as depicted in Fig. 1(b). In other words, we can relate harmonics of the reflected and transmitted electric fields to the incident fields using the same scattering coefficients as calculated from Eqs. (21) and (20), i.e., ${\bf E}_{\rm t}^{\rm F}=\overline{\overline{\tau}}^{\rm F}\cdot{\bf E}_{\rm i}^{\rm F}$, and ${\bf E}_{\rm t}^{\rm B}=\overline{\overline{\tau}}^{\rm B}\cdot{\bf E}_{\rm i}^{\rm B}$. Thus, we define the cost function in the optimization as

[TABLE]

In this way, the optimizer tries to reach unit forward transmission at the fundamental frequency $\omega_{0}$, while minimizing the backward transmission for the same frequency harmonic. We truncate the range of Floquet harmonics to $n\in[-N;N]$ with $N=10$. The optimized material parameters are listed in Supplementary Material.

The amplitudes of different transmitted harmonics for the optimized metasurface stack are shown in Fig. 2(a). The transmitted harmonics for the forward illumination are shown with blue stems. We can observe a spectrum with a Gaussian-type frequency distribution where the amplitude of the transmitted $\omega_{0}$ harmonic is very close to 1. On the contrary, for the backward illumination (red stems), the transmission is strongly suppressed at all frequencies, with the transmitted electric field at $\omega_{0}$ being equal to $0.021$ V/m. The value of the cost function represented in Eq. (23) after optimization is $0.0211$. Thus, the isolation ratio at the fundamental frequency reaches $33.56~\rm{dB}$. It is worth noting that, depending on the application, sideband frequency harmonics ($\omega_{n}$ with $n\neq 0$) possessing non-zero amplitudes in the forward illumination scenario may be undesirable. In this case, they can be easily suppressed by incorporating a conventional frequency band-pass filter into the metasurface stack.

3.2 Temporal circulation

To further demonstrate the versatility of the proposed stacked metasurfaces, we design a novel and exotic system. In microwave and optical engineering, one of the fundamental nonreciprocal components is the circulator—a device that operates at a single frequency but routes signals nonreciprocally among three distinct spatial ports in a unidirectional, cyclic manner [24]. In this section, we propose a temporal analogue of the circulator. Unlike conventional circulators, this device features a single spatial port but three distinct frequency ports, enabling the wave to circulate nonreciprocally across these frequencies in a manner analogous to spatial circulators. We term this device a temporal (or frequency) circulator. Potential applications include compact nonreciprocal light guiding, the realization of synthetic dimensions [31, 8], and dynamic spectral multiplexing for integrated photonic systems.

To ensure that there is only one spatial port, we incorporate a perfectly conducting plate (ground plane) behind the metasurface stack, as shown in the inset of Fig. 2(b). This way, we truncate half of the space below the plate and consider the normally incident and reflected waves to share the same single spatial port. For the three frequency (temporal) ports, we consider three frequencies coupled through the temporal modulations in the metasurfaces: $\omega_{0}$ and $\omega_{\pm 1}=\omega_{0}\pm\omega_{\rm m}$. The distance between the conducting plate and the closer metasurface is denoted as $D_{2}$. In the transfer-matrix method, we model the conducting plate as a parallel short circuit with the ABCD matrix $\overline{\overline{T}}_{\rm pl}$ being equal to $\overline{\overline{T}}_{\rm m}$ in Eq. (15) where $\overline{\overline{Y}}_{\mathrm{g}}$ is replaced by a diagonal matrix $\alpha\overline{\overline{I}}$ (here $\alpha\to\infty$).

For the operation of the temporal circulator, we define cost functions based on the desired reflection behavior across the three temporal harmonics. Each incident wave is assumed to have a $1~\text{V/m}$ electric-field amplitude. When the metasurface stack is excited by one harmonic (e.g., $\omega_{-1}$), the reflected wave should appear entirely at one of the other two considered harmonics (say, $\omega_{0}$) with the same $1~\text{V/m}$ amplitude, while the remaining harmonic ($\omega_{+1}$) is completely suppressed. This property ensures that waves at each input frequency are uniquely routed to another frequency port without amplitude loss. By cyclically permuting the input among the three harmonics, the same nonreciprocal mapping must hold. The cost functions that penalize deviations from this ideal circulation (unit-amplitude transfer to the designated port and zero leakage to the undesired port) are chosen as

[TABLE]

The total error function of the circulation optimization problem is defined as a summation of all these three functions, that is, $\Delta_{\rm circ}=\Delta_{1}+\Delta_{2}+\Delta_{3}$. The periodicity $d$ and radii $r_{0}$ are chosen to be the same as in the previous example.

Figures 2(b)–(d) show the reflection coefficients of waves at different frequency harmonics for the optimized metasurface stack for three different illumination scenarios (at $\omega_{0}$ and $\omega_{\pm 1}$). The optimized material parameters are listed in the Supplementary Material. Blue arrows in insets of the plots show the excitation frequency, while the light red arrows stand for the optimized reflected frequency harmonics. As we can see from Fig. 2(b), the incident $\omega_{-1}$ harmonic is reflected towards the $\omega_{0}$ temporal port with 1 V/m field amplitude, while the reflected $\omega_{+1}$ harmonic is suppressed to 0.3291 V/m. On the other hand, incident $\omega_{0}$ harmonic is reflected towards the $\omega_{+1}$ temporal port with 0.9991 V/m field amplitude, as is seen in Fig. 2(c) (note the large scale of the vertical axis in the figure). The reflected $\omega_{-1}$ harmonic reaches 0.1895 V/m. Finally, incident harmonic $\omega_{+1}$ is reflected predominantly to the $\omega_{-1}$ temporal port (0.9919 V/m) and minimally coupled to the $\omega_{0}$ port (0.001 V/m), as shown in Fig. 2(d). The value of the final cost function introduced as the sum of the error functions in Eq. (24) after optimization is $0.1739$. The isolation levels for the three incident scenarios mentioned above are $-9.90\rm{dB}$, $-14.90\rm{dB}$, and $-60\rm{dB}$, listed in the same order. It should be mentioned that the optimized metasurface stack is not matched and exhibits also non-zero “retro” reflections, that is, reflections to the same temporal port from which the incident wave arrived. These could be minimized by incorporating additional conditions in the cost function and increasing the number of metasurfaces in the stack. Non-zero reflected sideband harmonics beyond $\omega_{0}$ and $\omega_{\pm 1}$ can be suppressed in practical settings with bandpass frequency filters.

4 Proposed solution for waveguide realization

In this section, we introduce an efficient and straightforward solution for the microwave realization of the proposed stacked metasurfaces. Without loss of generality, in our illustration, we consider a stack consisting of two metasurfaces, as shown in Fig. 3(a). However, the method is generic and applies to an arbitrary number of metasurfaces. A linearly-polarized incident wave, propagating in the normal direction to the metasurfaces plane (with the wave vector along the $y$-axis), has its electric field component aligned parallel to the $z$-axis, while the magnetic field is perpendicular to the cylinders ($x$-axis). Naturally, the realization of the metasurface shown in Fig. 3(a) is very cumbersome as it requires many time-modulated cylindrical meta-atoms with lengths much larger than the operational wavelength.

To implement such structures in practice, one can consider an equivalent system inside a parallel-plate waveguide (PPW) and exploit the image theory principle [32]. Let us consider a PPW configuration shown in Fig. 3(b). Here, to represent each metasurface, we need only one cylindrical meta-atom with the same radius $r_{0}$ and a finite (sub-wavelength) height $h$ being placed between two conducting plates. The width of these plates in the $x$-direction must match the periodicity of the array, that is, $d$. The central region of the structure functions as a PPW supporting the TEM mode. Importantly, in the absence of side walls, the magnetic field at the PPW opening is oriented along the $x$-direction (with vanishing $y$-component), if the $h\ll d\ll\lambda_{0}$ [33]. Therefore, such a PPW configuration effectively emulates perfect electric and magnetic conductor walls separated by $h$ in the $z$-direction and by $d$ in the $x$-direction, respectively. This feature allows us to emulate the infinite geometry in Fig. 3(a) by a finite geometry shown in Fig. 3(b). The cascade of two metasurfaces is modeled by a cascade of two meta-atoms separated by a distance $D$. The height of the PPW can be chosen such that the waveguide characteristic impedance becomes matched to that of standard coaxial connectors, that is, $50\,\Omega$. For ensuring perfect mode conversion between the PPW and the coaxial connector (supporting the dominant cylindrical TEM mode), a tapering section can be designed, as shown in Fig. 3(b). The height and width of the tapering opening section should be configured to maintain constant impedance throughout the entire transition region. The equivalence between the scattering parameters of the free-space system in Fig. 3(a) and the PPW system in Fig. 3(b) is demonstrated in the Supplementary Material.

5 Conclusion

In this paper, a general theoretical framework for designing stacked time-varying metasurfaces is presented, along with several demonstrators of nonreciprocal devices. The plasma and collision frequencies are selected as the modulated material parameters, with material dispersion taken into account. Accurate circuit models are proposed to interpret the behavior of both individual cylindrical meta-atoms and their arrays. Establishing the equivalent circuit model enables straightforward computation of the transfer matrices for time-varying metasurfaces within a cascade. Finally, a practical realization method is introduced using a parallel-plate waveguide configuration. The proposed stacked metasurfaces offer a compact, versatile, and scalable route toward spatiotemporal platforms for next-generation photonic and microwave systems. The developed methodology is not limited to arrays of cylinders and can be used to cascades of any metasurfaces modeled by sheet impedances that are periodically modulated in time.

Further generalization of this approach can be achieved by incorporating spatial non-uniformity in the in-plane metasurface direction. In this respect, it is important to stress the fundamental advantages of the proposed time-modulated nonreciprocal metasurface structures compared to earlier developed solutions based on fully reflective metasurfaces (for example, as in [9]). It is well established that in order to fully control anomalous reflections into desired directions, it is necessary to carefully engineer distributions of high-order evanescent Floquet harmonics [34, 35]. This requirement, applicable also to nonreciprocal wave transformations, demands the use of dense subwavelength arrays of different time-modulated meta-atoms, which are difficult to optimize and realize. The fundamental reason for these difficulties is that the optimal anomalous reflectors must exhibit a non-local response. In contrast, full control of reflection-free transmission to any angle is possible to achieve using locally-responding arrays [34]. This means that the developed design method is also directly applicable to spatially non-uniform time-modulated arrays for transmission control. Since it is enough to optimize only the local response, each meta-atom can be optimized separately in an infinite array of identical meta-atoms, as is done in this paper for uniform arrays. Thus, we expect that this method will open a way to realize nonreciprocal metasurfaces that will be able to control also the direction of propagation of transmitted waves.

Acknowledgements

MM, ST, and VA acknowledge the Horizon 2020 Framework Programme (PULSE). ST and VA acknowledge the Research Council of Finland within the RCF-DoD Future Information Architecture for IoT initiative under grant no. 365679. X.W. acknowledges the Fundamental Research Funds for the Central Universities, China (project no. 3072024WD2603)