Switchable and Simultaneous Spatiotemporal Analog Computing

Optical wave-based computing has enabled the realization of real-time information processing in both space and time domains. In the past few years, analog computing has experienced rapid development but mostly for a single function. Motivated by parallel space-time computing and miniaturization, we show that reconfigurable graphene-based metasurfaces offer a promising path towards spatiotemporal computing with integrated functionalities by properly engineering both spatial- and temporal-frequency responses. This paper employs a tunable graphene-based metasurface to enable analog signal and image processing in both space and time by tuning the electrostatic bias. In the first part of the paper, we propose a switchable analog computing paradigm in which the proposed metasurface can switch among defined performances by selecting a proper external voltage for graphene monolayers. Spatial isotropic differentiation and edge detection in the spatial channel and first-order temporal differentiation and metasurface-based phaser with linear group-delay response in the temporal channel are demonstrated. In the second section of the paper, simultaneous and parallel spatiotemporal analog computing is demonstrated. The proposed metasurface processor has almost no static power consumption due to its floating-gate configuration. The spatial- and temporal-frequency transfer functions (TFs) are engineered by using a transmission line (TL) model, and the obtained results are validated with full-wave simulations. Our proposal will enable real-time parallel spatiotemporal analog signal and image processing.


Introduction
Digital signal processors are widely used to accomplish a large variety of computational tasks.Despite their flexibility, they come with several drawbacks, such as significant power consumption, restricted processing speed, and incompatibility with high frequency operation due to the technological limitations of current analog-to-digital converters 1 .The recent theoretical and manufacturing progress in the field of artificial photonic materials, e.g.photonic crystals or metamaterials, has inspired a return of the old paradigm of analog-based computing, by leveraging low loss structures that can process signals carried by light as it propagates through them.Optical analog signal processing could allow, in principle, for real-time, ultrafast, low energy consumption, and parallel processing 2,3 .
Recently proposed schemes for wave-based analog computing can be classified into two main categories, depending on whether they perform the processing operation in the spatial or temporal domain.
Spatial computing, on one hand, processes information encoded in the spatial dependency of an electromagnetic field, such as its intensity profile.Among the most popular techniques, the use of a 4-F correlator, which operates with two lenses performing Fourier transforms surrounding a spatial frequency filter, has found many applications in Fourier optics, and was recently miniaturized using metamaterials and metasurfaces.Another known method for spatial processing is the Transfer Function (TF) method 2 .
In the latter, the angular response of the system is designed in real space, such that it emulates a specific transfer function corresponding to a given linear mathematical operation.This method is potentially more compact as it avoids the use of components to take a Fourier transform and large propagation paths, leading to devices down or below the size of the wavelength, based e.g. on photonic crystal slabs 4 , plasmonic surfaces 5 , metasurfaces [6][7][8][9][10][11] , photonic spin Hall insulators 12 , inverse-designed metastructures 13 or topological wave insulators 14 .In optics, one-and two-dimensional image edge detection has been demonstrated using various platforms, such as optical metasurfaces [15][16][17][18] and surface plasmons 5 .On the other hand, temporal analog signal processing systems manipulate signals in the time domain with a dispersive structure, called a phaser 19 .Modern analog optical temporal computing dates back to the work of Pandian and Seraji 20 , in which the transient optical response of a fiber ring resonator is investigated and proposed for several applications, including optical pulse differentiation, integration and delay.In recent years, various architectures have been reported for temporal processing, such as differentiators [21][22][23] , integrators 23,24 , and equation solvers 25,26 .Despite all these advances, computing devices capable of performing both temporal and spatial operations, either sequentially or simultaneously, have been left largely unexplored.Spatiotemporal processors can exploit a much higher number of degrees of freedom, which could potentially be beneficial to analog computing systems, by enlarging the channel bandwidth and parallel operation capability.
For switching between spatial, temporal, or spatio-temporal operation modes, externally tunable electromagnetic properties, such as the amplitude and phase of the reflection, is ideal.This can be achieved using a tunable active metasurface 27,28 .Up to now, a variety of strategies have been proposed to design reconfigurable metasurfaces.Two-dimensional materials such as graphene, that exhibit a plethora of exceptional electromagnetic and photonic properties, have attracted tremendous attention as promising candidates for compact switchable devices.The relaxation time of the excited carriers in graphene is in the picosecond range, which is interesting for ultrafast wave manipulation 29 .The arbitrary control of graphene's complex surface conductivity can also be continuously tuned by manipulating its Fermi level by electric gating or photo-induced doping, which directly provides efficient real-time control of reflected/transmitted waves 30 .Accordingly, Fermi level control through external biasing or chemical doping has enabled the integration of THz devices with flexible substrates 31 .The strong interaction of graphene with electromagnetic fields has been leveraged in impressive applications such as wideband absorbers [32][33][34] , polarization rotators 35 , near field imagers 36 , nanoantennas 37 , biosensors 38 , and THz wave devices [39][40][41][42] .
In this work, we propose an appealing power-efficient opportunity to perform analog signal and image processing in both time and space domains without resorting to intricate solutions involving various devices and bulky Fourier lenses.Systematically speaking, here, we consider a 2 × 2 multiple-input and multiple-output (MIMO) graphene metasurface processor in which the two inputs/outputs correspond to temporal and spatial processing channels.We use the tunability of graphene to construct a planar dual space-time processor that can perform space and time-domain signal processing tasks sequentially or at the same time.The metasurface processor can accomplish an isotropic spatial differentiation operation for edge detection of a spatially encoded image.Simultaneously, in the temporal channel, the metasurface processor can serve as differentiator, or as a metasurface-based phaser for temporal pulse spreading when the input signal has a time variation.We confirm that the graphene-based metasurface can be tuned to enable analog signal processing in both space and time domains by changing the external potential bias.
In such a design, the graphene monolayer can capture the electrons tunneling from a positive external bias, but cannot release them after releasing the DC voltage because it is electrically isolated from the biasing electrodes.Hence, no extra power is required to maintain the graphene's conductivity at the desired value 43 .We realize the required amplitudes and phases of the transfer function (TF) associated with spatial and temporal operations by combining the surface impedance model of a graphene monolayer to a transmission line (TL) approach.Several proof-of-principle examples are presented to illustrate diverse wave-based signal processing functionalities like edge detection, complex spatial signal processing, temporal differentiation, and pulse spreading.

4/28
First, we theoretically investigate the spatial and temporal optical transfer functions (OTFs) for different operations under normal and oblique incident waves.The OTF Õ(k x , k y , Ω) ( Ω = ωω 0 , ω 0 is the central frequency) of a metasurface is the complex function that maps the incident electric field to the reflected/transmitted field.Let us consider a uniform metasurface located in the x -y plane (z = 0).The incident and reflected fields, expressed in time domain, are denoted as ψ in (x, y, t) and ψ out (x, y, t) in a laboratory frame, respectively, where x and y are defined in the coordinate system of the beam, as shown in Figure 1.The OTF is written as a matrix to account for changes in the polarization of the reflected/transmitted light field.
We consider two channels: the spatial channel and the temporal channel (see Figure 1).These channels are separated by distinct polarizations or illumination angles; therefore, we can consider two OTFs, Õs (k x , k y ) and Õt (Ω), corresponding to the two channels, in spatial and temporal domains, respectively.
In the paraxial regime, the incident beam has the form P in ψ in (x, y), where a 2-vector P in in the x-y plane describes the input polarization, and ψ in (x y) describes the scalar electric field distribution on the plane perpendicular to the beam propagation direction.In the spatial Fourier domain, each beam profile can be spectrally represented by a superposition of monochromatic plane waves, by Fourier transform: ψ in (x, y) = ψin k x , k y exp (jk x x + jk y y) dk x dk y .
In addition, a similar equation can be considered for the temporal domain, as below, Due to the tangential wavevector's continuity along with the interface, the incident spatial frequency component with (k x , k y ) in the incident plane only generates the output spatial frequency component with the same (k x , k y ) in the output plane.At each (k x , k y ), the output wave has an electric field of ψout p (k x , k y ) = Õs (k x , k y ).P in ψin (k x , k y ),

5/28
where Õs (k x , k y ) is 2×2 matrix (See Appendix A).The output beam passes through a polarizer selecting an output polarization P out , thereby an output electric field ψ out p = P out ψ out (x, y); here, ψ out , similar to ψ in , is the field distribution on the output plane and has a spatial Fourier transform ψout (k x , k y ) in the (k x , k y ) domain.Therefore, the relation between ψout and ψin is ψout (k x , k y ) = Õs (k x , k y ) ψin (k x , k y ) where Õs (k x , k y ) = P out † O s (k x , k y )P in .Our desired OTF is the one corresponding to isotropic differentiation operation because it is one of the most fundamental operations in mathematics, and has several applications in engineering and image processing 6 .The transfer function of isotropic differentiation is Õs (k x , k y ) = k 2 x + k 2 y , that must have the obvious property Õs (k x = 0, k y = 0) = 0.In our platform, one can achieve Õs (k x = 0, k y = 0) = 0 by choosing the proper input and output polarizations such that When the above equation is satisfied, the OTF in the vicinity of k x = k y = 0 has the below form (See Appendix A), In order to achieve two-dimensional homogeneous differentiation, the transfer function must have a rotationally invariant magnitude, and therefore γ 2 /γ 1 = ±i.The zeros of Õs (k x , k y ) carry topological charge ±1.It means that after passing the Gaussian beam through the spatial differentiation system, the Gaussian beam will form a vortex light beam 44 .If the reflection coefficient for p-or s-polarized wave (r p0 and r s0 ) is zero, the Equation 4 is satisfied and by setting β = pi/2, γ 2 and γ 1 can be calculated by (See Appendix A) 6/28 where P in = P in x , P in y T .Therefore, the output electric field distribution ψ out (x, y) will be proportional to , as desired.For the temporal channel, let us propose two distinct temporal processing: (i) performing the first-order temporal differentiation operation; and (ii), achieving temporal pulse spreading via a synthesized linear group-delay response at normal illumination (see Figure 1).In this channel, the input beam directly illuminates at normal illumination without using rotation in the output polarization.Therefore, the OTF for the temporal channel can be defined as For temporal differentiator, the OTF can be represented as Therefore, it is clear that the envelope of the temporally reflected/transmitted field with central frequency ω 0 has the field profile of where α is a constant value.
In any temporal analog signal processor, there is a phaser, i.e., a two-port component with a transfer function exhibiting a group delay versus frequency response, which may be designed to show the groupdelay (e.g., linear, quadratic, cubic, stepped, etc.) as core component 45 .One of the crucial application of phaser is temporal pulse spreading via designing the linear group-delay response.The pulse-spreading operation allows us to steer the amplitude envelopes of quasi-sinusoid EM waves as desired, which is one of the basic impacts of temporal analog computing.In fact, input signal traveling along such a phaser experiences time spreading since its different spectral components travel with different group velocities, they temporally rearranged 46 .By exploiting this temporal rearrangement, the various spectral components of a signal can be directly mapped onto the time domain and can then be processed in a real-time manner.
where a is the group-delay slope and b is constant.Consider the incident modulated Gauss pulse of duration T and bandwidth B, with central frequency ω 0 .As different spectral components of the Gauss pulse have different group delays when propagating through this phaser, the incident EM pulse spreads over the time sequence, resulting in a broader reflected pulse with a duration of 19 where ∆τ is the group-delay swing over the frequency band B, and C represents the spreading factor of the spatial phaser.Additionally, the peak power of the reflected EM pulse is diminished to P 0 /C.The schematic view of the proposed reconfigurable metasurface processor with the mentioned functionality in 8/28 the spatial and temporal domain is illustrated in Figure 1.

Graphene-Based Metasurface Design
, where k 0 is the free space propagation constant, ε d is permittivity of dielectric, and θ is the incident angle.The TL's characteristic impedance can be presented by Z s and Z p for the incidence wave polarized with s and p polarizations, respectively.We define characteristic impedances as Z s = η 0 / ε d -sin 2 θ , and Z p = η 0 ε d -sin 2 θ /ε d for s and p polarizations, respectively, where η 0 indicates free space impedance.According to the utilized method, the TL matrix is given by where h shows TL length and Z d has to be replaced with the characteristic impedance of the TL with the considered polarization.In addition, a shunt admittance Y is introduced by the matrix below: Then, the equivalent circuit matrix of the metasurface can be represented by Finally, the scattering matrix is calculated by where Z 0 is the characteristics impedance of free space for the specific polarization considered, namely Z 0,s = η 0 /cosθ for s-polarization and Z 0,p = η 0 cosθ for p-polarization.
In order to implement the system, we propose to use a reflective graphene-based metasurface capable of switching between different states, depending on the illuminating beam.Graphene complex conductivity 10/28 According to graphene's outstanding features, we have utilized it as a tunable platform in a metasurface processor to switch between predetermined computational states.To start the design and perform an accurate evaluation of the proposed structure, graphene is modeled as an infinitesimally thin sheet with surface impedance Z = 1/σ g , where σ g is the frequency-dependent complex conductivity of graphene.The surface conductivity of graphene including both intraband (σ intra ) and interband (σ inter ) transitions are governed by the well-known Kubo formula 48,49 σ g (ω, τ, µ c , T) = σ intra (ω, τ, µ c , T) + σ inter (ω, τ, µ c , T) , (17)   11/28 σ intra (ω, τ, µ c , T) = -j e 2 k B T πh 2 (ω-jτ -1 ) where e, h, and k B are constants corresponding to electron charge, the reduced Planck's constant, and the Boltzmann constant, respectively 48 .In the above equation, variables T, τ, and µ c correspond to the environmental temperature, relaxation time, and the chemical potential of graphene, and ω is the angular frequency 48 .Note that the above expression neglects the graphene's edge effects and considers that the Drude-like intraband contribution dominates, which are experimentally confirmed assumptions in the considered frequency range 50 .In explaining graphene's optical response, the conical band diagram is essential for defining light graphene interaction dynamics.Two types of band transitions are possible when a photon interacts with the graphene surface 51 .Depending on the relative positions of the Fermi level (E f ) and the energy of the incident photons, light absorption is either dominated by interband or intraband transitions, and the effects of these transitions are determined by Pauli blocking 52 .When the incident photon energy is lower than 2E f , intraband transitions become dominant, whereas in the opposite case, interband transitions dominate 53 .The interband conductivity is on the order of e 2 /h, and at frequencies below the THz regime and room temperatures, the interband term in complex conductivity is very small compared to the intraband term and usually can be ignored 48 .Moreover, graphene can be modeled as a thin dielectric.In this case the dielectric permittivity ε g is expressed as 54,55 where ∆ g is the graphene thickness, and ε 0 is the vacuum permittivity.In a simple model, we neglect the other parts of the floating-gate structure except for the graphene layer because their thickness is much smaller than a wavelength and their relative permittivities are similar to the ones of the nearby substrate.

12/28
Hence, the floating gate is modeled as a shunt impedance, and its value is equal to 1/σ g .In such a structure, the unpatterned graphene layer is a lossy medium that is represented through a series RL circuit in the demonstrated TL model.The frequency-dependent resistance and inductance are calculated by By optimizing the TLs and graphene layer parameters, we can achieve the required behavior of the structure.In this work, the environment temperature is considered to be T = 300 K.The proposed design consists of a fully covered graphene layer on top of a silicon substrate with a thickness of h and a metallic ground plane on the backside.In this composition, the relative permittivity and loss tangent of the silicon substrate are ε r = 11.9 and tan δ = 0.00025.The same structure consisting of a silicon substrate and an unpatterned graphene layer with different electrostatic bias is embedded upon the primary segment.We assumed the same TLs in both sections, but we could use different TLs with different materials and lengths if more degree of freedom were required by the target functionality.Indeed, by using the same TLs, we decreased the number of parameters to optimize and reduced the computational load of the optimization process.We numerically simulate the design in CST Studio Suite to get the reflection amplitude and phase.
The amplitude and phase control in the graphene metasurface is achieved via changes in its complex conductivity controlled by an external voltage on the floating gates, which changes the chemical potential.
In most prior research works, a simple capacitive structure consisting of graphene, an insulator, and a metallic electrode is used to tune graphene's conductivity.However, the disadvantage of this method is that an external driven voltage must be applied continually to sustain graphene's conductivity.Consequently, static power consumption is inevitable.In order to resolve this drawback, we propose a metasurface with nearly zero static power consumption based on a non-volatile floating-gate graphene structure as widely used in some non-volatile devices [56][57][58] .This structure consists of Si, SiO2, Al2O3, and graphene 13/28 monolayers, as illustrated in 43 .When a positive voltage is applied to the top Si layer, the electrons in the bottom Si layer can tunnel through the SiO2 and be captured by the graphene.Hence, the charge density and chemical potential of the graphene layer will be increased.On the contrary, when the reverse voltage is applied to the Si layer, the reverse electric field intensity in SiO2 and reverse tunneling current increase with the reverse voltage.Accordingly, the charge density and chemical potential of graphene are decreased 43 .The electrons can only tunnel through the SiO2 layer while they cannot tunnel through the Al2O3 layer because its thickness is much larger than that of SiO2 43 .By employing the proposed configuration, after removing the bias voltages, the graphene's charge density can remain unchanged for a long time since the graphene is electrically isolated from the Si layers.This means that no more power is needed to keep the graphene's charge density constant.
The selected biasing configuration for the metasurface is depicted schematically in Figure 2(b).In the graphene unpatterned layers, µ c can be tuned by adjusting the charge density of graphene n as given by 59 where v F = 10 6 m/sec is the Fermi velocity 60 .We employ a floating-gate structure to tune the charge density of the graphene.The floating-gate design is widely used in some non-volatile devices 61 .The proposed design is composed of Si, SiO2, Al2O3, and unpatterned graphene, as shown in Figure 2(b).
In the floating-gate design, when a voltage is applied on the top Si layer, the electrons in the bottom Si layer can tunnel through the SiO2 and been captured by the graphene, increasing its charge density.In this case, we assumed the voltage on the top Si layer is positive.On the contrary, when the reverse voltage is applied on the top Si layer, the electrons in graphene can tunnel through the SiO2 and been captured by the bottom Si layer, so the charge density is decreased 43 .In this design, the graphene monolayer is isolated electrically from the Si layers, which means after removing the voltage, the charge density of graphene can remain constant for an extended period of time.So, in this structure, no extra power is 15/28 needed to keep charge density consistent.The thickness of the Al2O3 and SiO2 are h Al2O3 = 20 nm, and h SiO2 = 10 nm, respectively.Therefore, the electrons can only tunnel through the SiO2 layer whereas cannot tunnel through the Al2O3 layer because its thickness is much larger compared to the SiO2 layer.
Also, the relative permittivity of the Al2O3 layer is ε Al2O3 = 9, and the relative permittivity of SiO2 is ε SiO2 = 3.9.The charge density of graphene is expressed as the integral of tunneling current J SiO2 in the SiO2 which is given by where t 0 is the duration of the voltage applied on the Si.The tunneling current according to the Fowler-Nordheim tunneling mechanism is given by 62 In the above equation, φ SiO 2 = 3.2 eV is the barrier height of SiO2, E SiO2 shows the electric field in SiO2, and m indicates effective mass of electron.The electric field in SiO2 is expressed by where the applied voltage to graphene layer is indicated by V g .Using the Kubo formula and voltagedependent property of floating-gate, the required external voltage for a specific value of chemical potential can be calculated.When the voltage on the top Si layer is positive, the value of σ g increases by voltage.
When the voltage is negative, the value of σ g decreases with the reverse voltage.Finally, after removing the voltage the, value of σ g remains constant, and no voltage is needed to sustain σ g .Hence, the floating-gate tuning design is preferred to the capacitor-based tuning methods, where the static power consumption is inevitable 63,64 .
As we discussed in the above section, we modeled reconfigurable metasurface by leveraging a TL approach to synthesize desired OTFs in spatial and temporal domains.In fact, we engineered the temporal-and spatial-frequency response of the proposed metasurface by fine-tuning the scattering parameters of Equation 16.We defined a cost function, F = ||OTF des (k x /k y , Ω) -OTF syn (k x /k y , Ω)||, and searched for a set of circuit parameters which ensures F → 0, where OTF des is our desired OTF and OTF syn is S 11 in Equation 16.A particle swarm optimization (PSO) algorithm 65 is adopted to minimize the value of F.

Switchable Analog Computing Performance
This section wants to investigate a switchable metasurface that can switch between determined operations in the desired channel (temporal or spatial).In this scenario, the relaxation time of graphene is assumed to be τ = 0.4 ps, and the dielectric thickness is h = 12 µm.The required chemical potential for graphene monolayers to achieve the desired functionality-isotropic spatial differentiation, temporal differentiation and linear group-delay response-is illustrated in Figure 2(d).We numerically simulate the design in CST Studio Suite to get the reflection amplitude and phase.The obtained results show that the graphene-based metasurface supports both spatial and temporal channels, which can be selected by applying a proper external electrostatic voltage.To explicitly present the properties of the synthesized metasurface processor, the spatial and temporal OTFs as a function of the incidence angle and temporal frequency are illustrated in  5.The amplitude and phase of OTF for (k x , 0) and (0, k y ) are plotted in Figure 3(c).This is evidence that two-dimensional edge detection can be achieved properly.The right side of Figure 3 shows the temporal OTFs associated with first-order differentiation operation and linear group delay response.Figure 3(d) exhibits the amplitude and phase of the synthesized OTF associated with first-order temporal differentiation.For ease of comparison, the ideal We now move to the phaser, which provides a desired group-delay response over a given frequency band 19 .For an incident EM wave that has various sinusoidal components, it has the form A n (t).e jω n t+φ n .( 27) Thus, the output EM wave has the form 66 From the above equation, we can observe that the amplitude envelopes of different sinusoidal components have different time delays depending on their frequencies.We now switch to the spatial channel at an oblique incident angle, at which the graphene-based metasurface performs isotropic spatial differentiation.We first consider the scenario where the input is a Gaussian beam, as shown in Figures 6(a (f).For the sake of comparison, we also plot the ideal normalized magnitude of a differentiated Gaussian

Simultaneous Spatiotemporal Analog Computing Performance
Here, we investigate another opportunity, namely to perform parallel analog spatiotemporal computing, with the proposed graphene-based metasurface.Our motivation is the potential enhancement of channel bandwidth by expanding the analog computing from the single spatial or temporal operation to parallel spatiotemporal operation 67,68 .
Consider the scenario in which temporal and spatial channels are simultaneously excited by spatial and temporal signals at different angles as well as frequencies.Similar to the previous section, we consider first-order spatial differentiation as a spatial operation and first-order temporal differentiation as well as linear group-delay response as temporal functions.In this scenario, graphene's relaxation time is considered as τ = 0.06 ps, and the substrates thickness are h = 18µm.The required chemical potential for graphene monolayers to achieve the desired functionality is µ c,1 = µ c,2 = 0.32 eV. Figure 8 provides the results of parallel computing including spatial and temporal transfer functions.Figures 8(a A complex signal which is a combination of Gaussian and sinusoidal functions (see Figure 9

Conclusion
To conclude, we exploited a novel switchable subwavelength architecture for optical analog computing of spatiotemporal and simultaneous computing.The proposed metasurface is compact without need of optical Fourier transform elements, where optical computation functions are directly achieved in the real space rather than the Fourier space.By leveraging the surface impedance model of graphene layers and using the TL approach, the computational metasurface for temporal and spatial processing with a proper spatial and temporal bandwidth was engineered.In the temporal channel, first-order differentiation and metasurface-based phaser for temporal pulse spreading application via synthesizing linear group-delay response are demonstrated.Besides, in the spatial channel, performing isotropic differentiation operation as well as edge detection are provided.According to the floating gate design utilized in this paper, the 21/28 external voltage can be removed after biasing the graphene layers according to the desired function.We have eventually validated the spatial and temporal differentiation with one-dimensional spatial/temporal signals and an image, with excellent agreement between full-wave numerical simulations and targetted operations.
Acknowledgements: A. Momeni and R. Fleury acknowledge funding from the Swiss National Science Foundation under the Eccellenza grant number 181232.

A Appendix: Calculation of Spatial Transfer Function
We define P in = (P in x , P in y ) T and P out = (P out x , P out y ) T as the normalized input and output polarizations in x-y plane, respectively.
The transfer function Õs (k x , k y ) = P out † O s (k x , k y )P in can be written as 44,69 Õs (k x , k y ) = where where the matrices V 1 and V 2 are originated from the rotations of coordinates and ζ = cot(θ 0 )k y /k 0 .Also, the r p and r s are the Fresnel coefficients for p-and s-polarized plane waves, respectively.According to Equation 4 in main text, and Equation 29 in (k x = 0, k y = 0) we have -r p0 P in x P out x * + r s0 P in y P out y * = 0, (31)

Figure 1 .
Figure 1.Schematic of the proposed spatiotemporal metasurface processor.We can switch between different functionality states by controlling an external voltage via an FPGA processor.

Figure 2 .
Figure 2. (a) Proposed TL model to achieve the desired reflection.(b) Schematic of the reconfigurable device and required layers for biasing the graphene monolayers.(c) Equivalent circuit model for the graphene-based computational metasurface.(d) Table of the required chemical potential of each graphene monolayer for a switchable scenario.
Metasurfaces are conventionally characterized by surface-averaged material parameters, i.e., polarizabilities, susceptibilities, or surface impedances.In this regard, we choose the surface impedance model, which represents a metasurface as an equivalent circuit model of specific configuration for realizing desired OTFs.The simple form of the proposed circuit model is illustrated in Figure2(a).It contains three shunt admittances Y 1 , Y 2 , and Y 3 and two TLs of arbitrary lengths.In this model, the propagation constant of the guided mode along the TL is k

Figure 3 .
Figure 3. Synthesized Spatial and temporal OTF for different operations in the spatial (a,b) and temporal (d-f) channels.(a) Amplitude and (b) phase distribution of the spatial transfer function.(c) Comparison of amplitude and phase distribution in the (k x , 0) and (0, k y ) planes.(d) Amplitude and phase of the OTF of the temporal differentiator.(e) Synthesized positive-triangular group delay and (f) amplitude and phase responses of the metasurface-based phaser.The synthesized and ideal OTFs are indicated with solid and dashed lines, respectively.Black and green data points are generated independently, using the circuit model of the graphene-based metasurface.

Figure 4 .
Figure 4. (a) The temporal envelope of the incident field is the input signal.(b) The ideal output we target, with its envelope representing the time derivative of the input signal.(c) Actual output, in very close agreement with the ideal output.

Figure 5 .
Figure 5. Temporal pulse spreading using a graphene-based metasurface processor.(a) Illustration of pulse spreading based on our metasurface-based phaser with linear group-delay response.(b) Temporal envelope of the input signal and (c) temporal envelope of the output signal, scaled horizontally and vertically by factors C and K, respectively.

Figure 6 .
Figure 6.Spatial differentiation on an input Gaussian beam and generation of a vortex beam.Intensity distribution of the (a) incident and (b) reflected fields.(c) Intensity distribution of the incident field along the dashed lines in (a).(d) Intensity distribution of the output field along the dashed lines in (b).The ideal and output signals have also been illustrated for the sake of comparison in (d).(e) and (f) The real part and phase of reflected fields, respectively.

Figure 3 .
Figure 3.The obtained results are exactly what we explained in Section 2 and the theoretical calculations, which are calculated based on the TL approach in Section 3 are in excellent agreement with simulated results.The left side of Figure 3 represents the spatial OTF for isotropic differentiation operation.The input polarization can be computed by Equations 6 and 7 and the fact that γ 2 /γ 1 = i.The operating frequency for the spatial channel is 1.58 THz and θ 0 = 47.5 • .Based on Equation 5 , we have the helical phase distribution of the spatial transfer function in Figure 3(b) as we expected from Equation 5.The amplitude and phase of

Figure 7 .
Figure 7. Edge detection of a 2-dimensional image by exploiting the proposed graphene-based metasurface processor.(a) Input and (b) the edge-detected images, respectively.

Figure 8 .
Figure 8. Synthesized Spatial and temporal OTF for parallel spatiotemporal computing.(a) Amplitude (blue) and phase (red) of the spatial transfer function associated with the spatial differentiation operation.(b) Amplitude and phase of the temporal transfer function associated with the temporal differentiation operation.(c) Amplitude and phase distribution of synthesized metasurface for linear group-dealy response.(d) Synthesized positive-triangular group delay.The synthesized and ideal OTFs are indicated with solid and dashed lines, respectively.Also, green and black cross lines are related to the results of TL model of the graphene-based metasurface.

Figure 5
shows the process to realize pulse spreading through a linear group-delay response metasurface-based phaser.An incident modulated Gauss pulse of duration T (the time range of 10% peak amplitude) and bandwidth B, with central frequency ω 0 , is radiated to the metasurface-based phaser.As we can see in Figure3(e), the group delay has a positive slope (approximatly 1.28 ps 2 /rad) from 0.7 to 1 THz, and a negative slope from 1 to 1.3 THz.Based on Equation12, the predicted spreading time is ∆τ = aB = 2.5 ps.As we can see in Figures5(b) and (c), the synthesized linear group delay leads to 2.5 ps spreading time duration of input temporal signal.The values of C and K = 1/ √ C are 1.25 and 0.89, respectively.
) and (c).The operator in Equation5, operating on an input Gaussian beam, produces a beam with a vortex-shaped beam that possesses a donut-shaped intensity profile (see Figure6(b)).The real part and phase of reflected fields are demonstrated in Figures6(e) and

Figure 9 .
Figure 9. (a) The complex input field profile is the combination of Gaussian and sinusoidal functions.We compare (b) the output signal envelope to (c) the exact derivative of the input signal envelope.The signals are normalized to one.(d) A Gaussian pulse is used as an input in the time domain.(e) and (f) The normalized output and ideal temporal signals, respectively.
) and(b)   show the spaial and temporal OTFs associated with first-order differentiation operation, respectively.The operating frequency for the temporal channel is 1.9 THz and θ 0 = 34.6 • .Also, the center frequency of the spatial channel is 0.65 THz and θ 0 = 34.6 • .Similar to the previous section, the ideal case and TL method results are also plotted in same figures.The bandwidth for spatial and temporal channels are |k/k 0 | < 0.2 and |Ω/ω 0 | < 0.1.Figure8(c) exhibits the amplitude and phase of synthesized OTF associated with linear group delay response of Figure8(d).From Figure8(d), the group delay has a positive slope from 6.45 to 6.7 THz and vise-versa a negative slop from 6.7 to 6.95 THz.
(a)), has been utilized as the beam profile of incidence.The corresponding results, including output signal and an ideal case, are illustrated in Figures9(b) and (c).In addition, we adopt a Gaussian pulse for temporal differentiation as the input signal (see Figure9(d)).The output signal and an ideal case, are demonstrated in Figures9(e) and (f).As can be seen, the metasurface successfully implements the spatio-temporal differentiation of the input signal.