Linear-frequency conversion with time-varying metasurfaces

Claude Amra, Ali Passian, Philippe Tchamitchian, Mauro Ettorre, Ahmed Alwakil, Juan Antonio Zapien, Paul Rouquette, Yannick Abautret, Myriam Zerrad

<figure><img src="image_1.png" /></figure>

To cite this version:

Claude Amra, Ali Passian, Philippe Tchamitchian, Mauro Ettorre, Ahmed Alwakil, et al. Linear-frequency conversion with time-varying metasurfaces. Physical Review Research, 2024, 6, pp.013002. 10.1103/PhysRevResearch.6.013002. insu-04479011

HAL Id: insu-04479011

https://insu.hal.science/insu-04479011

Submitted on 20 Apr 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
Linear-frequency conversion with time-varying metasurfaces

Claude Amra*

Aix Marseille Univ, CNRS, Centrale Marseille, Institut Fresnel, Faculté des Sciences - Campus Saint Jérôme,
Avenue Escadrille Normandie-Niémen, 13397 Marseille, France

Ali Passian†

Quantum Computing and Sensing Group, Computational Sciences and Engineering Division,
Oak Ridge National Laboratory, Oak Ridge Tennessee 37831, USA

Philippe Tchamitchian

Independent Mathematician, 33580 Neuffons, France

Mauro Ettorre‡ and Ahmed Alwakil §

Université de Rennes 1, Institut d'Electronique et des Technologies du numéRique (IETR), UMR 6164, Rennes, France

Juan Antonio Zapien||

Department of Materials Science and Engineering, City University of Hong Kong, Kowloon, Hong Kong SAR

Paul Rouquette¶

Aix Marseille Univ, CNRS, CNES, LAM UMR 7326, Marseille, France

Yannick Abautret# and Myriam Zerrad**

Aix Marseille Univ, CNRS, Centrale Marseille, Institut Fresnel, Marseille, France

<figure><img src="image_1.png" /></figure>

(Received 20 September 2023; accepted 21 November 2023; published 2 January 2024)

Frequency conversion is a hallmark of nonlinearity. The spectral manifestations, emergent within a system, can typically be attributed to a marked nonlinearity within the material properties, complex geometric configurations, and/or the unique functional form of interactions taking place in the constitutive subsystems. These phenomena, irrespective of their origins, have been harnessed and exploited in applications ranging from the generation of entangled photons, a cornerstone in quantum technologies, to nanomechanical frequency mixing, advancing subsurface scanning probe microscopy. Here, we propose a frequency conversion mechanism based on time-varying metasurfaces, an emerging frontier in metamaterial research. We show how temporal properties of metasurfaces can effectively emulate a nonlinear medium, thereby facilitating frequency conversion. The proposed material configuration has the potential not only to advance integrated photonics and quantum optics, but also to create opportunities in quantum sensing, quantum materials, and crucially quantum communications.

DOI: 10.1103/PhysRevResearch.6.013002

I. INTRODUCTION

Metamaterials, encompassing metasurfaces, are herald-
ing a new era in functional materials, exhibiting complex
and tailored behaviors [1]. Time-varying metasurfaces, as a
subclass of four-dimensional metamaterials [2] are gaining
considerable prominence in integrated photonics and quan-
tum devices [3]. Exploring metamaterials holds promise not
only for ushering in information processing, for example
via manipulating optical fields [4] but also as use cases for
emerging quantum computing optimization algorithms [5].
Frequency conversion typically harnessed through nonlinear
media [6], is an important instance, where single photons
participating in forming hybrid light-matter networks aid
in connecting shorter wavelengths employed in process-
ing and storage in qubit systems with sources and media
at telecommunications bands [7]. Frequency manipulation,

Published by the American Physical Society under the terms of the
Creative Commons Attribution 4.0 International license. Further
distribution of this work must maintain attribution to the author(s)
and the published article's title, journal citation, and DOI.

<footnote>* claude.amra@fresnel.fr</footnote>

<footnote>† passianan@ornl.gov</footnote>

<footnote>‡ mauro.ettorre@univ-rennes1.fr</footnote>

<footnote>§ ahmed.diaa.alwakil@gmail.com</footnote>

<footnote>|| apjazs@cityu.edu.hk</footnote>

<footnote>¶ paul.rouquette@lam.fr</footnote>

<footnote>\# abautret@fresnel.fr</footnote>

<footnote>**myriam.zerrad@fresnel.fr</footnote>
Linear-frequency conversion with time-varying metasurfaces

Claude Amra*

Aix Marseille Univ, CNRS, Centrale Marseille, Institut Fresnel, Faculté des Sciences - Campus Saint Jérôme,
Avenue Escadrille Normandie-Niémen, 13397 Marseille, France

Ali Passian†

Quantum Computing and Sensing Group, Computational Sciences and Engineering Division,
Oak Ridge National Laboratory, Oak Ridge Tennessee 37831, USA

Philippe Tchamitchian

Independent Mathematician, 33580 Neuffons, France

Mauro Ettorre‡ and Ahmed Alwakil §

Université de Rennes 1, Institut d'Electronique et des Technologies du numéRique (IETR), UMR 6164, Rennes, France

Juan Antonio Zapien||

Department of Materials Science and Engineering, City University of Hong Kong, Kowloon, Hong Kong SAR

Paul Rouquette¶

Aix Marseille Univ, CNRS, CNES, LAM UMR 7326, Marseille, France

Yannick Abautret# and Myriam Zerrad**

Aix Marseille Univ, CNRS, Centrale Marseille, Institut Fresnel, Marseille, France

<figure><img src="image_3.png" /></figure>

(Received 20 September 2023; accepted 21 November 2023; published 2 January 2024)

Frequency conversion is a hallmark of nonlinearity. The spectral manifestations, emergent within a system, can typically be attributed to a marked nonlinearity within the material properties, complex geometric configurations, and/or the unique functional form of interactions taking place in the constitutive subsystems. These phenomena, irrespective of their origins, have been harnessed and exploited in applications ranging from the generation of entangled photons, a cornerstone in quantum technologies, to nanomechanical frequency mixing, advancing subsurface scanning probe microscopy. Here, we propose a frequency conversion mechanism based on time-varying metasurfaces, an emerging frontier in metamaterial research. We show how temporal properties of metasurfaces can effectively emulate a nonlinear medium, thereby facilitating frequency conversion. The proposed material configuration has the potential not only to advance integrated photonics and quantum optics, but also to create opportunities in quantum sensing, quantum materials, and crucially quantum communications.

DOI: 10.1103/PhysRevResearch.6.013002

# I. INTRODUCTION

Metamaterials, encompassing metasurfaces, are heralding a new era in functional materials, exhibiting complex and tailored behaviors [1]. Time-varying metasurfaces, as a subclass of four-dimensional metamaterials [2] are gaining considerable prominence in integrated photonics and quantum devices [3]. Exploring metamaterials holds promise not only for ushering in information processing, for example via manipulating optical fields [4] but also as use cases for emerging quantum computing optimization algorithms [5]. Frequency conversion typically harnessed through nonlinear media [6], is an important instance, where single photons participating in forming hybrid light-matter networks aid in connecting shorter wavelengths employed in processing and storage in qubit systems with sources and media at telecommunications bands [7]. Frequency manipulation,

Published by the American Physical Society under the terms of the Creative Commons Attribution 4.0 International license. Further distribution of this work must maintain attribution to the author(s) and the published article's title, journal citation, and DOI.

<footnote>*claude.amra@fresnel.fr</footnote>

<footnote>†passianan@ornl.gov</footnote>

<footnote>‡mauro.ettorre@univ-rennes1.fr</footnote>

<footnote>§ahmed.diaa.alwakil@gmail.com</footnote>

<footnote>||apjazs@cityu.edu.hk</footnote>

<footnote>¶paul.rouquette@lam.fr</footnote>

<footnote>#abautret@fresnel.fr</footnote>

<footnote>**myriam.zerrad@fresnel.fr</footnote>
without necessarily invoking nonlinear materials, may be exploited, e.g., via nonlinear nanomechanical processes to achieve subsurface nanometrology [8,9]. In the present work on the interaction of light with the specific class of materials of temporal metasurfaces, to be more precise when reflecting the underlying physical processes, it seems prudent to consider a brief note on the spectral components generated when a pump frequency is introduced. The traditional understanding of nonlinear optical (or others, e.g., mechanical and optomechanical [9]) interactions, such as those observed in a Kerr medium [10], invokes the concept of harmonics or the fundamental frequency's integer multiples, appearing as distinct and widely spaced spectral lines in relation to the primary frequency. The orderly progression from the fundamental to the second harmonic, and so on, is a testament to the material's inherent nonlinearity and its response to the driving (e.g., electric, stress) field. Contrasting this with the time-varying metasurfaces subject to a pump frequency leading to frequency mixing that generates sidebands, which can be symmetrically spaced about the pump frequency but are not integer multiples of it. The proximity of these sidebands to the pump frequency is typically predicated on the modulation frequency being considerably smaller than the pump frequency. This spectral proximity may be seen as a distinctive signature of temporal metasurfaces, differentiating them from the well-separated harmonics due to traditional nonlinear interactions. However, as shown here, when the modulation frequency is of a magnitude that is not insignificant in comparison to the pump frequency, the resultant sidebands are no longer subtle shifts but are more spread out, potentially mimicking the spacing one would associate with traditional harmonics. Yet, despite this resemblance, the generation mechanism remains distinctly different. Unlike harmonics, which are a natural consequence of nonlinearity and are strictly tied as integer multiples to the original frequency, sidebands are a result of an external and controllable imposed modulation. Naturally, the distinction between harmonics and sidebands is not only about their spectral position but also the underlying mechanism. While harmonics in nonlinear optics originate directly from the intrinsic properties of the material, sidebands emerge via deliberate, time-dependent manipulation of the metasurface, that is, control versus intrinsic property, or tunability versus predetermined outcome. It is important to note that our approach, as delineated in this work, imposes no limitations on the frequency separation between the converted frequencies and the original incoming frequency. This lack of restriction highlights the versatility and broad applicability of our proposed methodology.

Metasurface advancements, paving the way for innovative manipulations of light and waves through time-modulated metasurfaces, promise revolutionary applications, particularly in secure communications, nonreciprocal devices, and dynamic control over wavefronts. Examples include the work of Taravati and Eleftheriades [11] who introduced a concept of linear-frequency conversion, providing new pathways for manipulating fields, and the work of Guo et al. [12] who explored nonreciprocal metasurfaces through space-time phase modulation, adding a new dimension to control principles. Similarly, the work of Shaltout et al. [13] exploring spatiotemporal control, offer intricate light manipulation strategies,

paralleled by Shaltout, Kildishev, and Shalaev's insights [14] into Lorentz nonreciprocity, highlighting temporal variations in metasurfaces. Sedeh, Salary, and Mosallaei [15] further this exploration by enabling time-varying optical vortices through time-modulated metasurfaces, with Sedeh et al. [16] providing new perspectives on topological space-time transitions. Salary and Mosallaei's work [17] integrates conducting oxide metasurfaces with time modulation for adaptive communications, a theme advanced by Wu and Grbic [18] through their work on serrodyne frequency translation. Sedeh et al. [19] also propose a novel approach for secure communication via graphene-based time-modulated metasurfaces. Conceptualizing spacetime metamaterials, Caloz and Deck-Léger [20] redefine perceptions of material properties in electromagnetic theory.

## II. INHOMOGENEOUS TEMPORAL CONVOLUTION

In the linear regime, the response of a system to an optical field $u_{in}$ is governed by the frequency $\omega$-dependent susceptibility, $\chi(\omega)$, of the system. Outside of the linear regime, frequency mixing can occur, and the system response becomes a function of the order of nonlinearity. As a prelude to our theory, it is instructive to consider, in the time $t$ domain, a perturbative expansion of the response function in terms of different orders of susceptibilities: for an input (scalar) field $u_{in}(t)$, the output field, in the absence of a permanent output, may be expressed as:

$$ u_{\text{out}}(t) = \chi^{(1)}(t) * u_{\text{in}}(t) + \chi^{(2)}(t, t') * [u_{\text{in}}(t) * u_{\text{in}}(t')] \\ + \chi^{(3)}(t, t', t'') * [u_{\text{in}}(t) * u_{\text{in}}(t') * u_{\text{in}}(t'')] \dots, \quad (1) $$

where $\chi^{(n)}$ is the nth-order susceptibility of the medium (alternatively one may write $u_{out}$ in the integral form, where the causality variables $t', t'', \dots$ are integrated over). The higher-order terms are known to lead to harmonic generation, four-wave mixing, spontaneous parametric down conversion, and other nonlinear phenomena. We may consider the convolution $\chi^{(1)}(t) * u_{in}(t)$ in (1), which characterizes the linear response of the system, as depicted in Fig. 1, such that in Fourier domain $u(t) \xrightarrow{\mathcal{F}} U(\omega)$, the output reads $U_{\text{out}}(\omega) = \chi(\omega)U_{\text{in}}(\omega)$. In this work, we will exploit the extension of this linear response to the more general form $\chi(t, t')$, that is, a susceptibility that accounts for externally induced memory effects.

Let $u = e, h$ select the electric $\mathbf{X}_e = \mathbf{E}$, or magnetic $\mathbf{X}_h = \mathbf{H}$ field. Then, for the response $\mathbf{P}_u$ to an input field $\mathbf{X}_u$, we write:

$$ \mathbf{P}_u = \chi_u * \mathbf{X}_u = \int \chi_u(t, t') \mathbf{X}_u(t - t') dt', \quad (2) $$

taking into account possible delayed response of the material to the field, and potentially noninstantaneous effects, if $\chi_u(t, t')$ is not a $\delta$ function of $t'$. Thus, the specificity of active metasurfaces (having physical properties that vary with time) lies in the inhomogeneous temporal convolution (2), which involves an additional dependence of the susceptibility versus the current time variable $t$.
without necessarily invoking nonlinear materials, may be exploited, e.g., via nonlinear nanomechanical processes to achieve subsurface nanometrology [8,9]. In the present work on the interaction of light with the specific class of materials of temporal metasurfaces, to be more precise when reflecting the underlying physical processes, it seems prudent to consider a brief note on the spectral components generated when a pump frequency is introduced. The traditional understanding of nonlinear optical (or others, e.g., mechanical and optomechanical [9]) interactions, such as those observed in a Kerr medium [10], invokes the concept of harmonics or the fundamental frequency's integer multiples, appearing as distinct and widely spaced spectral lines in relation to the primary frequency. The orderly progression from the fundamental to the second harmonic, and so on, is a testament to the material's inherent nonlinearity and its response to the driving (e.g., electric, stress) field. Contrasting this with the time-varying metasurfaces subject to a pump frequency leading to frequency mixing that generates sidebands, which can be symmetrically spaced about the pump frequency but are not integer multiples of it. The proximity of these sidebands to the pump frequency is typically predicated on the modulation frequency being considerably smaller than the pump frequency. This spectral proximity may be seen as a distinctive signature of temporal metasurfaces, differentiating them from the well-separated harmonics due to traditional nonlinear interactions. However, as shown here, when the modulation frequency is of a magnitude that is not insignificant in comparison to the pump frequency, the resultant sidebands are no longer subtle shifts but are more spread out, potentially mimicking the spacing one would associate with traditional harmonics. Yet, despite this resemblance, the generation mechanism remains distinctly different. Unlike harmonics, which are a natural consequence of nonlinearity and are strictly tied as integer multiples to the original frequency, sidebands are a result of an external and controllable imposed modulation. Naturally, the distinction between harmonics and sidebands is not only about their spectral position but also the underlying mechanism. While harmonics in nonlinear optics originate directly from the intrinsic properties of the material, sidebands emerge via deliberate, time-dependent manipulation of the metasurface, that is, control versus intrinsic property, or tunability versus predetermined outcome. It is important to note that our approach, as delineated in this work, imposes no limitations on the frequency separation between the converted frequencies and the original incoming frequency. This lack of restriction highlights the versatility and broad applicability of our proposed methodology.

Metasurface advancements, paving the way for innovative manipulations of light and waves through time-modulated metasurfaces, promise revolutionary applications, particularly in secure communications, nonreciprocal devices, and dynamic control over wavefronts. Examples include the work of Taravati and Eleftheriades [11] who introduced a concept of linear-frequency conversion, providing new pathways for manipulating fields, and the work of Guo et al. [12] who explored nonreciprocal metasurfaces through space-time phase modulation, adding a new dimension to control principles. Similarly, the work of Shaltout et al. [13] exploring spatiotemporal control, offer intricate light manipulation strategies,

paralleled by Shaltout, Kildishev, and Shalaev's insights [14] into Lorentz nonreciprocity, highlighting temporal variations in metasurfaces. Sedeh, Salary, and Mosallaei [15] further this exploration by enabling time-varying optical vortices through time-modulated metasurfaces, with Sedeh et al. [16] providing new perspectives on topological space-time transitions. Salary and Mosallaei's work [17] integrates conducting oxide metasurfaces with time modulation for adaptive communications, a theme advanced by Wu and Grbic [18] through their work on serrodyne frequency translation. Sedeh et al. [19] also propose a novel approach for secure communication via graphene-based time-modulated metasurfaces. Conceptualizing spacetime metamaterials, Caloz and Deck-Léger [20] redefine perceptions of material properties in electromagnetic theory.

## II. INHOMOGENEOUS TEMPORAL CONVOLUTION

In the linear regime, the response of a system to an optical field $u_{in}$ is governed by the frequency $\omega$-dependent susceptibility, $\chi(\omega)$, of the system. Outside of the linear regime, frequency mixing can occur, and the system response becomes a function of the order of nonlinearity. As a prelude to our theory, it is instructive to consider, in the time $t$ domain, a perturbative expansion of the response function in terms of different orders of susceptibilities: for an input (scalar) field $u_{in}(t)$, the output field, in the absence of a permanent output, may be expressed as:

$$ u_{\text{out}}(t) = \chi^{(1)}(t) * u_{\text{in}}(t) + \chi^{(2)}(t, t') * [u_{\text{in}}(t) * u_{\text{in}}(t')] \\ + \chi^{(3)}(t, t', t'') * [u_{\text{in}}(t) * u_{\text{in}}(t') * u_{\text{in}}(t'')] \dots, \quad (1) $$

where $\chi^{(n)}$ is the nth-order susceptibility of the medium (alternatively one may write $u_{out}$ in the integral form, where the causality variables $t', t'', \dots$ are integrated over). The higher-order terms are known to lead to harmonic generation, four-wave mixing, spontaneous parametric down conversion, and other nonlinear phenomena. We may consider the convolution $\chi^{(1)}(t) * u_{in}(t)$ in (1), which characterizes the linear response of the system, as depicted in Fig. 1, such that in Fourier domain $u(t) \xrightarrow{\mathcal{F}} U(\omega)$, the output reads $U_{\text{out}}(\omega) = \chi(\omega)U_{\text{in}}(\omega)$. In this work, we will exploit the extension of this linear response to the more general form $\chi(t, t')$, that is, a susceptibility that accounts for externally induced memory effects.

Let $u = e, h$ select the electric $\mathbf{X}_e = \mathbf{E}$, or magnetic $\mathbf{X}_h = \mathbf{H}$ field. Then, for the response $\mathbf{P}_u$ to an input field $\mathbf{X}_u$, we write:

$$ \mathbf{P}_u = \chi_u * \mathbf{X}_u = \int \chi_u(t, t') \mathbf{X}_u(t - t') dt', \quad (2) $$

taking into account possible delayed response of the material to the field, and potentially noninstantaneous effects, if $\chi_u(t, t')$ is not a $\delta$ function of $t'$. Thus, the specificity of active metasurfaces (having physical properties that vary with time) lies in the inhomogeneous temporal convolution (2), which involves an additional dependence of the susceptibility versus the current time variable $t$.
<figure><img src="image_2.png" /><figcaption>FIG. 1. Linear-frequency conversion via a metasurface with memory. The active metasurface ($z = 0$, thickness $\delta z$) generates a continuum of frequencies $\{\omega\}$ when illuminated by a harmonic of pulsation $\omega_1$. The blue arrow indicates the interaction $e(t)$ with the external environment, which changes the susceptibility $\chi$ of the metasurface at high frequency over time. Annotated are the indices $n$ of the media, excitation, and scattered fields $E$ and their transforms, and (in white) the reduction of the continuum to a single frequency $\omega_1$ in reflection, and a single frequency $\omega_2$ in transmission.</figcaption></figure>

Prior to setting up the metasurface equations, a brief clarification of the necessity and implications of $\chi(t) \rightarrow \chi(t, t')$ transition, that is, introducing a temporally generalized susceptibility, seems logical. We first note that the susceptibility $\chi(t')$ represents the metamaterial's response to an incoming field, accounting for memory effects linked to the inertia of the material's charges and the retardation effects due to the finite time for the polarization to respond to field changes. Here, inertia refers to the resistance of the charge carriers to (motional) changes, e.g., the delay before electrons reach a new equilibrium position in response to the electric field. This delay or memory of past interactions is embedded in the susceptibility $\chi(t')$ and is what imparts the material with dispersive and absorptive properties, observable in the form of phase shifts and attenuation in the transmitted light. However, when considering a metasurface under an additional time-varying external stimulus, such as electrical, optical, mechanical, or thermal excitations, the system's response at any given time $t$ becomes influenced by the external control exerted up to that point. This is not merely a retardation effect; it is a dynamic modulation of the metasurface's properties over time. Therefore, we expand the conventional susceptibility to $\chi(t, t')$ to account for these additional memory effects, where $\chi(t, t')$ signifies that the response of the material at time $t$ now depends on both the historical electric field and the time-evolving state of the material due to the external excitation. This expanded formalism $\chi(t) \rightarrow \chi(t, t')$ is crucial for capturing the response of the metasurface given its properties are no longer static but evolving, leading to a time-variant interaction with the incident field. This evolution introduces a form of temporal inhomogeneity, akin to spatial inhomogeneity in traditional metasurfaces. Just as spatially varying features can cause the incident light to experience different local environments, the temporally varying properties mean that light at different times $\{t\}$ encounters effectively

different materials. Thus, our work draws a distinction between the classical memory effects encapsulated by $\chi(t')$ (characterized by dispersion and absorption due to inertia and retardation) and the additional time-dependent modulations described by $\chi(t, t')$. The adoption of $\chi(t, t')$ in our work is not merely a formal generalization but a necessary framework to model the interactions that define the adaptive, tunable metasurfaces.

Our objective is to explore the conversion process $\omega \xrightarrow{\chi} \omega'$ for photons transmitting through a metasurface characterized by $\chi(t, t')$ and to elucidate its spectral versatility. Let us consider a metasurface positioned at $z = 0$ (Fig. 1) at the separation of two nondispersive transparent media, with real refractive indices $n_1$ and $n_2$. The thickness of the metasurface ($\delta z$ in Fig. 1), is small compared to $\lambda$, the wavelengths of the radiation. The physical properties of this metasurface are assumed to vary at high frequency versus time $t$, thanks to the interaction with an external channel $e(t)$. The component is illuminated by collimated, monochromatic radiation at frequency $\omega_1$. A continuum of frequencies $\{\omega\}$ is then generated by reflection and transmission, in the specular directions. Let $P_e$ and $P_h$ denote the electric and magnetic polarizations of the metasurface. Then, Maxwell's equations (classical, macroscopic) can be written in the space-time domain and in the sense of distributions, as:

$$rotE = -\frac{\partial B}{\partial t} + M\delta(z), \quad (3)$$

$$rotH = \frac{\partial D}{\partial t} + J\delta(z), \quad (4)$$

with:

$$M = -\mu_0 \frac{\partial P_h}{\partial t}, \quad P_h = \chi_h * H, \quad (5)$$

$$J = \frac{\partial P_e}{\partial t}, \quad P_e = \varepsilon_0 \chi_e * E, \quad (6)$$

where, $M$ and $J$ denote magnetic and electric currents proportional to the time derivatives of the $P_h$ and $P_e$ polarizations of the material. In the following for notational clarity, we adhere to the following Fourier transform definitions:

$$\begin{align}
\tilde{f}(\omega, t') &= \mathcal{F}_t\{f(t, t')\}(\omega) = \left(\frac{1}{2\pi}\right) \int_{-\infty}^{\infty} f(t, t')e^{i\omega t'} dt, \nonumber \\
\tilde{f}(t, \omega') &= \mathcal{F}_{t'}\{f(t, t')\}(\omega') = \left(\frac{1}{2\pi}\right) \int_{-\infty}^{\infty} f(t, t')e^{i\omega t'} dt', \nonumber \\
\hat{f}(\omega, \omega') &= \mathcal{F}_{t,t'}\{f(t, t')\}(\omega, \omega') \nonumber \\
&= \left(\frac{1}{4\pi^2}\right) \int_{-\infty}^{\infty} \int_{-\infty}^{\infty} f(t, t')e^{i\omega t'} e^{i\omega t'} dt dt'. \tag{7}
\end{align}$$

The specificity of the so-called active metasurface (whose physical properties vary with time) lies in the inhomogeneous temporal convolution noted *, namely:

$$\begin{align}
\chi_u * X_u &= \int \chi_u(t, t') X_u(t - t') dt', \nonumber \\
u &= e, h, \text{ and } X = E, H, \tag{8}
\end{align}$$

where the susceptibility $\chi_u(t, t')$ is a function of two time variables $t$ and $t'$. The electromagnetic field $X$ is calculated at time
<figure><img src="image_4.png" /><figcaption>FIG. 1. Linear-frequency conversion via a metasurface with memory. The active metasurface ($z = 0$, thickness $\delta z$) generates a continuum of frequencies $\{\omega\}$ when illuminated by a harmonic of pulsation $\omega_1$. The blue arrow indicates the interaction $e(t)$ with the external environment, which changes the susceptibility $\chi$ of the metasurface at high frequency over time. Annotated are the indices $n$ of the media, excitation, and scattered fields $E$ and their transforms, and (in white) the reduction of the continuum to a single frequency $\omega_1$ in reflection, and a single frequency $\omega_2$ in transmission.</figcaption></figure>

Prior to setting up the metasurface equations, a brief clarification of the necessity and implications of $\chi(t) \to \chi(t, t')$ transition, that is, introducing a temporally generalized susceptibility, seems logical. We first note that the susceptibility $\chi(t')$ represents the metamaterial's response to an incoming field, accounting for memory effects linked to the inertia of the material's charges and the retardation effects due to the finite time for the polarization to respond to field changes. Here, inertia refers to the resistance of the charge carriers to (motional) changes, e.g., the delay before electrons reach a new equilibrium position in response to the electric field. This delay or memory of past interactions is embedded in the susceptibility $\chi(t')$ and is what imparts the material with dispersive and absorptive properties, observable in the form of phase shifts and attenuation in the transmitted light. However, when considering a metasurface under an additional time-varying external stimulus, such as electrical, optical, mechanical, or thermal excitations, the system's response at any given time $t$ becomes influenced by the external control exerted up to that point. This is not merely a retardation effect; it is a dynamic modulation of the metasurface's properties over time. Therefore, we expand the conventional susceptibility to $\chi(t, t')$ to account for these additional memory effects, where $\chi(t, t')$ signifies that the response of the material at time $t$ now depends on both the historical electric field and the time-evolving state of the material due to the external excitation. This expanded formalism $\chi(t) \to \chi(t, t')$ is crucial for capturing the response of the metasurface given its properties are no longer static but evolving, leading to a time-variant interaction with the incident field. This evolution introduces a form of temporal inhomogeneity, akin to spatial inhomogeneity in traditional metasurfaces. Just as spatially varying features can cause the incident light to experience different local environments, the temporally varying properties mean that light at different times $\{t\}$ encounters effectively

different materials. Thus, our work draws a distinction between the classical memory effects encapsulated by $\chi(t')$ (characterized by dispersion and absorption due to inertia and retardation) and the additional time-dependent modulations described by $\chi(t, t')$. The adoption of $\chi(t, t')$ in our work is not merely a formal generalization but a necessary framework to model the interactions that define the adaptive, tunable metasurfaces.

Our objective is to explore the conversion process $\omega \xrightarrow{\chi} \omega'$ for photons transmitting through a metasurface characterized by $\chi(t, t')$ and to elucidate its spectral versatility. Let us consider a metasurface positioned at $z = 0$ (Fig. 1) at the separation of two nondispersive transparent media, with real refractive indices $n_1$ and $n_2$. The thickness of the metasurface ($\delta z$ in Fig. 1), is small compared to $\lambda$, the wavelengths of the radiation. The physical properties of this metasurface are assumed to vary at high frequency versus time $t$, thanks to the interaction with an external channel $e(t)$. The component is illuminated by collimated, monochromatic radiation at frequency $\omega_1$. A continuum of frequencies $\{\omega\}$ is then generated by reflection and transmission, in the specular directions. Let $P_e$ and $P_h$ denote the electric and magnetic polarizations of the metasurface. Then, Maxwell's equations (classical, macroscopic) can be written in the space-time domain and in the sense of distributions, as:

$$rotE = -\frac{\partial B}{\partial t} + M\delta(z), \quad (3)$$

$$rotH = \frac{\partial D}{\partial t} + J\delta(z), \quad (4)$$

with:

$$M = -\mu_0 \frac{\partial P_h}{\partial t}, \quad P_h = \chi_h * H, \quad (5)$$

$$J = \frac{\partial P_e}{\partial t}, \quad P_e = \varepsilon_0 \chi_e * E, \quad (6)$$

where, $M$ and $J$ denote magnetic and electric currents proportional to the time derivatives of the $P_h$ and $P_e$ polarizations of the material. In the following for notational clarity, we adhere to the following Fourier transform definitions:

$$ \begin{aligned} \tilde{f}(\omega, t') &= \mathcal{F}_t\{f(t, t')\}(\omega) = \left(\frac{1}{2\pi}\right) \int_{-\infty}^{\infty} f(t, t')e^{i\omega t'} dt, \\ \tilde{f}(t, \omega') &= \mathcal{F}_{t'}\{f(t, t')\}(\omega') = \left(\frac{1}{2\pi}\right) \int_{-\infty}^{\infty} f(t, t')e^{i\omega t'} dt', \\ \hat{f}(\omega, \omega') &= \mathcal{F}_{t,t'}\{f(t, t')\}(\omega, \omega') \\ &= \left(\frac{1}{4\pi^2}\right) \int_{-\infty}^{\infty} \int_{-\infty}^{\infty} f(t, t')e^{i\omega t'} e^{i\omega t'} dt dt'. \end{aligned} \quad (7) $$

The specificity of the so-called active metasurface (whose physical properties vary with time) lies in the inhomogeneous temporal convolution noted *, namely:

$$ \begin{aligned} \chi_u * X_u &= \int \chi_u(t, t') X_u(t - t') dt', \\ u &= e, h, \text{ and } X = E, H, \end{aligned} \quad (8) $$

where the susceptibility $\chi_u(t, t')$ is a function of two time variables $t$ and $t'$. The electromagnetic field $X$ is calculated at time
$t$, while the variable $t'$ (the causality variable) allows inertial phenomena to be taken into account. The (causality) variable $t'$ is understood to account for intrinsic inertial phenomena [$\chi_u(t, t') = 0, \forall t' < 0$].

We note here, by comparison with the classical case of so-called passive surfaces, that time $t$ intervenes in the susceptibility. This means that the weighting of inertia in Eq. (8) depends not only on the delay $t'$, but also on the time $t$ because of an external action acting in parallel on the material at high frequency. Under these conditions, harmonics of the type $e^{-i\omega_1 t}$ are no longer eigenvalues of the polarization, since we have for example for $u = h$, that is, $X_h = H$:

$$ 
\begin{aligned}
H(t) = e^{-i\omega_1 t} &\Rightarrow P_h(t) = \chi_h * H \\
&= e^{-i\omega_1 t} \int \chi_h(t, t') e^{i\omega_1 t'} dt' = 2\pi e^{-i\omega_1 t} \tilde{\chi}_h(t, \omega_i),
\end{aligned}
\quad (9) $$

where the Fourier transform of $\chi_h(t, t')$ with respect to the variable $t'$ is denoted $\tilde{\chi}_h(t, \omega)$. We thus see in Eq. (9) that the response to a harmonic is no longer a harmonic, given the time dependence of $\tilde{\chi}_h(t, \omega_i)$, which is characteristic of the generated continuum. The classical case of a passive surface can be found by writing $\chi_h(t, t') = \chi_h(t')$, that is, $\tilde{\chi}_h(t, \omega_i) = \tilde{\chi}_h(\omega_i)$ in the Fourier plane, so that the harmonics once again become eigenvalues of the polarization, without continuum generation.

### III. FIELD CALCULATIONS

The aim of this work is to demonstrate that both the electrical and magnetic susceptibilities $\chi_{u=e,h}(t, t')$ of the metasurface may be found such that this continuum {$\omega$} transforms into a single predefined frequency $\omega_2$, that is, achieving conversion in a way similar to that of nonlinear optics. To calculate the fields, let $\tilde{X}$ and $\tilde{X}'$ stand for the Fourier transforms of $X$ with respect to $t$, and $t'$, respectively. We obtain, for the monochromatic incident field in the superstrate ($z < 0$) at frequency $\omega_1$ (omitting the spatial dependence for simplicity):

$$ \frac{E_1^+(t, \omega_1)}{A_1^+(\omega_1)} = e^{-i\omega_1 t} \Rightarrow \frac{\tilde{E}_1^+(\omega, \omega_1)}{A_1^+(\omega_1)} = \delta(\omega - \omega_1), \quad z < 0. \quad (10) $$

As for the reflected field, it carries a priori a frequency continuum in the superstrate in the form:

$$ E_1^-(t, \omega_1) = \int \tilde{E}_1^-(\omega, \omega_1) e^{-i\omega t} d\omega. \quad (11) $$

Now, noting that the transmitted field also carries a continuum {$\omega$}, we introduce a reflection $r$ and transmission $t$ coefficient and write:

$$ \tilde{E}_1^-(\omega, \omega_1) = r(\omega, \omega_1)A_1^+(\omega_1), \quad z < 0, \quad (12) $$

$$ \tilde{E}_2^+(\omega, \omega_1) = t(\omega, \omega_1)A_1^+(\omega_1), \quad z > 0. \quad (13) $$

For simplicity, we will restrict the treatment to the transverse electric polarization mode of the incident field. In the presence of surface currents, the tangential discontinuities $z \times \delta X$ at surface $z = 0$ read:

$$ z \times \delta E = M = -\mu_0 \frac{\partial P_h}{\partial t}, \quad \text{and} \quad z \times \delta H = J = \frac{\partial P_e}{\partial t}, \quad (14) $$

and thus in the Fourier plane:

$$ z \times \delta \tilde{E} = \tilde{M} = i\omega\mu_0\tilde{P}_h, \quad \text{and} \quad z \times \delta \tilde{H} = \tilde{J} = -i\omega\tilde{P}_e. \quad (15) $$

We now need to develop these expressions by taking into account the analytical form of the fields in the surrounding media with indices $n_1$ (the superstrate) and $n_2$ (the substrate). Since these media are assumed to be homogeneous, linear, and isotropic, the fields have been written in the form of a frequency packet, see Eq. (11). We may omit the space variable $\rho = (x, y, z) = (r, z)$ in the argument of the field. In fact the spatial variation of the field is a linear combination of terms of the form

$$ e^{i[\sigma \cdot r \pm \alpha_i(\sigma, \omega)z]}, \quad \text{where} \quad \sigma = k_i \sin \theta_i (1, 0), $$

is the media-independent spatial pulsation, with

$$ \alpha_i = \sqrt{k_i^2 - \sigma^2} = k_i \cos \theta_i, \quad \text{and} \quad k_i = 2\pi n_i / \lambda. $$

The wavelength (in vacuum) is denoted by $\lambda$ and the direction of the wave in the plane of incidence $xz$ is denoted $\theta_i$ in medium $i$. In the following, for simplicity, we will continue omitting this spatial variation, which is justified since the discontinuity will be written specifically at $z = 0$, and the $xy$ variation remains the same in all media. Using the notion of an effective index ($\tilde{n}$) to relate the tangential components of the Fourier transforms of the electric and magnetic fields in a homogeneous medium, i.e., in TE polarization, we establish:

$$ 
\begin{aligned}
\tilde{H}_i^{\pm} &= \pm \tilde{n}_i z \times \tilde{E}_i^{\pm}, \\
\tilde{n}_i &= \frac{\alpha_i}{\omega \mu_i} = \left( \frac{1}{\eta_0 \mu_r^i} \right) n_i \cos \theta_i, \\
\eta_0 &= \sqrt{\frac{\mu_0}{\varepsilon_0}},
\end{aligned}
\quad (16) $$

where $\mu_r^i$ is the relative permeability and $\eta_0$ is the vacuum impedance. The + sign must be taken for progressive waves, and the - sign for retrograde waves. Finally, Eqs. (14) and (15), combined with Eq. (16), give for TE polarization:

$$ \frac{\tilde{M}}{A_1^+(\omega_1)} = \delta(\omega - \omega_1) - t(\omega, \omega_1) + r(\omega, \omega_1), \quad (17) $$

$$ \frac{\tilde{J}}{A_1^+(\omega_1)} = \tilde{n}_1(\omega)[\delta(\omega - \omega_1) - r(\omega, \omega_1)] - \tilde{n}_2(\omega)t(\omega, \omega_1). \quad (18) $$

The next step is to develop the currents $M$ and $J$. By the virtue of the convolution (2), we write: $\tilde{J}(\omega) = -\varepsilon_0 \mathcal{I}_e$, and $\tilde{M}(\omega) = \mu_0 \mathcal{I}_h$, with

$$ \mathcal{I}_u = 2\pi i\omega \int \hat{\chi}_u(\omega', \omega - \omega') \tilde{X}_u(\omega - \omega') d\omega', \quad (19) $$

where $u = e, h$, and $\hat{\chi}_u$ denotes the double Fourier transform of the susceptibility. A question that now arises is what fields to choose in these convolutions that are integrals of discontinuous quantities taken at $z = 0$. To overcome this difficulty, we invoke a suitable mean field $X_{ua} = [X_{u1}(t) + X_{u2}(t)]/2$, and thus obtain:

$$ \frac{2\tilde{E}_a}{A_1^+(\omega_1)} = r(\omega, \omega_1) + t(\omega, \omega_1) + \delta(\omega - \omega_1), \quad (20) $$
$t$, while the variable $t'$ (the causality variable) allows inertial phenomena to be taken into account. The (causality) variable $t'$ is understood to account for intrinsic inertial phenomena [$\chi_u(t, t') = 0, \forall t' < 0$].

We note here, by comparison with the classical case of so-called passive surfaces, that time $t$ intervenes in the susceptibility. This means that the weighting of inertia in Eq. (8) depends not only on the delay $t'$, but also on the time $t$ because of an external action acting in parallel on the material at high frequency. Under these conditions, harmonics of the type $e^{-i\omega_1 t}$ are no longer eigenvalues of the polarization, since we have for example for $u = h$, that is, $X_h = H$:

$$ 
\begin{aligned}
H(t) = e^{-i\omega_1 t} &\Rightarrow P_h(t) = \chi_h * H \\
&= e^{-i\omega_1 t} \int \chi_h(t, t') e^{i\omega_1 t'} dt' = 2\pi e^{-i\omega_1 t} \tilde{\chi}_h(t, \omega_i),
\end{aligned}
\quad (9) $$

where the Fourier transform of $\chi_h(t, t')$ with respect to the variable $t'$ is denoted $\tilde{\chi}_h(t, \omega)$. We thus see in Eq. (9) that the response to a harmonic is no longer a harmonic, given the time dependence of $\tilde{\chi}_h(t, \omega_i)$, which is characteristic of the generated continuum. The classical case of a passive surface can be found by writing $\chi_h(t, t') = \chi_h(t')$, that is, $\tilde{\chi}_h(t, \omega_i) = \tilde{\chi}_h(\omega_i)$ in the Fourier plane, so that the harmonics once again become eigenvalues of the polarization, without continuum generation.

### III. FIELD CALCULATIONS

The aim of this work is to demonstrate that both the electrical and magnetic susceptibilities $\chi_{u=e,h}(t, t')$ of the metasurface may be found such that this continuum $\{\omega\}$ transforms into a single predefined frequency $\omega_2$, that is, achieving conversion in a way similar to that of nonlinear optics. To calculate the fields, let $\tilde{X}$ and $\tilde{X}'$ stand for the Fourier transforms of $X$ with respect to $t$, and $t'$, respectively. We obtain, for the monochromatic incident field in the superstrate ($z < 0$) at frequency $\omega_1$ (omitting the spatial dependence for simplicity):

$$ \frac{E_1^+(t, \omega_1)}{A_1^+(\omega_1)} = e^{-i\omega_1 t} \Rightarrow \frac{\tilde{E}_1^+(\omega, \omega_1)}{A_1^+(\omega_1)} = \delta(\omega - \omega_1), \quad z < 0. \quad (10) $$

As for the reflected field, it carries a priori a frequency continuum in the superstrate in the form:

$$ E_1^-(t, \omega_1) = \int \tilde{E}_1^-(\omega, \omega_1) e^{-i\omega t} d\omega. \quad (11) $$

Now, noting that the transmitted field also carries a continuum $\{\omega\}$, we introduce a reflection $r$ and transmission $t$ coefficient and write:

$$ \tilde{E}_1^-(\omega, \omega_1) = r(\omega, \omega_1)A_1^+(\omega_1), \quad z < 0, \quad (12) $$

$$ \tilde{E}_2^+(\omega, \omega_1) = t(\omega, \omega_1)A_1^+(\omega_1), \quad z > 0. \quad (13) $$

For simplicity, we will restrict the treatment to the transverse electric polarization mode of the incident field. In the presence of surface currents, the tangential discontinuities $z \times \delta X$ at surface $z = 0$ read:

$$ z \times \delta E = M = -\mu_0 \frac{\partial P_h}{\partial t}, \quad \text{and} \quad z \times \delta H = J = \frac{\partial P_e}{\partial t}, \quad (14) $$

and thus in the Fourier plane:

$$ z \times \delta \tilde{E} = \tilde{M} = i\omega\mu_0\tilde{P}_h, \quad \text{and} \quad z \times \delta \tilde{H} = \tilde{J} = -i\omega\tilde{P}_e. \quad (15) $$

We now need to develop these expressions by taking into account the analytical form of the fields in the surrounding media with indices $n_1$ (the superstrate) and $n_2$ (the substrate). Since these media are assumed to be homogeneous, linear, and isotropic, the fields have been written in the form of a frequency packet, see Eq. (11). We may omit the space variable $\rho = (x, y, z) = (r, z)$ in the argument of the field. In fact the spatial variation of the field is a linear combination of terms of the form

$$ e^{i[\sigma \cdot r \pm \alpha_i(\sigma, \omega)z]}, \quad \text{where} \quad \sigma = k_i \sin \theta_i (1, 0), $$

is the media-independent spatial pulsation, with

$$ \alpha_i = \sqrt{k_i^2 - \sigma^2} = k_i \cos \theta_i, \quad \text{and} \quad k_i = 2\pi n_i / \lambda. $$

The wavelength (in vacuum) is denoted by $\lambda$ and the direction of the wave in the plane of incidence $xz$ is denoted $\theta_i$ in medium $i$. In the following, for simplicity, we will continue omitting this spatial variation, which is justified since the discontinuity will be written specifically at $z = 0$, and the $xy$ variation remains the same in all media. Using the notion of an effective index ($\tilde{n}$) to relate the tangential components of the Fourier transforms of the electric and magnetic fields in a homogeneous medium, i.e., in TE polarization, we establish:

$$ 
\begin{aligned}
\tilde{H}_i^{\pm} &= \pm \tilde{n}_i z \times \tilde{E}_i^{\pm}, \\
\tilde{n}_i &= \frac{\alpha_i}{\omega \mu_i} = \left( \frac{1}{\eta_0 \mu_r^i} \right) n_i \cos \theta_i, \\
\eta_0 &= \sqrt{\frac{\mu_0}{\varepsilon_0}},
\end{aligned}
\quad (16) $$

where $\mu_r^i$ is the relative permeability and $\eta_0$ is the vacuum impedance. The $+$ sign must be taken for progressive waves, and the $-$ sign for retrograde waves. Finally, Eqs. (14) and (15), combined with Eq. (16), give for TE polarization:

$$ \frac{\tilde{M}}{A_1^+(\omega_1)} = \delta(\omega - \omega_1) - t(\omega, \omega_1) + r(\omega, \omega_1), \quad (17) $$

$$ \frac{\tilde{J}}{A_1^+(\omega_1)} = \tilde{n}_1(\omega)[\delta(\omega - \omega_1) - r(\omega, \omega_1)] - \tilde{n}_2(\omega)t(\omega, \omega_1). \quad (18) $$

The next step is to develop the currents $M$ and $J$. By the virtue of the convolution (2), we write: $\tilde{J}(\omega) = -\varepsilon_0 \mathcal{I}_e$, and $\tilde{M}(\omega) = \mu_0 \mathcal{I}_h$, with

$$ \mathcal{I}_u = 2\pi i\omega \int \tilde{\chi}_u(\omega', \omega - \omega') \tilde{X}_u(\omega - \omega') d\omega', \quad (19) $$

where $u = e, h$, and $\tilde{\chi}_u$ denotes the double Fourier transform of the susceptibility. A question that now arises is what fields to choose in these convolutions that are integrals of discontinuous quantities taken at $z = 0$. To overcome this difficulty, we invoke a suitable mean field $X_{ua} = [X_{u1}(t) + X_{u2}(t)]/2$, and thus obtain:

$$ \frac{2\tilde{E}_a}{A_1^+(\omega_1)} = r(\omega, \omega_1) + t(\omega, \omega_1) + \delta(\omega - \omega_1), \quad (20) $$
$$ \frac{2\tilde{H}_a}{A_1^+(\omega_1)} = \tilde{n}_1 r(\omega, \omega_1) - \tilde{n}_2 t(\omega, \omega_1) - \tilde{n}_1 \delta(\omega - \omega_1). \quad (21) $$

These average fields will be introduced in Eq. (19), which will allow the development of the currents in Eqs. (17) and (18). It will be assumed from here that the index dispersion of the substrate and superstrate is negligible [$\tilde{n}_i(\omega) = \tilde{n}_i$], especially since our study will be reduced to a moderate spectral window.

## IV. FREQUENCY CONVERSION

To achieve a linear-frequency conversion, we must meet the conditions to reduce the continuum. Our frequency conversion objective to limit the continuum to discrete frequencies requires the reflection and transmission coefficients to be written in the form of Dirac distributions. We will therefore seek here a condition such that the reflected field is reduced to a harmonic of frequency $\omega_1$ (identical to the incident field), and the transmitted field to an arbitrary harmonic of frequency $\omega_2$ (see Fig. 1), or:

$$ r(\omega, \omega_1) = r_0(\omega_1, \omega_1)\delta(\omega - \omega_1), \quad (22) $$

$$ t(\omega, \omega_1) = t_0(\omega_1, \omega_2)\delta(\omega - \omega_2). \quad (23) $$

This allows the development of the average field in the form of Dirac distributions, which simplifies the inhomogeneous convolution product given for the polarization currents in (19). However, this simplification also requires the susceptibility to be written in the form of Dirac distributions (which will be discussed in Sec. V), that is:

$$ \hat{\chi}_u(\omega - \omega_1, \omega_1) = \hat{\chi}_u^{11}\delta(\omega - \omega_1) + \hat{\chi}_u^{12}\delta(\omega - \omega_2), \quad (24) $$

$$ \hat{\chi}_u(\omega - \omega_2, \omega_2) = \hat{\chi}_u^{21}\delta(\omega - \omega_1) + \hat{\chi}_u^{22}\delta(\omega - \omega_2), \quad (25) $$

using which, Eq. (19), Eq. (17), and Eq. (18), we get:

$$ \frac{\tilde{J}(\omega)}{i\pi\epsilon_0 A_1^+} = -(1 + r_0)(\Omega_e^{11} + \Omega_e^{12}) - t_0(\Omega_e^{21} + \Omega_e^{22}), \quad (26) $$

$$ \frac{\tilde{M}(\omega)}{i\pi \mu_0 A_1^+} = \tilde{n}_1(r_0 - 1)(\Omega_h^{11} + \Omega_h^{12}) - \tilde{n}_2 t_0 (\Omega_h^{21} + \Omega_h^{22}), \quad (27) $$

where

$$ \Omega_u^{11} = \omega_1 \hat{\chi}_u^{11} \delta(\omega - \omega_1), \quad \Omega_u^{12} = \omega_2 \hat{\chi}_u^{12} \delta(\omega - \omega_2), \quad (28) $$

$$ \Omega_u^{21} = \omega_1 \hat{\chi}_u^{21} \delta(\omega - \omega_1), \quad \Omega_u^{22} = \omega_2 \hat{\chi}_u^{22} \delta(\omega - \omega_2). \quad (29) $$

It then remains to identify the coefficients, which weigh the Dirac distributions after identification of Eq. (26), Eq. (27), and Eq. (17), Eq. (18), that is:

$$ \tilde{n}_1(r_0 - 1)[\omega_1 \hat{\chi}_h^{11} \delta(\omega - \omega_1) + \omega_2 \hat{\chi}_h^{12} \delta(\omega - \omega_2)] - \tilde{n}_2 t_0 [\omega_1 \hat{\chi}_h^{21} \delta(\omega - \omega_1) + \omega_2 \hat{\chi}_h^{22} \delta(\omega - \omega_2)] \\ = \frac{1}{i\pi \mu_0} \{(1 + r_0)\delta(\omega - \omega_1) - t_0 \delta(\omega - \omega_2)\}, \quad (30) $$

$$ (1 + r_0)[\omega_1 \hat{\chi}_e^{11} \delta(\omega - \omega_1) + \omega_2 \hat{\chi}_e^{12} \delta(\omega - \omega_2)] + t_0[\omega_1 \hat{\chi}_e^{21} \delta(\omega - \omega_1) + \omega_2 \hat{\chi}_e^{22} \delta(\omega - \omega_2)] \\ = \frac{1}{i\pi \epsilon_0} \{\tilde{n}_1(r_0 - 1)\delta(\omega - \omega_1) + \tilde{n}_2 t_0 \delta(\omega - \omega_2)\}. \quad (31) $$

Each time, for the magnetic aspect of Eq. (30), and for the electrical aspect of Eq. (31), we have two equations involving four susceptibility coefficients. These equations will make it possible to express the reflection and transmission coefficients as a function of the susceptibility parameters, and also to highlight the constraints on these parameters. Following Eqs. (30) and (31), we get a system of four equations with two unknowns ($r_0, t_0$) and eight susceptibility coefficients. With the magnetic parameters of Eq. (30), we first obtain:

$$ r_0 \left\{ \tilde{n}_1 \hat{\chi}_h^{11} - \frac{1}{\gamma \omega_1} \right\} + t_0 \left\{ -\tilde{n}_2 \hat{\chi}_h^{21} \right\} = \left\{ \tilde{n}_1 \hat{\chi}_h^{11} + \frac{1}{\gamma \omega_1} \right\}, \quad (32) $$

$$ r_0 \{\tilde{n}_1 \hat{\chi}_h^{12}\} + t_0 \left\{ -\tilde{n}_2 \hat{\chi}_h^{22} + \frac{1}{\gamma \omega_2} \right\} = \{\tilde{n}_1 \hat{\chi}_h^{12}\}, \quad (33) $$

and, following Eq. (31) with the electric parameters:

$$ r_0 \left\{ \hat{\chi}_e^{11} - \frac{\tilde{n}_1}{\gamma' \omega_1} \right\} + t_0 \left\{ \hat{\chi}_e^{21} \right\} = - \left\{ \hat{\chi}_e^{11} + \frac{\tilde{n}_1}{\gamma' \omega_1} \right\}, \quad (34) $$

$$ r_0 \{\hat{\chi}_e^{12}\} + t_0 \left\{ \hat{\chi}_e^{22} - \frac{\tilde{n}_2}{\gamma' \omega_2} \right\} = - \hat{\chi}_e^{12}, \quad (35) $$

where

$$ (\gamma, \gamma') = i\pi(\mu_0, \epsilon_0). \quad (36) $$

At this point, all that remains is to identify the image space of the reflection coefficient ($r_0 = r_e = r_h$) and transmission coefficient ($t_0 = t_e = t_h$). Note here that ($r_e, t_e$) are calculated from the electric parameters in Eq. (34), Eq. (35), while ($r_h, t_h$) are calculated from the magnetic parameters in Eqs. (32) and (33). Because these solutions must be identical, we get a relationship between the electric and magnetic parameters, that is:

$$ \hat{\chi}_h^{22} = \frac{1}{i\pi \mu_0 \omega_2 \tilde{n}} + \left( \frac{r_e - 1}{t_e} \right) \hat{\chi}_h^{12}, \quad (37) $$

$$ \hat{\chi}_h^{11} = \frac{1}{i\pi \mu_0 \omega_1 \tilde{n}} \left( \frac{1+r_e}{r_e-1} \right) + \left( \frac{t_e}{r_e-1} \right) \hat{\chi}_h^{21}, \quad (38) $$
$$ \frac{2\tilde{H}_a}{A_1^+(\omega_1)} = \tilde{n}_1 r(\omega, \omega_1) - \tilde{n}_2 t(\omega, \omega_1) - \tilde{n}_1 \delta(\omega - \omega_1). \quad (21) $$

These average fields will be introduced in Eq. (19), which will allow the development of the currents in Eqs. (17) and (18). It will be assumed from here that the index dispersion of the substrate and superstrate is negligible [$\tilde{n}_i(\omega) = \tilde{n}_i$], especially since our study will be reduced to a moderate spectral window.

## IV. FREQUENCY CONVERSION

To achieve a linear-frequency conversion, we must meet the conditions to reduce the continuum. Our frequency conversion objective to limit the continuum to discrete frequencies requires the reflection and transmission coefficients to be written in the form of Dirac distributions. We will therefore seek here a condition such that the reflected field is reduced to a harmonic of frequency $\omega_1$ (identical to the incident field), and the transmitted field to an arbitrary harmonic of frequency $\omega_2$ (see Fig. 1), or:

$$ r(\omega, \omega_1) = r_0(\omega_1, \omega_1)\delta(\omega - \omega_1), \quad (22) $$

$$ t(\omega, \omega_1) = t_0(\omega_1, \omega_2)\delta(\omega - \omega_2). \quad (23) $$

This allows the development of the average field in the form of Dirac distributions, which simplifies the inhomogeneous convolution product given for the polarization currents in (19). However, this simplification also requires the susceptibility to be written in the form of Dirac distributions (which will be discussed in Sec. V), that is:

$$ \hat{\chi}_u(\omega - \omega_1, \omega_1) = \hat{\chi}_u^{11}\delta(\omega - \omega_1) + \hat{\chi}_u^{12}\delta(\omega - \omega_2), \quad (24) $$

$$ \hat{\chi}_u(\omega - \omega_2, \omega_2) = \hat{\chi}_u^{21}\delta(\omega - \omega_1) + \hat{\chi}_u^{22}\delta(\omega - \omega_2), \quad (25) $$

using which, Eq. (19), Eq. (17), and Eq. (18), we get:

$$ \frac{\tilde{J}(\omega)}{i\pi\epsilon_0 A_1^+} = -(1 + r_0)(\Omega_e^{11} + \Omega_e^{12}) - t_0(\Omega_e^{21} + \Omega_e^{22}), \quad (26) $$

$$ \frac{\tilde{M}(\omega)}{i\pi \mu_0 A_1^+} = \tilde{n}_1(r_0 - 1)(\Omega_h^{11} + \Omega_h^{12}) - \tilde{n}_2 t_0 (\Omega_h^{21} + \Omega_h^{22}), \quad (27) $$

where

$$ \Omega_u^{11} = \omega_1 \hat{\chi}_u^{11} \delta(\omega - \omega_1), \quad \Omega_u^{12} = \omega_2 \hat{\chi}_u^{12} \delta(\omega - \omega_2), \quad (28) $$

$$ \Omega_u^{21} = \omega_1 \hat{\chi}_u^{21} \delta(\omega - \omega_1), \quad \Omega_u^{22} = \omega_2 \hat{\chi}_u^{22} \delta(\omega - \omega_2). \quad (29) $$

It then remains to identify the coefficients, which weigh the Dirac distributions after identification of Eq. (26), Eq. (27), and Eq. (17), Eq. (18), that is:

$$ \tilde{n}_1(r_0 - 1)[\omega_1 \hat{\chi}_h^{11} \delta(\omega - \omega_1) + \omega_2 \hat{\chi}_h^{12} \delta(\omega - \omega_2)] - \tilde{n}_2 t_0 [\omega_1 \hat{\chi}_h^{21} \delta(\omega - \omega_1) + \omega_2 \hat{\chi}_h^{22} \delta(\omega - \omega_2)] \\ = \frac{1}{i\pi \mu_0} \{(1 + r_0)\delta(\omega - \omega_1) - t_0 \delta(\omega - \omega_2)\}, \quad (30) $$

$$ (1 + r_0)[\omega_1 \hat{\chi}_e^{11} \delta(\omega - \omega_1) + \omega_2 \hat{\chi}_e^{12} \delta(\omega - \omega_2)] + t_0[\omega_1 \hat{\chi}_e^{21} \delta(\omega - \omega_1) + \omega_2 \hat{\chi}_e^{22} \delta(\omega - \omega_2)] \\ = \frac{1}{i\pi \epsilon_0} \{\tilde{n}_1(r_0 - 1)\delta(\omega - \omega_1) + \tilde{n}_2 t_0 \delta(\omega - \omega_2)\}. \quad (31) $$

Each time, for the magnetic aspect of Eq. (30), and for the electrical aspect of Eq. (31), we have two equations involving four susceptibility coefficients. These equations will make it possible to express the reflection and transmission coefficients as a function of the susceptibility parameters, and also to highlight the constraints on these parameters. Following Eqs. (30) and (31), we get a system of four equations with two unknowns ($r_0, t_0$) and eight susceptibility coefficients. With the magnetic parameters of Eq. (30), we first obtain:

$$ r_0 \left\{ \tilde{n}_1 \hat{\chi}_h^{11} - \frac{1}{\gamma \omega_1} \right\} + t_0 \left\{ -\tilde{n}_2 \hat{\chi}_h^{21} \right\} = \left\{ \tilde{n}_1 \hat{\chi}_h^{11} + \frac{1}{\gamma \omega_1} \right\}, \quad (32) $$

$$ r_0 \{ \tilde{n}_1 \hat{\chi}_h^{12} \} + t_0 \{ -\tilde{n}_2 \hat{\chi}_h^{22} + \frac{1}{\gamma \omega_2} \} = \{ \tilde{n}_1 \hat{\chi}_h^{12} \}, \quad (33) $$

and, following Eq. (31) with the electric parameters:

$$ r_0 \left\{ \hat{\chi}_e^{11} - \frac{\tilde{n}_1}{\gamma' \omega_1} \right\} + t_0 \left\{ \hat{\chi}_e^{21} \right\} = - \left\{ \hat{\chi}_e^{11} + \frac{\tilde{n}_1}{\gamma' \omega_1} \right\}, \quad (34) $$

$$ r_0 \{ \hat{\chi}_e^{12} \} + t_0 \{ \hat{\chi}_e^{22} - \frac{\tilde{n}_2}{\gamma' \omega_2} \} = - \hat{\chi}_e^{12}, \quad (35) $$

where

$$ (\gamma, \gamma') = i\pi(\mu_0, \epsilon_0). \quad (36) $$

At this point, all that remains is to identify the image space of the reflection coefficient ($r_0 = r_e = r_h$) and transmission coefficient ($t_0 = t_e = t_h$). Note here that ($r_e, t_e$) are calculated from the electric parameters in Eq. (34), Eq. (35), while ($r_h, t_h$) are calculated from the magnetic parameters in Eqs. (32) and (33). Because these solutions must be identical, we get a relationship between the electric and magnetic parameters, that is:

$$ \hat{\chi}_h^{22} = \frac{1}{i\pi \mu_0 \omega_2 \tilde{n}} + \left( \frac{r_e - 1}{t_e} \right) \hat{\chi}_h^{12}, \quad (37) $$

$$ \hat{\chi}_h^{11} = \frac{1}{i\pi \mu_0 \omega_1 \tilde{n}} \left( \frac{1+r_e}{r_e-1} + \left( \frac{t_e}{r_e-1} \right) \hat{\chi}_h^{21} \right), \quad (38) $$
where $(r_e, t_e)$ are calculated from the four electric coefficients, and:

$$ \hat{\chi}_e^{22} = \frac{\tilde{n}}{i\pi\epsilon_0\omega_2} - \left(\frac{1+r_h}{t_h}\right) \hat{\chi}_e^{12}, \quad (39) $$

$$ \hat{\chi}_e^{11} = \frac{-\tilde{n}}{i\pi\epsilon_0\omega_1} \left( \frac{1-r_h}{1+r_h} \right) - \left( \frac{t_h}{1+r_h} \right) \hat{\chi}_e^{21}, \quad (40) $$

where $(r_h, t_h)$ are calculated from the four magnetic coefficients. This shows that the number of free parameters is reduced to six. More precisely, $(r_e, t_e)$ can be calculated from Eq. (34), Eq. (35) with four electric coefficients, and the resulting value of $(r_e, t_e)$ reduces the number of free magnetic parameters to two [see Eqs. (37) and (38)]. In a similar way, $(r_h, t_h)$ can be calculated from Eqs. (32), and (33) with four magnetic coefficients, and the resulting value of $(r_h, t_h)$ reduces the number of free electric parameters to two [see Eqs. (39) and (40)]. Having said that, we also need to know in which range of values it is interesting to vary the parameters to observe notable variations in the pair $(r_0, t_0)$. To do this, we solve Eqs. (32)–(35) while highlighting the dimension of the susceptibility, as:

$$ r_h = \frac{\left(\hat{\chi}_h^{11} + \frac{1}{\tilde{n}\gamma\omega_1}\right)\left(-\hat{\chi}_h^{22} + \frac{1}{\tilde{n}\gamma\omega_2}\right) + \hat{\chi}_h^{12}\hat{\chi}_h^{21}}{\left(\hat{\chi}_h^{11} - \frac{1}{\tilde{n}\gamma\omega_1}\right)\left(-\hat{\chi}_h^{22} + \frac{1}{\tilde{n}\gamma\omega_2}\right) + \hat{\chi}_h^{12}\hat{\chi}_h^{21}}, \quad (41) $$

$$ t_h = \frac{\frac{1}{\tilde{n}\gamma\omega_1}\hat{\chi}_h^{12}}{-\hat{\chi}_h^{12}\hat{\chi}_h^{21} + \left(\hat{\chi}_h^{22} - \frac{1}{\tilde{n}\gamma\omega_2}\right)\left(\hat{\chi}_h^{11} - \frac{1}{\tilde{n}\gamma\omega_1}\right)}, \quad (42) $$

$$ r_e = \frac{-\left(\hat{\chi}_e^{11} + \frac{\tilde{n}}{\gamma'\omega_1}\right)\left(-\hat{\chi}_e^{22} + \frac{\tilde{n}}{\gamma'\omega_2}\right) - \hat{\chi}_e^{12}\hat{\chi}_e^{21}}{\left(\hat{\chi}_e^{11} - \frac{\tilde{n}}{\gamma'\omega_1}\right)\left(-\hat{\chi}_e^{22} + \frac{\tilde{n}}{\gamma'\omega_2}\right) + \hat{\chi}_e^{12}\hat{\chi}_e^{21}}, \quad (43) $$

$$ t_e = \frac{\frac{-2\tilde{n}}{\gamma'\omega_1}\hat{\chi}_e^{12}}{-\hat{\chi}_e^{12}\hat{\chi}_e^{21} + \left(\hat{\chi}_e^{22} - \frac{\tilde{n}}{\gamma'\omega_2}\right)\left(\hat{\chi}_e^{11} - \frac{\tilde{n}}{\gamma'\omega_1}\right)}. \quad (44) $$

With $\theta$ defined in Fig. 1, we note for the magnetic coefficients in Eqs. (41) and (42), that:

$$ \frac{1}{\tilde{n}\gamma\omega_i} = -\frac{i}{2\pi^2} \frac{\mu_r}{n \cos \theta} \lambda_i, \quad (45) $$

and for the electric coefficients in Eqs. (43) and (44):

$$ \frac{\tilde{n}}{\gamma'\omega_i} = -\frac{i}{2\pi^2} \frac{n \cos \theta}{\mu_r} \lambda_i. \quad (46) $$

Furthermore, for the links between the susceptibility parameters, we obtain:

$$ \hat{\chi}_h^{22} = \frac{-i}{2\pi^2} \frac{\mu_r}{n \cos \theta} \lambda_2 + \left(\frac{r_e - 1}{t_e}\right) \hat{\chi}_h^{12}, \quad (47) $$

$$ \hat{\chi}_h^{11} = \frac{-i}{2\pi^2} \frac{\mu_r}{n \cos \theta} \lambda_1 \left( \frac{1+r_e}{r_e-1} \right) + \left( \frac{t_e}{r_e-1} \right) \hat{\chi}_h^{21}, \quad (48) $$

$$ \hat{\chi}_e^{22} = \frac{-i}{2\pi^2} \frac{n \cos \theta}{\mu_r} \lambda_2 - \left( \frac{1+r_h}{t_h} \right) \hat{\chi}_e^{12}, \quad (49) $$

$$ \hat{\chi}_e^{11} = \frac{i}{2\pi^2} \frac{n \cos \theta}{\mu_r} \lambda_1 \left( \frac{1-r_h}{1+r_h} \right) - \left( \frac{t_h}{1+r_h} \right) \hat{\chi}_e^{21}. \quad (50) $$

For all these equations, the order of magnitude of the coefficients is the wavelength $\lambda_1$ or $\lambda_2$, at least for $r \neq \pm 1$ and $t \neq 0$. Consequently, for any numerical calculation, it will be

interesting to zoom in to the vicinity of $|\hat{\chi}_h^{ij}| \approx \lambda_{1,2}$, or more precisely between $0.1\lambda_{1,2}$ and $10\lambda_{1,2}$.

## V. DISCUSSIONS

Analysis of the previous equations shows that the space of susceptibility parameters, which lead to the same value of $(r_0, t_0)$ is of dimension 4. To show that, one can write the equations at a given couple $(r_0, t_0)$, with the susceptibility parameters as the unknowns. Therefore, apart from the couple $(r_0 = \pm 1, t_0 = 0)$, any pair of value $(r_0, t_0)$ can be achieved using the susceptibility coefficients. However, the arbitrary values that can take $(r_0, t_0)$ do not disagree with the energy balance that we discuss now. Contrary to usual, the values $R_h = |r_h|^2$ and $T_h = |t_h|^2$ are not bounded by unity. This is because the reflected ($\phi_r$) and transmitted ($\phi_t$) fluxes have been conventionally defined as if we were illuminating a passive component from the outside, that is, as $\phi_r = R\phi_i$ and $\phi_t = T\phi_i$, with $\phi_i$ being the incident flux. This does not work anymore with an active metasurface, since the reflected and transmitted fluxes are provided by the sum of the incident flux $\phi_i$ and the external power $F_{\text{ext}}$ supplied to the metasurface. In other words, the energy balance is written as:

$$ F_{\text{ext}} + \phi_i = \phi_r + \phi_t = \phi_i (R_h + T_h). \quad (51) $$

In the classical case (no external action), we have $F_{\text{ext}} = 0$ and we find $R_h + T_h = 1$, bearing in mind that absorption is zero because of the negligible thickness of the metasurface. Here, in stating the negligible absorption, we do not ignore the possibility of other surface modes that might be supported by a metasurface. While this should be interesting to explore, it does not present any obstacle to achieving the presented conversion objectives. By devising a metasurface that does not exhibit interfering modes (unless beneficial for the conversion), we may either avoid excitation of such modes or design a metasurface that forbids those modes. The interaction dynamics of light with matter can be intensified under certain resonant conditions, even with negligible material thickness. While the metasurface's thinness does reduce volumetric absorption, the exact amount is also contingent on specific resonant modes, if present at all. We emphasize that the design of our metasurface can be tailored to minimize such resonances that could potentially hinder the desired frequency conversion. Given the tunable nature of metasurface design, this requirement does not seem to present a prohibitive factor. The situation is not unlike noble metal nanosystems that support resonant excitation of plasmons [21], or unlike quasicrystals and all-dielectric stacks [22] that support Bloch waves or sharp resonances. In all such systems, while material volume is small, strong but tunable resonant excitations can be tailored. Following Eq. (51), with the active metasurface the two fluxes ($\phi_r, \phi_t$) must be normalized by the total energy ($F_{\text{ext}} + \phi_i$), which supplied, and this leads to the introduction of the normalized coefficients:

$$ (R_N, T_N) = \frac{1}{R_h + T_h} (R_h, T_h). \quad (52) $$

Preliminary numerical calculations of the fluxes for the simplified case of $\hat{\chi}_e^{ij} = \hat{\chi}_e$ are visualized in Figs. 2 and 3. This calculation is carried out for an illumination wavelength
where $(r_e, t_e)$ are calculated from the four electric coefficients, and:

$$ \hat{\chi}_e^{22} = \frac{\tilde{n}}{i\pi \varepsilon_0 \omega_2} - \left(\frac{1+r_h}{t_h}\right) \hat{\chi}_e^{12}, \quad (39) $$

$$ \hat{\chi}_e^{11} = \frac{-\tilde{n}}{i\pi \varepsilon_0 \omega_1} \left( \frac{1-r_h}{1+r_h} \right) - \left( \frac{t_h}{1+r_h} \right) \hat{\chi}_e^{21}, \quad (40) $$

where $(r_h, t_h)$ are calculated from the four magnetic coefficients. This shows that the number of free parameters is reduced to six. More precisely, $(r_e, t_e)$ can be calculated from Eq. (34), Eq. (35) with four electric coefficients, and the resulting value of $(r_e, t_e)$ reduces the number of free magnetic parameters to two [see Eqs. (37) and (38)]. In a similar way, $(r_h, t_h)$ can be calculated from Eqs. (32), and (33) with four magnetic coefficients, and the resulting value of $(r_h, t_h)$ reduces the number of free electric parameters to two [see Eqs. (39) and (40)]. Having said that, we also need to know in which range of values it is interesting to vary the parameters to observe notable variations in the pair $(r_0, t_0)$. To do this, we solve Eqs. (32)–(35) while highlighting the dimension of the susceptibility, as:

$$ r_h = \frac{\left(\hat{\chi}_h^{11} + \frac{1}{\tilde{n}\gamma\omega_1}\right)\left(-\hat{\chi}_h^{22} + \frac{1}{\tilde{n}\gamma\omega_2}\right) + \hat{\chi}_h^{12}\hat{\chi}_h^{21}}{\left(\hat{\chi}_h^{11} - \frac{1}{\tilde{n}\gamma\omega_1}\right)\left(-\hat{\chi}_h^{22} + \frac{1}{\tilde{n}\gamma\omega_2}\right) + \hat{\chi}_h^{12}\hat{\chi}_h^{21}}, \quad (41) $$

$$ t_h = \frac{\frac{1}{\tilde{n}\gamma\omega_1}\hat{\chi}_h^{12}}{-\hat{\chi}_h^{12}\hat{\chi}_h^{21} + \left(\hat{\chi}_h^{22} - \frac{1}{\tilde{n}\gamma\omega_2}\right)\left(\hat{\chi}_h^{11} - \frac{1}{\tilde{n}\gamma\omega_1}\right)}, \quad (42) $$

$$ r_e = \frac{-\left(\hat{\chi}_e^{11} + \frac{\tilde{n}}{\gamma'\omega_1}\right)\left(-\hat{\chi}_e^{22} + \frac{\tilde{n}}{\gamma'\omega_2}\right) - \hat{\chi}_e^{12}\hat{\chi}_e^{21}}{\left(\hat{\chi}_e^{11} - \frac{\tilde{n}}{\gamma'\omega_1}\right)\left(-\hat{\chi}_e^{22} + \frac{\tilde{n}}{\gamma'\omega_2}\right) + \hat{\chi}_e^{12}\hat{\chi}_e^{21}}, \quad (43) $$

$$ t_e = \frac{\frac{-2\tilde{n}}{\gamma'\omega_1}\hat{\chi}_e^{12}}{-\hat{\chi}_e^{12}\hat{\chi}_e^{21} + \left(\hat{\chi}_e^{22} - \frac{\tilde{n}}{\gamma'\omega_2}\right)\left(\hat{\chi}_e^{11} - \frac{\tilde{n}}{\gamma'\omega_1}\right)}. \quad (44) $$

With $\theta$ defined in Fig. 1, we note for the magnetic coefficients in Eqs. (41) and (42), that:

$$ \frac{1}{\tilde{n}\gamma\omega_i} = -\frac{i}{2\pi^2} \frac{\mu_r}{n \cos \theta} \lambda_i, \quad (45) $$

and for the electric coefficients in Eqs. (43) and (44):

$$ \frac{\tilde{n}}{\gamma'\omega_i} = -\frac{i}{2\pi^2} \frac{n \cos \theta}{\mu_r} \lambda_i. \quad (46) $$

Furthermore, for the links between the susceptibility parameters, we obtain:

$$ \hat{\chi}_h^{22} = \frac{-i}{2\pi^2} \frac{\mu_r}{n \cos \theta} \lambda_2 + \left(\frac{r_e - 1}{t_e}\right) \hat{\chi}_h^{12}, \quad (47) $$

$$ \hat{\chi}_h^{11} = \frac{-i}{2\pi^2} \frac{\mu_r}{n \cos \theta} \lambda_1 \left( \frac{1+r_e}{r_e-1} \right) + \left( \frac{t_e}{r_e-1} \right) \hat{\chi}_h^{21}, \quad (48) $$

$$ \hat{\chi}_e^{22} = \frac{-i}{2\pi^2} \frac{n \cos \theta}{\mu_r} \lambda_2 - \left( \frac{1+r_h}{t_h} \right) \hat{\chi}_e^{12}, \quad (49) $$

$$ \hat{\chi}_e^{11} = \frac{i}{2\pi^2} \frac{n \cos \theta}{\mu_r} \lambda_1 \left( \frac{1-r_h}{1+r_h} \right) - \left( \frac{t_h}{1+r_h} \right) \hat{\chi}_e^{21}. \quad (50) $$

For all these equations, the order of magnitude of the coefficients is the wavelength $\lambda_1$ or $\lambda_2$, at least for $r \neq \pm 1$ and $t \neq 0$. Consequently, for any numerical calculation, it will be

interesting to zoom in to the vicinity of $|\hat{\chi}_h^{ij}| \approx \lambda_{1,2}$, or more precisely between $0.1\lambda_{1,2}$ and $10\lambda_{1,2}$.

## V. DISCUSSIONS

Analysis of the previous equations shows that the space of susceptibility parameters, which lead to the same value of $(r_0, t_0)$ is of dimension 4. To show that, one can write the equations at a given couple $(r_0, t_0)$, with the susceptibility parameters as the unknowns. Therefore, apart from the couple $(r_0 = \pm 1, t_0 = 0)$, any pair of value $(r_0, t_0)$ can be achieved using the susceptibility coefficients. However, the arbitrary values that can take $(r_0, t_0)$ do not disagree with the energy balance that we discuss now. Contrary to usual, the values $R_h = |r_h|^2$ and $T_h = |t_h|^2$ are not bounded by unity. This is because the reflected ($\phi_r$) and transmitted ($\phi_t$) fluxes have been conventionally defined as if we were illuminating a passive component from the outside, that is, as $\phi_r = R\phi_i$ and $\phi_t = T\phi_i$, with $\phi_i$ being the incident flux. This does not work anymore with an active metasurface, since the reflected and transmitted fluxes are provided by the sum of the incident flux $\phi_i$ and the external power $F_{\text{ext}}$ supplied to the metasurface. In other words, the energy balance is written as:

$$ F_{\text{ext}} + \phi_i = \phi_r + \phi_t = \phi_i (R_h + T_h). \quad (51) $$

In the classical case (no external action), we have $F_{\text{ext}} = 0$ and we find $R_h + T_h = 1$, bearing in mind that absorption is zero because of the negligible thickness of the metasurface. Here, in stating the negligible absorption, we do not ignore the possibility of other surface modes that might be supported by a metasurface. While this should be interesting to explore, it does not present any obstacle to achieving the presented conversion objectives. By devising a metasurface that does not exhibit interfering modes (unless beneficial for the conversion), we may either avoid excitation of such modes or design a metasurface that forbids those modes. The interaction dynamics of light with matter can be intensified under certain resonant conditions, even with negligible material thickness. While the metasurface’s thinness does reduce volumetric absorption, the exact amount is also contingent on specific resonant modes, if present at all. We emphasize that the design of our metasurface can be tailored to minimize such resonances that could potentially hinder the desired frequency conversion. Given the tunable nature of metasurface design, this requirement does not seem to present a prohibitive factor. The situation is not unlike noble metal nanosystems that support resonant excitation of plasmons [21], or unlike quasicrystals and all-dielectric stacks [22] that support Bloch waves or sharp resonances. In all such systems, while material volume is small, strong but tunable resonant excitations can be tailored. Following Eq. (51), with the active metasurface the two fluxes ($\phi_r, \phi_t$) must be normalized by the total energy ($F_{\text{ext}} + \phi_i$), which supplied, and this leads to the introduction of the normalized coefficients:

$$ (R_N, T_N) = \frac{1}{R_h + T_h} (R_h, T_h). \quad (52) $$

Preliminary numerical calculations of the fluxes for the simplified case of $\hat{\chi}_e^{ij} = \hat{\chi}_e$ are visualized in Figs. 2 and 3. This calculation is carried out for an illumination wavelength
<figure><img src="image_3.png" /><figcaption>FIG. 2. Normalized fluxes relative to the wavelength ratio $\lambda_2/\lambda_1$ simulated for an incident illumination wavelength $\lambda_1 = 1000$ nm. All susceptibility parameters are identical. At the condition where $\lambda_2 = \lambda_1$, indicating no frequency conversion, the transmission is at its theoretical maximum, signifying the complete passage of light with no reflection due to the matching media on either side of the metasurface.</figcaption></figure>

$\lambda_1 = 1000$ nm, while the conversion wavelength is $\lambda_2$. We observe in Fig. 2 at $\lambda_2 = \lambda_1$ (no frequency conversion), the normalized transmission is equal to 1, given that the extreme media are identical. From this value, conversion efficiency decreases for $\lambda_2 < \lambda_1$, then stabilizes at 80%. On the other side ($\lambda_2 > \lambda_1$), the efficiency decreases to zero. However, we must keep in mind that this is a specific result of equal

susceptibility parameters (all equal to $\lambda_1$ in Fig. 2), while we should explore the four electric coefficients. In Fig. 3, the conversion efficiency is plotted versus the two wavelengths. We first observe that the straight line $\lambda_1 = \lambda_2$ is in accordance with the previous remark ($T = 1, R = 0$). We also observe that the range of conversion wavelengths $\lambda_2$ by transmission is reduced at small $\lambda_1$ wavelengths.

Here, it is worth noting that the frequency $\omega_2$ is understood to arise in the system via the modulation of the metasurface properties in time. This could be analogous to an external pump in a pump-probe setup, where the metasurface undergoes dynamic changes influenced by a separate control beam or electromagnetic influence, effectively introducing $\omega_2$. Relevant temporal control schemes have been reported, e.g., based on the concept of thermoplasmonics [21] (in which the dielectric function of a thin film undergoes a modulation), further underscoring the feasibility of this approach.

So far we have not discussed the temporal form of the susceptibility, the relevance, and existence of which must be verified. Starting from (24)–(25), we can immediately write, after an inverse Fourier transformation:

$$ \tilde{\chi}_u(t, \omega_1) = \hat{\chi}_u^{11}(\omega_1, \omega_1) + e^{-i(\omega_2 - \omega_1)t} \hat{\chi}_u^{12}(\omega_1, \omega_2), \quad (53) $$

$$ \tilde{\chi}_u(t, \omega_2) = e^{i(\omega_2 - \omega_1)t} \hat{\chi}_u^{21}(\omega_2, \omega_1) + \hat{\chi}_u^{22}(\omega_2, \omega_2), \quad (54) $$

implying that the time shape of the Fourier transform of the susceptibility with respect to $t'$, i.e., $\tilde{\chi}_u(t, \omega)$, must be imposed at the two frequencies $\omega_i$ involved in the conversion. We note that this condition is only imposed at these two frequencies. In this context, we therefore need to check that this condition is realistic, i.e., compatible with causality. In fact $\chi(t, t')$ is causal versus $t'$, so that $\tilde{\chi}_u(t, \omega)$ is the Fourier transform of a causal function. In this context, we impose two values $\tilde{\chi}_u(t, \omega_i)$ on the Fourier transform of a causal function. We

<figure><img src="image_4.png" /><figcaption>FIG. 3. Normalized fluxes for reflection and transmission contoured against the incident wavelength $\lambda_1$ and conversion wavelength $\lambda_2$. The efficiency of the metasurface is represented by the color gradient, with yellow for high and blue for low efficiency. The diagonal where $\lambda_1 = \lambda_2$ marks the condition of no frequency conversion with maximum transmission ($T = 1$) and no reflection ($R = 0$), denoting the state of the metasurface where the incident and converted wavelengths are identical. The maps highlight how the conversion efficiency diminishes when $\lambda_2 > \lambda_1$ and reveal a restricted conversion range for shorter incident wavelengths.</figcaption></figure>
<figure><img src="image_5.png" /><figcaption>FIG. 2. Normalized fluxes relative to the wavelength ratio $\lambda_2/\lambda_1$ simulated for an incident illumination wavelength $\lambda_1 = 1000$ nm. All susceptibility parameters are identical. At the condition where $\lambda_2 = \lambda_1$, indicating no frequency conversion, the transmission is at its theoretical maximum, signifying the complete passage of light with no reflection due to the matching media on either side of the metasurface.</figcaption></figure>

$\lambda_1 = 1000$ nm, while the conversion wavelength is $\lambda_2$. We observe in Fig. 2 at $\lambda_2 = \lambda_1$ (no frequency conversion), the normalized transmission is equal to 1, given that the extreme media are identical. From this value, conversion efficiency decreases for $\lambda_2 < \lambda_1$, then stabilizes at 80%. On the other side ($\lambda_2 > \lambda_1$), the efficiency decreases to zero. However, we must keep in mind that this is a specific result of equal

susceptibility parameters (all equal to $\lambda_1$ in Fig. 2), while we should explore the four electric coefficients. In Fig. 3, the conversion efficiency is plotted versus the two wavelengths. We first observe that the straight line $\lambda_1 = \lambda_2$ is in accordance with the previous remark ($T = 1, R = 0$). We also observe that the range of conversion wavelengths $\lambda_2$ by transmission is reduced at small $\lambda_1$ wavelengths.

Here, it is worth noting that the frequency $\omega_2$ is understood to arise in the system via the modulation of the metasurface properties in time. This could be analogous to an external pump in a pump-probe setup, where the metasurface undergoes dynamic changes influenced by a separate control beam or electromagnetic influence, effectively introducing $\omega_2$. Relevant temporal control schemes have been reported, e.g., based on the concept of thermoplasmonics [21] (in which the dielectric function of a thin film undergoes a modulation), further underscoring the feasibility of this approach.

So far we have not discussed the temporal form of the susceptibility, the relevance, and existence of which must be verified. Starting from (24)–(25), we can immediately write, after an inverse Fourier transformation:

$$ \tilde{\chi}_u(t, \omega_1) = \hat{\chi}_u^{11}(\omega_1, \omega_1) + e^{-i(\omega_2 - \omega_1)t} \hat{\chi}_u^{12}(\omega_1, \omega_2), \quad (53) $$

$$ \tilde{\chi}_u(t, \omega_2) = e^{i(\omega_2 - \omega_1)t} \hat{\chi}_u^{21}(\omega_2, \omega_1) + \hat{\chi}_u^{22}(\omega_2, \omega_2), \quad (54) $$

implying that the time shape of the Fourier transform of the susceptibility with respect to $t'$, i.e., $\tilde{\chi}_u(t, \omega)$, must be imposed at the two frequencies $\omega_i$ involved in the conversion. We note that this condition is only imposed at these two frequencies. In this context, we therefore need to check that this condition is realistic, i.e., compatible with causality. In fact $\chi(t, t')$ is causal versus $t'$, so that $\tilde{\chi}_u(t, \omega)$ is the Fourier transform of a causal function. In this context, we impose two values $\tilde{\chi}_u(t, \omega_i)$ on the Fourier transform of a causal function. We

<figure><img src="image_6.png" /><figcaption>FIG. 3. Normalized fluxes for reflection and transmission contoured against the incident wavelength $\lambda_1$ and conversion wavelength $\lambda_2$. The efficiency of the metasurface is represented by the color gradient, with yellow for high and blue for low efficiency. The diagonal where $\lambda_1 = \lambda_2$ marks the condition of no frequency conversion with maximum transmission ($T = 1$) and no reflection ($R = 0$), denoting the state of the metasurface where the incident and converted wavelengths are identical. The maps highlight how the conversion efficiency diminishes when $\lambda_2 > \lambda_1$ and reveal a restricted conversion range for shorter incident wavelengths.</figcaption></figure>
therefore need to ensure that the image space spanned by the Fourier transform of a causal function can include two arbitrary values. We can show that this is indeed the case, and that the result generalizes to a discrete set of frequencies $\omega_i$. It is not useful to give the general proof here, as we are now going to explicitly describe such a function in the case of two frequencies.

Indeed, until now the general form of the susceptibility has not been made explicit in the initial time space $(t, t')$, although this is essential for future numerical implementation. With reference to Appendix A, we obtain an expression of the form:

$$ \frac{\chi_u(t, t')}{A(t')} = f_u(t') + g_u(t') \cos[(\omega_2 - \omega_1)t] + q_u(t') \sin[(\omega_2 - \omega_1)t], \quad (55) $$

where $A(t')$ is a causal envelope given by:

$$ A(t') = \left(\frac{\xi}{2}\right) t'^2 e^{-\frac{t'}{\tau}} H(t'), \quad (56) $$

where, the 1/2 factor is for normalization, $\xi$ a constant scaling factor for the quadratic growth, $H(t')$ is the Heaviside distribution for the activation of the envelop, and $\tau$ is a characteristic decay time. Each function $v_u = f_u$, $g_u$ and $q_u$ (recall that $u = e$ or $h$) is defined with four parameters ($\lambda_{iu}^v, \mu_{iu}^v$) as:

$$ v_u(t') = \lambda_{1u}^v \cos \omega_1 t' + \mu_{1u}^v \sin \omega_1 t' + \lambda_{2u}^v \cos \omega_2 t' + \mu_{2u}^v \sin \omega_2 t'. \quad (57) $$

Finally, these parameters are given by the following matrix relationships:

$$ \begin{aligned} w_{2u}^v &= [G_2 - H_2(G_1)^{-1}H_1]^{-1} [2\beta_{2u}^v - 2H_2(G_1)^{-1}\beta_{1u}^v], \\ w_{1u}^v &= 2(G_1)^{-1}\beta_{1u}^v - (G_1)^{-1}H_1 w_{2u}^v, \end{aligned} \quad (58) $$

with: $w_{iu}^v = (\lambda_{iu}^v, \mu_{iu}^v)$, and $\beta_{iu}^v = [\Re(\alpha_{iu}^v), \Im(\alpha_{iu}^v)]$, and $\alpha_{iu}^v$ being the complex coefficients:

$$ \begin{aligned} (\alpha_{1u}^f, \alpha_{2u}^f) &= (\hat{\chi}_{u}^{11}, \hat{\chi}_{u}^{22}), \\ (\alpha_{1u}^g, \alpha_{2u}^g) &= (\hat{\chi}_{u}^{12}, \hat{\chi}_{u}^{21}), \\ (\alpha_{1u}^q, \alpha_{2u}^q) &= i(-\hat{\chi}_{u}^{12}, \hat{\chi}_{u}^{21}), \end{aligned} $$

while the four matrices $G_i$, $H_i$ are given in Appendix A. Thus, we have obtained the analytical expression for the process $\omega_1 \xrightarrow{\chi(t,t')} \omega_2$, that is, a susceptibility that gives rise to a frequency conversion in transmission.

In addressing the practical implications of our study, it is crucial to discuss the feasibility and practical aspects of temporal metasurfaces. The considerations are multifaceted, encompassing whether there is a need for a side band distant from the pump frequency, and how closely the modulating frequency should match the pump's order. Despite these factors, we do not encounter prohibitive theoretical constraints or practical limitations, even when the metasurface operates within high-frequency domains such as visible or near-infrared spectra. Indeed, the feasibility of relevant high-frequency modulation schemes has been substantiated in various photonic systems [23]. Thus, the scope of applicability for temporal metasurfaces, as advanced in our study,

remains viable and broad. In implementing the described conversion, our central assumption has meant that the response of the material at a particular moment is determined not only by the interactions up to that moment (the ordinary memory) but also by the way in which the material's properties have been altered up to that point. Metaphorically, it is as if the material's history is being rewritten in real time, which changes the instruction that the memory supplies about how the material should respond to current stimuli.

## VI. CONCLUSION

In conclusion, while in the classical sense, frequency conversion typically arises due to nonlinear interactions, we have shown here that linear-frequency conversion is feasible in a metasurface without resorting to nonlinear interactions in materials. By invoking only linear equations, we described how single-frequency conversion may be achieved in transmission through a metasurface, similar to parametric processes in time-varying systems. As captured by the introduced generalized time-inhomogeneous convolution product, the time variation itself via a susceptibility that features external memory effects, introduces an effective nonlinearity. The introduced conversion was shown to respect causality and energy conservation. Our theoretical prediction paves the way to the experimental realization of a powerful frequency conversion and control platform. The obtained solutions suggest that time-varying planar metasurfaces offer a viable path to mimicking nonlinear devices while providing optimal linear frequency conversion from one medium to another. This is remarkable, since in the proposed approach, no extraordinary assumptions were made for the type of materials or range of frequencies envisioned. All materials were chosen to be isotropic, homogeneous, and nonchiral. Within this framework, the generalized susceptibilities, represented by two-variable functions $\chi_u(t, t')$ and first-Fourier transform $\tilde{\chi}_u(t, \omega)$ rendered the system to be time and frequency dispersive. Such double dispersion furnished sufficient additional degrees of freedom to provide the opportunity for frequency conversion. The calculations indicate no limit on the number of arbitrary frequencies that can be converted allowing, albeit at the cost of added analytical complexity, for the process $(\omega_1, \omega_2, \dots) \rightarrow (\omega'_1, \omega'_2, \dots)$ by reflection or transmission. Under the assumption of perfect monochromaticity of the incident field, we also have full coherence among all the fields. As with any analytical prediction, numerical solutions are indispensable for validation and in further investigation. Tackling the derived equations for obtaining numerical solutions, beyond the simplified results shown in Figs. 2 and 3, is forthcoming. Extending our findings into the quantum regime is also being considered for future direction.

## ACKNOWLEDGMENTS

We acknowledge CNRS, Aix Marseille University, and Ecole Centrale for their constant support. A.P. acknowledges support from the Laboratory Directed Research and Development Program at Oak Ridge National Laboratory (ORNL). ORNL is managed by UT-Battelle, LLC, for the US Department of Energy under Contract No. DE-AC05-00OR22725.
therefore need to ensure that the image space spanned by the Fourier transform of a causal function can include two arbitrary values. We can show that this is indeed the case, and that the result generalizes to a discrete set of frequencies $\omega_i$. It is not useful to give the general proof here, as we are now going to explicitly describe such a function in the case of two frequencies.

Indeed, until now the general form of the susceptibility has not been made explicit in the initial time space $(t, t')$, although this is essential for future numerical implementation. With reference to Appendix A, we obtain an expression of the form:

$$ \frac{\chi_u(t, t')}{A(t')} = f_u(t') + g_u(t') \cos[(\omega_2 - \omega_1)t] + q_u(t') \sin[(\omega_2 - \omega_1)t], \quad (55) $$

where $A(t')$ is a causal envelope given by:

$$ A(t') = \left(\frac{\xi}{2}\right) t'^2 e^{-\frac{t'}{\tau}} H(t'), \quad (56) $$

where, the 1/2 factor is for normalization, $\xi$ a constant scaling factor for the quadratic growth, $H(t')$ is the Heaviside distribution for the activation of the envelop, and $\tau$ is a characteristic decay time. Each function $v_u = f_u$, $g_u$ and $q_u$ (recall that $u = e$ or $h$) is defined with four parameters ($\lambda_{iu}^v, \mu_{iu}^v$) as:

$$ v_u(t') = \lambda_{1u}^v \cos \omega_1 t' + \mu_{1u}^v \sin \omega_1 t' + \lambda_{2u}^v \cos \omega_2 t' + \mu_{2u}^v \sin \omega_2 t'. \quad (57) $$

Finally, these parameters are given by the following matrix relationships:

$$ \begin{aligned} w_{2u}^v &= [G_2 - H_2(G_1)^{-1}H_1]^{-1} [2\beta_{2u}^v - 2H_2(G_1)^{-1}\beta_{1u}^v], \\ w_{1u}^v &= 2(G_1)^{-1}\beta_{1u}^v - (G_1)^{-1}H_1 w_{2u}^v, \end{aligned} \quad (58) $$

with: $w_{iu}^v = (\lambda_{iu}^v, \mu_{iu}^v)$, and $\beta_{iu}^v = [\Re(\alpha_{iu}^v), \Im(\alpha_{iu}^v)]$, and $\alpha_{iu}^v$ being the complex coefficients:

$$ \begin{aligned} (\alpha_{1u}^f, \alpha_{2u}^f) &= (\hat{\chi}_{u}^{11}, \hat{\chi}_{u}^{22}), \\ (\alpha_{1u}^g, \alpha_{2u}^g) &= (\hat{\chi}_{u}^{12}, \hat{\chi}_{u}^{21}), \\ (\alpha_{1u}^q, \alpha_{2u}^q) &= i(-\hat{\chi}_{u}^{12}, \hat{\chi}_{u}^{21}), \end{aligned} $$

while the four matrices $G_i$, $H_i$ are given in Appendix A. Thus, we have obtained the analytical expression for the process $\omega_1 \xrightarrow{\chi(t,t')} \omega_2$, that is, a susceptibility that gives rise to a frequency conversion in transmission.

In addressing the practical implications of our study, it is crucial to discuss the feasibility and practical aspects of temporal metasurfaces. The considerations are multifaceted, encompassing whether there is a need for a side band distant from the pump frequency, and how closely the modulating frequency should match the pump's order. Despite these factors, we do not encounter prohibitive theoretical constraints or practical limitations, even when the metasurface operates within high-frequency domains such as visible or near-infrared spectra. Indeed, the feasibility of relevant high-frequency modulation schemes has been substantiated in various photonic systems [23]. Thus, the scope of applicability for temporal metasurfaces, as advanced in our study,

remains viable and broad. In implementing the described conversion, our central assumption has meant that the response of the material at a particular moment is determined not only by the interactions up to that moment (the ordinary memory) but also by the way in which the material's properties have been altered up to that point. Metaphorically, it is as if the material's history is being rewritten in real time, which changes the instruction that the memory supplies about how the material should respond to current stimuli.

## VI. CONCLUSION

In conclusion, while in the classical sense, frequency conversion typically arises due to nonlinear interactions, we have shown here that linear-frequency conversion is feasible in a metasurface without resorting to nonlinear interactions in materials. By invoking only linear equations, we described how single-frequency conversion may be achieved in transmission through a metasurface, similar to parametric processes in time-varying systems. As captured by the introduced generalized time-inhomogeneous convolution product, the time variation itself via a susceptibility that features external memory effects, introduces an effective nonlinearity. The introduced conversion was shown to respect causality and energy conservation. Our theoretical prediction paves the way to the experimental realization of a powerful frequency conversion and control platform. The obtained solutions suggest that time-varying planar metasurfaces offer a viable path to mimicking nonlinear devices while providing optimal linear frequency conversion from one medium to another. This is remarkable, since in the proposed approach, no extraordinary assumptions were made for the type of materials or range of frequencies envisioned. All materials were chosen to be isotropic, homogeneous, and nonchiral. Within this framework, the generalized susceptibilities, represented by two-variable functions $\chi_u(t, t')$ and first-Fourier transform $\tilde{\chi}_u(t, \omega)$ rendered the system to be time and frequency dispersive. Such double dispersion furnished sufficient additional degrees of freedom to provide the opportunity for frequency conversion. The calculations indicate no limit on the number of arbitrary frequencies that can be converted allowing, albeit at the cost of added analytical complexity, for the process $(\omega_1, \omega_2, \dots) \rightarrow (\omega'_1, \omega'_2, \dots)$ by reflection or transmission. Under the assumption of perfect monochromaticity of the incident field, we also have full coherence among all the fields. As with any analytical prediction, numerical solutions are indispensable for validation and in further investigation. Tackling the derived equations for obtaining numerical solutions, beyond the simplified results shown in Figs. 2 and 3, is forthcoming. Extending our findings into the quantum regime is also being considered for future direction.

## ACKNOWLEDGMENTS

We acknowledge CNRS, Aix Marseille University, and Ecole Centrale for their constant support. A.P. acknowledges support from the Laboratory Directed Research and Development Program at Oak Ridge National Laboratory (ORNL). ORNL is managed by UT-Battelle, LLC, for the US Department of Energy under Contract No. DE-AC05-00OR22725.
APPENDIX: TIME SHAPE OF SUSCEPTIBILITY AND PROOF OF CAUSALITY

In fact, $\chi(t, t')$ is causal versus $t'$, so that $\tilde{\chi}_u(t, \omega)$ is the Fourier transform of a causal function. In this context, conditions Eqs. (53) and (54) mean imposing two values $\tilde{\chi}_u(t, \omega_i)$ on the Fourier transform of a causal function. We therefore need to ensure that the image space spanned by a causal function can include two arbitrary values. We can show that this is indeed the case and that the result generalizes to a discrete set of frequencies $\omega_i$. It is not useful to give the general proof here, as we are in fact going to explicit such a function in the case of two frequencies. To that end, the relations given by Eqs. (53) and (54) are first rewritten as:

$$ \tilde{\chi}_u(t, \omega_1) = \hat{\chi}_u^{11} + \hat{\chi}_u^{12} \cos[(\omega_2 - \omega_1)t] - i\hat{\chi}_u^{12} \sin[(\omega_2 - \omega_1)t], \quad (A1) $$

$$ \tilde{\chi}_u(t, \omega_2) = \hat{\chi}_u^{22} + \hat{\chi}_u^{21} \cos[(\omega_2 - \omega_1)t] + i\hat{\chi}_u^{21} \sin[(\omega_2 - \omega_1)t]. \quad (A2) $$

Then we consider three real causal functions $f_u(t')$, $g_u(t')$, $q_u(t')$ that satisfy, with $u = e$, or $h$:

$$ \tilde{f}_u(\omega_1) = \hat{\chi}_u^{11} \quad \text{and} \quad \tilde{f}_u(\omega_2) = \hat{\chi}_u^{22}, \quad (A3) $$

$$ \tilde{g}_u(\omega_1) = \hat{\chi}_u^{12} \quad \text{and} \quad \tilde{g}_u(\omega_2) = \hat{\chi}_u^{21}, \quad (A4) $$

$$ \tilde{q}_u(\omega_1) = -i\hat{\chi}_u^{12} \quad \text{and} \quad \tilde{q}_u(\omega_2) = i\hat{\chi}_u^{21}. \quad (A5) $$

The next step is to build a causal susceptibility in the form:

$$ \chi_u(t, t') = f_u(t') + g_u(t') \cos[(\omega_2 - \omega_1)t] + q_u(t') \sin[(\omega_2 - \omega_1)t]. \quad (A6) $$

That is, after Fourier transform versus $t'$:

$$ \tilde{\chi}_u(t, \omega) = \tilde{f}_u(\omega) + \tilde{g}_u(\omega) \cos[(\omega_2 - \omega_1)t] + \tilde{q}_u(\omega) \sin[(\omega_2 - \omega_1)t]. \quad (A7) $$

At this step, one can check that Eqs. (A3), (A4), (A5) allow the susceptibility $\chi_u(t, t')$ to satisfy relations Eqs. (A1) and (A2), which can be summarized as:

$$ \tilde{v}_u(\omega_1) = \alpha_{1u}^{v} = \alpha_{1u}^{v} \quad \text{and} \quad \tilde{v}_u(\omega_2) = \alpha_{2u}^{v} = \alpha_{2u}^{v}, \quad (A8) $$

with $v_u = f_u$, $g_u$, or $q_u$, and $\alpha_{iu}^v$ being the complex coefficients:

$$ \alpha_{1u}^f = \hat{\chi}_u^{11} \quad \text{and} \quad \alpha_{2u}^f = \hat{\chi}_u^{22}, \quad (A9) $$

$$ \alpha_{1u}^g = \hat{\chi}_u^{12} \quad \text{and} \quad \alpha_{2u}^g = \hat{\chi}_u^{21}, \quad (A10) $$

$$ \alpha_{1u}^q = -i\hat{\chi}_u^{12} \quad \text{and} \quad \alpha_{2u}^q = i\hat{\chi}_u^{21}. \quad (A11) $$

These $\alpha_{1u}^v$ coefficients being imposed, additional degrees of freedom are required to satisfy Eqs. (A3), (A4), (A5) and lead to write each function $f_u$, $g_u$, and $q_u$ as:

$$ v_u(t') = A(t')[\lambda_{1u}^v \cos(\omega_1 t') + \mu_{1u}^v \sin(\omega_1 t') + \lambda_{2u}^v \cos(\omega_2 t') + \mu_{2u}^v \sin(\omega_2 t')], \quad (A12) $$

with $A(t')$ a common causal envelope, and $\lambda_{iu}^v$, $\mu_{iu}^v$ real parameters. Then it remains to relate, for each $u$ and $v$ value, the four parameters $(\lambda_{u,i=1,2}^v, \mu_{u,i=1,2}^v)$ to $(\alpha_{1u}^v, \alpha_{2u}^v)$. After tedious calculation, following Eqs. (A8)–(A12) we get:

$$ 2\alpha_{1u}^v = (\lambda_{1u}^v + i\mu_{1u}^v)\tilde{A}(0) + (\lambda_{1u}^v - i\mu_{1u}^v)\tilde{A}(2\omega_1) + (\lambda_{2u}^v + i\mu_{2u}^v)\tilde{A}(\omega_1 - \omega_2) + (\lambda_{2u}^v - i\mu_{2u}^v)\tilde{A}(\omega_1 + \omega_2), \quad (A13) $$

$$ 2\alpha_{2u}^v = (\lambda_{1u}^v + i\mu_{1u}^v)\tilde{A}(\omega_2 - \omega_1) + (\lambda_{1u}^v - i\mu_{1u}^v)A(\omega_2 + \omega_1) + (\lambda_{2u}^v + i\mu_{2u}^v)\tilde{A}(0) + (\lambda_{2u}^v - i\mu_{2u}^v)\tilde{A}(2\omega_2). \quad (A14) $$

In the complex plane, Eqs. (A13)–(A14) give a system of four equations with four unknowns, that we can write in a matrix form:

$$ G_1 w_{1u}^v + H_1 w_{2u}^v = 2\beta_{1u}^v, \quad (A15) $$

$$ H_2 w_{1u}^v + G_2 w_{2u}^v = 2\beta_{2u}^v, \quad (A16) $$

with the vectors:

$$ w_{iu}^v = (\lambda_{iu}^v, \mu_{iu}^v), \quad \text{and} \quad \beta_{iu}^v = [\Re(\alpha_{iu}^v), \Im(\alpha_{iu}^v)]. \quad (A17) $$

and the matrices:

$$ G_1 = \begin{bmatrix} \Re[\tilde{A}(0) + \tilde{A}(2\omega_1)] & -\Im[\tilde{A}(0) - \tilde{A}(2\omega_1)] \\ \Im[\tilde{A}(0) + \tilde{A}(2\omega_1)] & -\Im[\tilde{A}(\omega_1 - \omega_2) - \tilde{A}(\omega_1 + \omega_2)] \end{bmatrix}, \quad (A18) $$
APPENDIX: TIME SHAPE OF SUSCEPTIBILITY AND PROOF OF CAUSALITY

In fact, $\chi(t, t')$ is causal versus $t'$, so that $\tilde{\chi}_u(t, \omega)$ is the Fourier transform of a causal function. In this context, conditions Eqs. (53) and (54) mean imposing two values $\tilde{\chi}_u(t, \omega_i)$ on the Fourier transform of a causal function. We therefore need to ensure that the image space spanned by a causal function can include two arbitrary values. We can show that this is indeed the case and that the result generalizes to a discrete set of frequencies $\omega_i$. It is not useful to give the general proof here, as we are in fact going to explicit such a function in the case of two frequencies. To that end, the relations given by Eqs. (53) and (54) are first rewritten as:

$$ \tilde{\chi}_u(t, \omega_1) = \hat{\chi}_u^{11} + \hat{\chi}_u^{12} \cos[(\omega_2 - \omega_1)t] - i\hat{\chi}_u^{12} \sin[(\omega_2 - \omega_1)t], \quad (A1) $$

$$ \tilde{\chi}_u(t, \omega_2) = \hat{\chi}_u^{22} + \hat{\chi}_u^{21} \cos[(\omega_2 - \omega_1)t] + i\hat{\chi}_u^{21} \sin[(\omega_2 - \omega_1)t]. \quad (A2) $$

Then we consider three real causal functions $f_u(t')$, $g_u(t')$, $q_u(t')$ that satisfy, with $u = e$, or $h$:

$$ \tilde{f}_u(\omega_1) = \hat{\chi}_u^{11} \quad \text{and} \quad \tilde{f}_u(\omega_2) = \hat{\chi}_u^{22}, \quad (A3) $$

$$ \tilde{g}_u(\omega_1) = \hat{\chi}_u^{12} \quad \text{and} \quad \tilde{g}_u(\omega_2) = \hat{\chi}_u^{21}, \quad (A4) $$

$$ \tilde{q}_u(\omega_1) = -i\hat{\chi}_u^{12} \quad \text{and} \quad \tilde{q}_u(\omega_2) = i\hat{\chi}_u^{21}. \quad (A5) $$

The next step is to build a causal susceptibility in the form:

$$ \chi_u(t, t') = f_u(t') + g_u(t') \cos[(\omega_2 - \omega_1)t] + q_u(t') \sin[(\omega_2 - \omega_1)t]. \quad (A6) $$

That is, after Fourier transform versus $t'$:

$$ \tilde{\chi}_u(t, \omega) = \tilde{f}_u(\omega) + \tilde{g}_u(\omega) \cos[(\omega_2 - \omega_1)t] + \tilde{q}_u(\omega) \sin[(\omega_2 - \omega_1)t]. \quad (A7) $$

At this step, one can check that Eqs. (A3), (A4), (A5) allow the susceptibility $\chi_u(t, t')$ to satisfy relations Eqs. (A1) and (A2), which can be summarized as:

$$ \tilde{v}_u(\omega_1) = \alpha_{1u}^{v_u} = \alpha_{1u}^{v} \quad \text{and} \quad \tilde{v}_u(\omega_2) = \alpha_{2u}^{v_u} = \alpha_{2u}^{v}, \quad (A8) $$

with $v_u = f_u$, $g_u$, or $q_u$, and $\alpha_{iu}^v$ being the complex coefficients:

$$ \alpha_{1u}^f = \hat{\chi}_u^{11} \quad \text{and} \quad \alpha_{2u}^f = \hat{\chi}_u^{22}, \quad (A9) $$

$$ \alpha_{1u}^g = \hat{\chi}_u^{12} \quad \text{and} \quad \alpha_{2u}^g = \hat{\chi}_u^{21}, \quad (A10) $$

$$ \alpha_{1u}^q = -i\hat{\chi}_u^{12} \quad \text{and} \quad \alpha_{2u}^q = i\hat{\chi}_u^{21}. \quad (A11) $$

These $\alpha_{1u}^v$ coefficients being imposed, additional degrees of freedom are required to satisfy Eqs. (A3), (A4), (A5) and lead to write each function $f_u$, $g_u$, and $q_u$ as:

$$ v_u(t') = A(t')[\lambda_{1u}^v \cos(\omega_1 t') + \mu_{1u}^v \sin(\omega_1 t') + \lambda_{2u}^v \cos(\omega_2 t') + \mu_{2u}^v \sin(\omega_2 t')], \quad (A12) $$

with $A(t')$ a common causal envelope, and $\lambda_{iu}^v$, $\mu_{iu}^v$ real parameters. Then it remains to relate, for each $u$ and $v$ value, the four parameters $(\lambda_{u,i=1,2}^v, \mu_{u,i=1,2}^v)$ to $(\alpha_{1u}^v, \alpha_{2u}^v)$. After tedious calculation, following Eqs. (A8)–(A12) we get:

$$ 2\alpha_{1u}^v = (\lambda_{1u}^v + i\mu_{1u}^v)\tilde{A}(0) + (\lambda_{1u}^v - i\mu_{1u}^v)\tilde{A}(2\omega_1) + (\lambda_{2u}^v + i\mu_{2u}^v)\tilde{A}(\omega_1 - \omega_2) + (\lambda_{2u}^v - i\mu_{2u}^v)\tilde{A}(\omega_1 + \omega_2), \quad (A13) $$

$$ 2\alpha_{2u}^v = (\lambda_{1u}^v + i\mu_{1u}^v)\tilde{A}(\omega_2 - \omega_1) + (\lambda_{1u}^v - i\mu_{1u}^v)A(\omega_2 + \omega_1) + (\lambda_{2u}^v + i\mu_{2u}^v)\tilde{A}(0) + (\lambda_{2u}^v - i\mu_{2u}^v)\tilde{A}(2\omega_2). \quad (A14) $$

In the complex plane, Eqs. (A13)–(A14) give a system of four equations with four unknowns, that we can write in a matrix form:

$$ G_1 w_{1u}^v + H_1 w_{2u}^v = 2\beta_{1u}^v, \quad (A15) $$

$$ H_2 w_{1u}^v + G_2 w_{2u}^v = 2\beta_{2u}^v, \quad (A16) $$

with the vectors:

$$ w_{iu}^v = (\lambda_{iu}^v, \mu_{iu}^v), \quad \text{and} \quad \beta_{iu}^v = [\Re(\alpha_{iu}^v), \Im(\alpha_{iu}^v)]. \quad (A17) $$

and the matrices:

$$ G_1 = \begin{bmatrix} \Re[\tilde{A}(0) + \tilde{A}(2\omega_1)] & -\Im[\tilde{A}(0) - \tilde{A}(2\omega_1)] \\ \Im[\tilde{A}(0) + \tilde{A}(2\omega_1)] & -\Im[\tilde{A}(\omega_1 - \omega_2) - \tilde{A}(\omega_1 + \omega_2)] \end{bmatrix}, \quad (A18) $$
$$H_1 = \begin{bmatrix} \Re[\tilde{A}(\omega_1 - \omega_2) + \tilde{A}(\omega_1 + \omega_2)] & -\Im[\tilde{A}(\omega_1 - \omega_2) - \tilde{A}(\omega_1 + \omega_2)] \\ \Im[\tilde{A}(\omega_1 - \omega_2) + \tilde{A}(\omega_1 + \omega_2)] & \Re[\tilde{A}(\omega_1 - \omega_2) - \tilde{A}(\omega_1 + \omega_2)] \end{bmatrix}, \qquad (A19)$$

$$H_2 = \begin{bmatrix} \Re[\tilde{A}(\omega_2 - \omega_1) + \tilde{A}(\omega_2 + \omega_1)] & -\Im[\tilde{A}(\omega_2 - \omega_1) - \tilde{A}(\omega_2 + \omega_1)] \\ \Im[\tilde{A}(\omega_2 - \omega_1) + \tilde{A}(\omega_2 + \omega_1)] & \Re[\tilde{A}(\omega_2 - \omega_1) - \tilde{A}(\omega_2 + \omega_1)]\chi_e \end{bmatrix}, \qquad (A20)$$

$$G_2 = \begin{bmatrix} \Re[\tilde{A}(0) + \tilde{A}(2\omega_2)] & -\Im[\tilde{A}(0) - \tilde{A}(2\omega_2)] \\ \Im[\tilde{A}(0) + \tilde{A}(2\omega_2)] & \Re[\tilde{A}(0) - \tilde{A}(2\omega_2)] \end{bmatrix}. \qquad (A21)$$

Note in Eqs. (A18)–(A21) that $G_1$ is symmetric because of Hermiticity ($\Re[\tilde{A}(0)] = 0$), and that $H_2$ is the transposed matrix of $H_1$ since $\Re(\tilde{A})$ is even and $\Im(\tilde{A})$ is odd. Finally, the susceptibility is fully identified since the $(\lambda_{iu}^v, \mu_{iu}^v)$ coefficients are given by:

$$w_{2u}^v = \{G_2 - H_2(G_1)^{-1}H_1\}^{-1}\{2\beta_{2u}^v - 2H_2(G_1)^{-1}\beta_{1u}^v\}, \qquad (A22)$$

$$w_{1u}^v = 2(G_1)^{-1}\beta_{1u}^v - (G_1)^{-1}H_1 w_{2u}^v. \qquad (A23)$$

In Eqs. (A22), (A23) the matrices whose inverses are calculated are symmetrical. We therefore have the analytical expression for a susceptibility $\chi_{u=e,h}(t, t')$ giving rise to a frequency conversion from $\omega_1$ to $\omega_2$ in transmission. Note that the envelope can be chosen as:

$$A(t') = \left(\frac{\xi}{2}\right) t'^2 e^{-t'/\tau} H(t') \implies \tilde{A}(\omega) = \frac{\xi}{\left(\frac{1}{\tau} - i\omega\right)^3}. \qquad (A24)$$

[1] M. Kadic, G. W. Milton, M. van Hecke, and M. Wegener, 3d metamaterials, Nat. Rev. Phys. **1**, 198 (2019).

[2] N. Engheta, Four-dimensional optics using time-varying metamaterials, Science **379**, 1190 (2023).

[3] T. Stav, A. Faerman, E. Maguid, D. Oren, V. Kleiner, E. Hasman, and M. Segev, Quantum entanglement of the spin and orbital angular momentum of photons using metamaterials, Science **361**, 1101 (2018).

[4] N. Engheta, Circuits with light at nanoscales: Optical nanocircuits inspired by metamaterials, Science **317**, 1698 (2007).

[5] S. Kim, W. Shang, S. Moon, T. Pastega, E. Lee, and T. Luo, High-performance transparent radiative cooler designed by quantum computing, ACS Energy Lett. **7**, 4134 (2022).

[6] R. Tyumenev, J. Hammer, N. Joly, P. S. J. Russell, and D. Novoa, Tunable and state-preserving frequency conversion of single photons in hydrogen, Science **376**, 621 (2022).

[7] F. Kaiser, P. Vergyris, A. Martin, D. Aktas, M. P. D. Micheli, O. Alibart, and S. Tanzilli, Quantum optical frequency up-conversion for polarisation entangled qubits: Towards interconnected quantum information devices, Opt. Express **27**, 25603 (2019).

[8] L. Tetard, A. Passian, S. Eslami, N. Jalili, R. H. Farahi, and T. Thundat, Virtual resonance and frequency difference generation by van der waals interaction, Phys. Rev. Lett. **106**, 180801 (2011).

[9] A. Farokh Payam and A. Passian, Imaging beyond the surface region: Probing hidden materials via atomic force microscopy, Sci. Adv. **9**, eadg8292 (2023).

[10] S. K. O. Soman, A tutorial on fiber kerr nonlinearity effect and its compensation in optical communication systems, J. Opt. **23**, 123502 (2021).

[11] S. Taravati and G. V. Eleftheriades, Pure and linear frequency-conversion temporal metasurface, Phys. Rev. Appl. **15**, 064011 (2021).

[12] X. Guo, Y. Ding, Y. Duan, and X. Ni, Nonreciprocal metasurface with space-time phase modulation, Light Sci. Appl. **8**, 123 (2019).

[13] A. M. Shaltout, V. M. Shalaev, and M. L. Brongersma, Spatiotemporal light control with active metasurfaces, Science **364**, eaat3100 (2019).

[14] A. Shaltout, A. Kildishev, and V. Shalaev, Time-varying metasurfaces and lorentz non-reciprocity, Opt. Mater. Express **5**, 2459 (2015).

[15] H. B. Sedeh, M. M. Salary, and H. Mosallaei, Time-varying optical vortices enabled by time-modulated metasurfaces, Nanophoton. **9**, 2957 (2020).

[16] H. Barati Sedeh, M. M. Salary, and H. Mosallaei, Topological space-time photonic transitions in angular-momentum-biased metasurfaces, Adv. Opt. Mater. **8**, 2000075 (2020).

[17] M. M. Salary and H. Mosallaei, Time-modulated conducting oxide metasurfaces for adaptive multiple access optical communication, IEEE Trans. Antennas Propag. **68**, 1628 (2020).

[18] Z. Wu and A. Grbic, Serrodyne frequency translation using time-modulated metasurfaces, IEEE Trans. Antennas Propag. **68**, 1599 (2020).

[19] H. Barati Sedeh, M. M. Salary, and H. Mosallaei, Active multiple access secure communication enabled by graphene-based time-modulated metasurfaces, IEEE Trans. Antennas Propag. **70**, 664 (2022).

[20] C. Caloz and Z.-L. Deck-Léger, Spacetime metamaterials' Part I: General concepts, IEEE Trans. Antennas Propag. **68**, 1569 (2020).

[21] A. Passian, A. Lereu, E. Arakawa, A. Wig, T. Thundat, and T. Ferrell, Modulation of multiple photon energies by use of surface plasmons, Opt. Lett. **30**, 41 (2005).

[22] D. Niu, M. Zerrad, A. Lereu, A. Moreau, J. Lumeau, J. A. Zapien, A. Passian, V. Aubry, and C. Amra, Excitation of bloch surface waves in zero-admittance multilayers for high-sensitivity sensor applications, Phys. Rev. Appl. **13**, 054064 (2020).

[23] D. Pacifici, H. J. Lezec, and H. A. Atwater, All-optical modulation by plasmonic excitation of cdse quantum dots, Nat. Photon. **1**, 402 (2007).
$$H_1 = \begin{bmatrix} \Re[\tilde{A}(\omega_1 - \omega_2) + \tilde{A}(\omega_1 + \omega_2)] & -\Im[\tilde{A}(\omega_1 - \omega_2) - \tilde{A}(\omega_1 + \omega_2)] \\ \Im[\tilde{A}(\omega_1 - \omega_2) + \tilde{A}(\omega_1 + \omega_2)] & \Re[\tilde{A}(\omega_1 - \omega_2) - \tilde{A}(\omega_1 + \omega_2)] \end{bmatrix}, \quad (A19)$$

$$H_2 = \begin{bmatrix} \Re[\tilde{A}(\omega_2 - \omega_1) + \tilde{A}(\omega_2 + \omega_1)] & -\Im[\tilde{A}(\omega_2 - \omega_1) - \tilde{A}(\omega_2 + \omega_1)] \\ \Im[\tilde{A}(\omega_2 - \omega_1) + \tilde{A}(\omega_2 + \omega_1)] & \Re[\tilde{A}(\omega_2 - \omega_1) - \tilde{A}(\omega_2 + \omega_1)]\chi_e \end{bmatrix}, \quad (A20)$$

$$G_2 = \begin{bmatrix} \Re[\tilde{A}(0) + \tilde{A}(2\omega_2)] & -\Im[\tilde{A}(0) - \tilde{A}(2\omega_2)] \\ \Im[\tilde{A}(0) + \tilde{A}(2\omega_2)] & \Re[\tilde{A}(0) - \tilde{A}(2\omega_2)] \end{bmatrix}. \quad (A21)$$

Note in Eqs. (A18)–(A21) that $G_1$ is symmetric because of Hermiticity ($\Re[\tilde{A}(0)] = 0$), and that $H_2$ is the transposed matrix of $H_1$ since $\Re(\tilde{A})$ is even and $\Im(\tilde{A})$ is odd. Finally, the susceptibility is fully identified since the $(\lambda_{iu}^v, \mu_{iu}^v)$ coefficients are given by:

$$w_{2u}^v = \{G_2 - H_2(G_1)^{-1}H_1\}^{-1} \{2\beta_{2u}^v - 2H_2(G_1)^{-1}\beta_{1u}^v\}, \quad (A22)$$

$$w_{1u}^v = 2(G_1)^{-1}\beta_{1u}^v - (G_1)^{-1}H_1 w_{2u}^v. \quad (A23)$$

In Eqs. (A22), (A23) the matrices whose inverses are calculated are symmetrical. We therefore have the analytical expression for a susceptibility $\chi_{u=e,h}(t, t')$ giving rise to a frequency conversion from $\omega_1$ to $\omega_2$ in transmission. Note that the envelope can be chosen as:

$$A(t') = \left(\frac{\xi}{2}\right) t'^2 e^{-t'/\tau} H(t') \implies \tilde{A}(\omega) = \frac{\xi}{\left(\frac{1}{\tau} - i\omega\right)^3}. \quad (A24)$$

[1] M. Kadic, G. W. Milton, M. van Hecke, and M. Wegener, 3d metamaterials, Nat. Rev. Phys. **1**, 198 (2019).

[2] N. Engheta, Four-dimensional optics using time-varying metamaterials, Science **379**, 1190 (2023).

[3] T. Stav, A. Faerman, E. Maguid, D. Oren, V. Kleiner, E. Hasman, and M. Segev, Quantum entanglement of the spin and orbital angular momentum of photons using metamaterials, Science **361**, 1101 (2018).

[4] N. Engheta, Circuits with light at nanoscales: Optical nanocircuits inspired by metamaterials, Science **317**, 1698 (2007).

[5] S. Kim, W. Shang, S. Moon, T. Pastega, E. Lee, and T. Luo, High-performance transparent radiative cooler designed by quantum computing, ACS Energy Lett. **7**, 4134 (2022).

[6] R. Tyumenev, J. Hammer, N. Joly, P. S. J. Russell, and D. Novoa, Tunable and state-preserving frequency conversion of single photons in hydrogen, Science **376**, 621 (2022).

[7] F. Kaiser, P. Vergyris, A. Martin, D. Aktas, M. P. D. Micheli, O. Alibart, and S. Tanzilli, Quantum optical frequency up-conversion for polarisation entangled qubits: Towards interconnected quantum information devices, Opt. Express **27**, 25603 (2019).

[8] L. Tetard, A. Passian, S. Eslami, N. Jalili, R. H. Farahi, and T. Thundat, Virtual resonance and frequency difference generation by van der waals interaction, Phys. Rev. Lett. **106**, 180801 (2011).

[9] A. Farokh Payam and A. Passian, Imaging beyond the surface region: Probing hidden materials via atomic force microscopy, Sci. Adv. **9**, eadg8292 (2023).

[10] S. K. O. Soman, A tutorial on fiber kerr nonlinearity effect and its compensation in optical communication systems, J. Opt. **23**, 123502 (2021).

[11] S. Taravati and G. V. Eleftheriades, Pure and linear frequency-conversion temporal metasurface, Phys. Rev. Appl. **15**, 064011 (2021).

[12] X. Guo, Y. Ding, Y. Duan, and X. Ni, Nonreciprocal metasurface with space-time phase modulation, Light Sci. Appl. **8**, 123 (2019).

[13] A. M. Shaltout, V. M. Shalaev, and M. L. Brongersma, Spatiotemporal light control with active metasurfaces, Science **364**, eaat3100 (2019).

[14] A. Shaltout, A. Kildishev, and V. Shalaev, Time-varying metasurfaces and lorentz non-reciprocity, Opt. Mater. Express **5**, 2459 (2015).

[15] H. B. Sedeh, M. M. Salary, and H. Mosallaei, Time-varying optical vortices enabled by time-modulated metasurfaces, Nanophoton. **9**, 2957 (2020).

[16] H. Barati Sedeh, M. M. Salary, and H. Mosallaei, Topological space-time photonic transitions in angular-momentum-biased metasurfaces, Adv. Opt. Mater. **8**, 2000075 (2020).

[17] M. M. Salary and H. Mosallaei, Time-modulated conducting oxide metasurfaces for adaptive multiple access optical communication, IEEE Trans. Antennas Propag. **68**, 1628 (2020).

[18] Z. Wu and A. Grbic, Serrodyne frequency translation using time-modulated metasurfaces, IEEE Trans. Antennas Propag. **68**, 1599 (2020).

[19] H. Barati Sedeh, M. M. Salary, and H. Mosallaei, Active multiple access secure communication enabled by graphene-based time-modulated metasurfaces, IEEE Trans. Antennas Propag. **70**, 664 (2022).

[20] C. Caloz and Z.-L. Deck-Léger, Spacetime metamaterials' Part I: General concepts, IEEE Trans. Antennas Propag. **68**, 1569 (2020).

[21] A. Passian, A. Lereu, E. Arakawa, A. Wig, T. Thundat, and T. Ferrell, Modulation of multiple photon energies by use of surface plasmons, Opt. Lett. **30**, 41 (2005).

[22] D. Niu, M. Zerrad, A. Lereu, A. Moreau, J. Lumeau, J. A. Zapien, A. Passian, V. Aubry, and C. Amra, Excitation of bloch surface waves in zero-admittance multilayers for high-sensitivity sensor applications, Phys. Rev. Appl. **13**, 054064 (2020).

[23] D. Pacifici, H. J. Lezec, and H. A. Atwater, All-optical modulation by plasmonic excitation of cdse quantum dots, Nat. Photon. **1**, 402 (2007).