Analytical formulation of spatiotemporal modulated graphene-based waveguides using Floquet-Bloch theory

In this work, an analytical model to study graphene-based spatiotemporal modulated structures is developed and verified through comparison with full wave numerical simulations. Graphene is an ideal material for realizing spatiotemporal modulated structures at high frequencies of THz and optics. In this analysis, the electromagnetic response of studied structures is expressed in terms of weighted Floquet-Bloch modes supported by the structure, while graphene is modeled by a spatiotemporal modulated surface current that imposes certain boundary conditions on the modes. The developed analytical technique is a comprehensive tool and can be used for accurate modeling of different kinds of spatiotemporal devices including lossy, guided, and leaky wave structures. To demonstrate the accuracy of the model, two plasmonic waveguides with space and time modulated graphene conductivity are analyzed and their interband and intraband transition between modes are thoroughly investigated. Using the developed analytical model, spatiotemporal modulation phenomena such as mode conversion, wave amplification and nonreciprocal response are explored and discussed for the studied structures.

www.nature.com/scientificreports/three fundamental harmonics are considered in lossless transmission line equations, which limits this method to be used in some applications which require to consider more harmonics 47,48 .To overcome the inherent limitations of CMT, Floquet-Bloch theory was proposed for analyzing structures with spatiotemporal modulated constitutive parameters (permeability and permittivity), which provides complete answers 6,31,49 .
In this work, we extend the use of Floquet-Bloch theory to accurately analyze spatiotemporal modulated graphene-based integrated structures, in order to provide a complete and accurate response and handle both guided and leaky modes precisely, for any modulation parameter.Unlike many previous studies concerned with spatiotemporal modulation of bulk permittivity and/or permeability, here we are engaged with spatiotemporal modulated surface conductivity of the graphene layer.As concrete examples, two typical plasmonic waveguides with both spatially and temporally modulated graphene conductivities are studied to investigate interband transition between different modes and intraband transition in the single mode waveguide.The dependence of mode conversion to the parameters of modulated wave such as direction and modulation depth is thoroughly studied.The accuracy of the proposed analytical method is verified through comparison with the numerical full wave results achieved by using COMSOL Multiphysics software.
We envisage that our study of graphene-based plasmonic waveguides paves the way towards a kind of transceiver front-end.The mere spatial modulation of graphene conductivity enables conversion of a graphene-based plasmonic waveguide to a leaky wave antenna.It follows that the time modulation of graphene conductivity not only changes the frequency of guided wave, but also the frequency of radiated wave.Moreover, the modulation parameters can be changed, to engineer the transmit and receive patterns separately in view of splitting the transmit and receive ports, and to enable the dynamic beam steering.In other words, the spatiotemporal modulation of graphene conductivity allows us to propose a multifunctional device.
The structure of the paper is as follows.First, the analytical model for analysis of spatiotemporal graphenebased structures is introduced and closed-form formulations are developed for two different graphene-based waveguides.Then, the accuracy of the developed formulas is verified through comparison with the results of previously reported works and also numerical results obtained using full wave simulation in COMSOL Multiphysics.
Using Floquet-Bloch theorem, we develop a semi-analytical method for fast and accurate modeling of spatiotemporal graphene-based waveguides.While full wave numerical simulation of the structure is rather complex due to the coupling of many harmonics of the fundamental spatial and temporal frequencies, the simple semianalytical method determines the electromagnetic fields with almost no computational cost.The semi-analytic method is applicable to various structures with spatiotemporal modulated reactance surfaces which are of prime interest in view of nonreciprocal devices with no external static magnetic fields.

Parallel plate graphene-based waveguides (PPGW)
The investigated graphene-based waveguiding structure is shown in Fig. 1.The structure includes two laterally infinite sheets of graphene, separated by a dielectric of relative permittivity ǫ r2 and height h s .The bottom graphene layer is fabricated on a substrate of relative permittivity ǫ r1 and a height much larger than the wave- length.To maintain the symmetry of the structure, the top graphene layer is fabricated on the same substrate.This waveguide supports two fundamental even and odd TM plasmonic modes 54 .Here, we assume that the top graphene layer is periodically modulated both in time and space, while the bottom layer has no modulation.The modulation can be done electronically by applying corresponding voltages to the top graphene layer [55][56][57] or optically by using optical pumping 58 .More conveniently, the modulation can be achieved by a guided wave supported by another quasi-TEM waveguide (with a frequency much smaller than the main frequency of PPGW).This modulation directly affects the graphene chemical potential which in turn modifies the graphene surface conductivity 35,59,60 .
The surface conductivity of graphene controlled by its chemical potential in the linear regime and in the absence of magnetic bias can be described by Kubo's formula including two terms describing interband and intraband carrier transitions 35,59,60 : where e is the electron charge, ℏ is the reduced Planck constant,k B is Boltzmann constant, µ c is the chemical potential, τ is the relaxation time of graphene, T is the temperature and ω is the angular frequency.
At THz frequencies where ωℏ < 2|µ c | , the interband term is negligible 54,61 .Furthermore, as all values of modulated chemical potential in this study are considered far from Dirac point, we assume that |µ c | ≫ k B T , which allows to employ Drude form of surface conductivity for the intraband term resulting in: As the frequency variation around the working frequency are small, it can be assumed that the frequency dispersion of Drude conductivity can be neglected.Thus, spatiotemporal modulation of chemical potential affects the graphene conductivity, directly.We can express this periodic modulation as a Fourier series: where ω m is the time modulation frequency, β m is the spatial modulation frequency, and σ q is the conductivity coefficient corresponding to each harmonic.Particularly, the fundamental harmonic (σ 0 ) is the average surface conductivity of the static graphene.
As Kubo's formula is valid under local equilibrium condition 62 , graphene conductivity (Eq.( 2)) can follow the spatial changes of chemical potential even if the spatial modulation period is as small as few nanometers or even less 13,63,64 .Besides, in high quality graphene with weak disorder and hence low collision frequency, if the modulation frequency of time varying chemical potential is much smaller than the main frequency and plasma resonance frequency, the graphene conductivity can follow the temporal modulation of chemical potential as well [65][66][67] .Therefore, we can effectively use time modulation frequencies up to THz band.
Having a spatiotemporal modulated graphene layer at the top, makes our structure be periodic both in space and time, thus according to the Floquet-Bloch theorem 68,69 , the electromagnetic response of the system would also be periodic.Therefore, the electromagnetic fields supported by the structure can be represented in terms of Floquet-Bloch modes.Indeed, through the propagation process, the fundamental Floquet mode may convert to other higher order modes and vice versa.Therefore, the amplitude of each harmonic does not remain constant.To account for this effect, and to consider the amplitude changes of each Floquet mode through propagation, we consider the propagation of all harmonics instead of only the fundamental one.The final response is then a combination of all harmonic responses.Therefore, the time domain plasmonic magnetic field, H y of TM waves supported by the waveguide can be expressed in the series form of time harmonics as: where ω n = ω 0 + nω m is the frequency of n th time Floquet harmonic, with ω 0 being the excitation frequency.H y (ω n , x, z) is the harmonic response of magnetic field at this frequency and can be expressed in three different domains of inside dielectric, above the top graphene layer and below the bottom graphene layer as: where A n (ω i ), B n (ω i ), C n (ω i ), D n (ω i ) are the amplitude coefficients of n th harmonic, and β(ω i ) is the propaga- tion constant at the frequency of ω i , while α 1n,i , α 2n,i are decaying factors of plasmonic waves in the above or www.nature.com/scientificreports/below dielectric with the relative permittivity of ǫ r1 and middle dielectric with the relative permittivity of ǫ r2 and are calculated as: is the wave vector of space Floquet harmonic and K 0i = ω i c .The upper-case sign is used when the propagating and modulation waves are in the same direction (codirectional modulation) while the lower-case sign is used when the propagating and modulation waves are in the opposite direction (contradirectional modulation).The coefficients a i in Eq. ( 5), are calculated using the boundary condition at x = 0 where only the fundamental harmonic is excited ( H(ω 0 , 0, h s 2 ) = 1 , H(ω n , 0, h s 2 ) =0 for n = 0).In this formulation, the conversion phenomenon which happens when different Floquet harmonics propagate inside the guiding structure, has been modeled through the series of space harmonics of Eq. ( 5) and also in the calculation of decay factors shown in Eqs.(6a), (6b).The electric field of the supported modes is then easily calculated using Eq. ( 5), and Maxwell's equations.Now that the formulation of the fields is achieved in the three regions, we need to apply boundary conditions at interfaces of z = ± h s 2 .To apply appropriate boundary conditions on the frequency domain tangential electric and magnetic field, the phasors E x and H y at each frequency of ω n are considered.Furthermore, since we have graphene layers at both interfaces, its impact on the discontinuity of tangential component of magnetic fields should be considered in the form of surface current as: where n is the unit vector normal to the graphene sheet and − → E t is the tangential electric phasor and − − → H t1,2 are magnetic field phasors of both sides of graphene sheet.It should also be noted that applying boundary conditions on both graphene sheets make the equations for different harmonics to be decoupled, so we drop i indices in the following equations for simplicity.Applying boundary conditions result in: Combing these equations results in a recursive relation for each harmonic as: For the convergence condition to be satisfied for these infinite set of linear homogenous equations of infinite unknowns |A n | s should be uniformly decreasing.Thus, the number of equations can be truncated to 2N + 1 (from −N to N ) and remaining equations can be represented in the matrix form of [K][A] = 0 , where [A] is the vector matrix of A n s, and [K] is the square matrix of coefficient with elements of: Finally, the dispersion equation can be achieved by applying the det(K) = 0 which is the condition to have nontrivial solution of the above equations.After finding β(ω i ) , using the dispersion equation, all A n s can be obtained in terms of the amplitude of the fundamental harmonic, A 0 .
For the case of sinusoidally modulated graphene, which is one of the most practical cases used in different applications 1,2,5,7,8,31 , the conductivity of the graphene can be written as: where M is the modulation depth and σ 0 is the average conductivity.In this case, the recursive formulation of Eq. ( 9) will result in the following equations: The unknown coefficients A n can be found in term of A 0 , by continued fractions as: By substituting A 1 and A −1 from Eq. ( 13) into Eq. (12a) for n = 0 , the dispersion equation is achieved as:

Single mode graphene-based waveguide (SMGW)
The second structure which is a single mode graphene plasmonic waveguide (SMGW) is shown in Fig. 2.This graphene-based waveguide is composed of a laterally infinite graphene layer above a grounded dielectric substrate.
Assuming that the conductivity of the graphene layer is modulated as Eq. ( 3), the analysis would be very similar to the previous section with Floquet-Bloch series of electric and magnetic fields, except for the boundary condition dictated by the metallic layer at the bottom of the structure.So, it yields the recursive relation as: consequently, the elements of coefficient matrix [K] will be (11) www.nature.com/scientificreports/ In the case of sinusoidally modulation, the coefficients of the three-term recursive of Eq. (12a) will be

Results and discussion
To investigate the accuracy of our Floquet-Bloch theory-based approach, here we compare our analytical results with full wave numerical simulation results obtained via COMSOL Multiphysics.PPGW and SMGW structures shown in Figs. 1 and 2 are studied.The dimensions and parameters characterizing the systems are listed in Table 1.The modulation of the graphene conductivity is assumed to be sinusoidal as in Eq. (11).Concerning PPGW, we reproduced Fig. 2a and b of Ref. 5 to examine our analytical approach.
A remark on the modulation depth ( M ) is in order.As a first approximation, graphene chemical potential can be written as 70 µ c = ℏv F πη V g + V offset , where ℏ is the reduced Planck constant, v F ∼ 10 6 m s is the Fermi velocity, V offset is the offset voltage caused by natural doping, V g is bias voltage and η can be estimated using a parallel plate capacitor model.For a typical chemical potential µ c0 = 0.5eV and a typical 100nm distance between the graphene layer and the ground plane, one finds V offset = 87.5V .Variation of bias voltages about ±18.4V results in chemical potentials between 0.45 and 0.55eV .It follows that a modulation depth M = 0.1 is not out of reach.
First, we consider PPGW-case1 with dimensions and parameters shown in Table 1.In previous works 5 , it has been shown that this structure can be used as an isolator: when the propagating and modulation waves are codirectional, most of the power will remain in the excited frequency of f 0 , while in the case of contra-direction, power of fundamental harmonic with the frequency of f 0 will be mostly coupled to the wave with frequency of f 0 + f m , and a non-remarkable power at the frequency of f 0 will reach to the end of the waveguide 5 .Here, we use our analytical model to investigate and study this phenomenon.The results of this analysis are shown in Fig. 3 and also Table 2.
Figure 3 illustrates the real part of longitudinal electric field E x , calculated using magnetic field of Eq. ( 5) and Maxwell's equations, of the waves propagating inside the modulated parallel plate waveguide for both codirectional (Fig. 3a-c) and contra-directional (Fig. 3d-f) cases.The non-reciprocal response of the spatiotemporal modulated PPGW (its directional response) is apparent from the results of Fig. 3. Excited Odd mode at f 0 is substantially coupled to the Even mode at f 0 + f m in contra-directional modulation (see Fig. 3e), while no significant mode conversion occurs in the codirectional case (see Fig. 3b,c).As wave propagates through spatiotemporal modulated structure and experiences frequency and propagation constant conversions, the power exchange between harmonics occurs periodically.Coherence length is a typical parameter to estimate the length of the device in which the most of the power is delivered to other harmonics from the fundamental one.Power exchange and coherence length can be calculated with the proposed analytical model.As shown in Fig. 3d, in the contra-directional modulation, the amplitude of the fundamental harmonic is at its lowest level at the coherence length of L c = 1.86 = 56µm and most of the power is transfered to the 1th harmonic.After this length, the power will return to the fundamental mode gradually and this process is repeated periodically.This determines an optimum length when designing isolators using these spatiotemporal structures.The results shown in Fig. 3 confirm the findings reported in Ref. 5 .
To have more clarification on the results shown in Fig. 3, more details are illustrated in Table 2.In this table, the middle results achieved in the analytical modelling are presented.One of the results presented in Table 2, is the effective refractive index of the spatiotemporal modulated waveguide, n eff , which is calculated using the dispersion relation presented in Eq. ( 14).This value has been calculated and illustrated both for contra-directional and codirectional cases.This table also illustrates the normalized amplitude coefficients, A n for ten different har- monics; n = −5 : 5 .The results of Table 2 exhibit a non-reciprocal response for the modulated PPWG waveguide as also reported in Ref. 5 .As shown in this table, in contra-directional modulation, the amplitude of the first harmonic, A 1 is much bigger than any other harmonics, while in the codirectional case, none of the harmonics has a significant amplitude, meaning that in this case no modulation happens and the power remains at the excitation frequency of f 0 as also illustrated in Fig. 3. ( 16) Table 1.The parameters used in the structures of Figs. 1, 2. By stronger modulation (larger modulation depth, M ), it is expected that more power be transferred at a shorter coherence length.To verify this prediction and also to test our analytical model for more examples, we study another parallel plate waveguide, this time with higher modulation depth.In this study, M has been increased to 0.2, and other parameters are kept the same as shown in Table 1 (PPGW-case1).The results of this study are shown in Fig. 4.This figure illustrates both 2D and 1D results.As shown in this figure, the coherence length would be equal to L c = 0.9 which is much shorter than the case with the lower modulation depth we studied earlier.This confirms the prediction that states by increasing the modulation depth, one can design isolators with shorter lengths.The results of Fig. 4, also indicates the wave amplification due to the power transfer from both fundamental mode of signal wave and modulating wave to the desired mode.Travelling wave parametric amplification is one of the first and well known properties of spatiotemporal waveguides already reported in previous works 71,72 .With the purpose of amplification, it is better to choose the waveguide length equal to the coherence length.
To verify the accuracy of our analytical model further, different dielectric materials are considered below and above graphene layers with dimensions and parameters shown in Table 1 (PPGW-case2).Here we compare our analytical results with the numerical full wave results achieved by performing simulation in COMSOL Multiphysics.This comparison is shown in Fig. 5, in which the magnetic field at the modulated graphene sheet (the upper layer) is illustrated.In the numerical simulation, coupled physics are defined for each frequency harmonic and modulated graphene is modeled using space variable surface current density in each physics.More detailed information of simulation setup is presented in section "Method".As shown in Fig. 5, a great agreement  www.nature.com/scientificreports/www.nature.com/scientificreports/ is observed between our analytical model and the numerical results achieved from COMSOL, confirming the accuracy of the developed model.Furthermore, our analytical model is about four times faster than the numerical simulation, using the same computational resources.Now we study SMGW structure shown in Fig. 2 with the geometrical and physical parameters listed in Table 1.The results of this study are shown in Figs. 6, 7, and also Table 3.In the single mode structure, spatiotemporal modulation can cause intraband transition which changes the frequency and propagation constant along the  www.nature.com/scientificreports/dispersion diagram of one mode.Figure 6 illustrates 2D results of the first three harmonics.As shown in this figure, coupling is still dependent to modulation direction and in the codirectional case, most of the power of the fundamental harmonic exchanges to other harmonics at the coherence length of L c = 1.6 = 160µm .However, this conversion to counterpart harmonics is almost equal and multiple harmonics are excited (as also illustrated in Table 3).It is owing to the fact that with this small modulation frequency, we remain close to the main frequency on the dispersion diagram and the wave experiences little dispersion, thus it is inevitable to excite other harmonics.In order to have more accurate results in Figs. 6 and 7, more harmonics have been contributed in the series of Eq. ( 5) for SMGW than PPGW. Figure 7 illustrates the magnetic field of three first harmonics on the graphene layer and compares it with numerical results achieved in COMSOL.Here, also we observe a very good agreement between the analytical and numerical results verifying the accuracy of the developed analytical model.By increasing the modulation depth, the coherence length can be shortened but simultaneously more power will be transferred to higher order harmonics.In contra-directional case, the results show no significant power coupling and the reduction of the amplitude of fundamental mode is mostly due to plasmonic wave attenuation.
In order to have an accurate comparison between our analytical results and the simulation results, modulation parameters are chosen so that no leaky mode is excited in lower order harmonics.However, our analysis is comprehensive and can be unconditionally used for leaky waves as well.
The PPGW structure (Fig. 1) could be fabricated in the following way 73,74 : The bottom substrate with relative permittivity ǫ r1 is deposited on a silicon wafer.Surface planarization or spin-coating of an extra thin layer improves the smoothness of the substrate surface.The chemical vapor deposition (CVD)-grown graphene on copper is then transferred to the smooth substrate.The middle layer with relative permittivity ǫ r2 can be deposited via atomic layer deposition (ALD) or plasma-enhanced chemical vapor deposition (PECVD) methods.The spincoating process can be utilized for polymeric materials.The second graphene layer is transferred to the middle layer.The top layer with relative permittivity ǫ r1 is deposited.A similar procedure can be employed to fabricate SMGW (Fig. 2): Via ALD or PECVD methods, a dielectric layer with relative permittivity ǫ r is deposited on a back metalized silicon wafer.The graphene is then transferred to the substrate.

Conclusion
Graphene-based spatiotemporal modulated structures are promising for different kinds of applications such as isolators, parametric amplifiers and leaky wave antennas.We have presented a Floquet-Bloch theory-based semi-analytical method for analysis of structures involving spatiotemporal modulated surface conductivities.While full wave numerical simulation of modulated structures is rather complex due to the coupling of many harmonics of the fundamental spatial and temporal frequencies, the semi-analytical method determines the electromagnetic fields with almost no computational cost.
We have studied the behavior of two different structures, namely parallel plate graphene-based waveguide (PPGW) and plasmonic single mode graphene-based waveguide (SMGW).The results are verified through comparison with full wave numerical simulations.We show extraordinary spatiotemporal effects such as nonreciprocity, mode conversion and amplification for these structures.We envisage that our study paves the way towards a kind of transceiver front-end.Indeed, the spatial modulation of graphene conductivity enables conversion of a graphene-based plasmonic waveguide to a leaky wave antenna.The time modulation of graphene conductivity allows to engineer the transmit and receive patterns separately and to enable the dynamic beam steering.In other words, the spatiotemporal modulation of graphene conductivity allows us to propose a multifunctional device.

Methods
The numerical results are obtained using full wave simulation by COMSOL Multiphysics.Particularly, three Frequency Domain Physics interfaces for three frequencies f 0 , f 0 + f m , f 0 − f m are coupled together.In the first Frequency Domain Physics interface, the left (right) boundary, associated with codirectional (contra-directional) case mentioned before, is assigned to be user defined port.This allows excitation of the structure with the appropriate plasmonic wave.Other external boundaries are defined as scattering boundary conditions.www.nature.com/scientificreports/Material properties at three frequencies f 0 , f 0 + f m , f 0 − f m are properly taken into account.The surface cur- rent density of the graphene layer depends on the space-time modulated conductivity of the graphene and the tangential component of the electric field.This naturally couples Frequency Domain Physics interfaces (for example, an electric field term of frequency f 0 − f m and a conductivity term of frequency f m result in a surface current density term of frequency f 0 ).

Figure 1 .
Figure 1.Graphene-based plasmonic waveguide structure.Parallel plate waveguide consists of two graphene layers of infinite width separated by a dielectric of relative permittivity ǫ r2 and height h s .

Figure 2 .
Figure 2. Single mode graphene-based plasmonic waveguide made of a graphene layer on top of metal backed substrate of relative permittivity of ǫ r and height h s .

Figure 3 .
Figure 3. Real part of normalized electric field E x of modulated PPGW-case1.(a-c) codirectional case and (d-f) contra-directional case of fundamental harmonic at f 0 and first harmonics at f 0 ± f m (the structure is excited at the frequency of f 0 ).

Figure 4 .
Figure 4. Real part of 1D normalized H y (a-c) and 2D normalized E x (d-f) for fundamental harmonic(a,d) and two first harmonics (b,c,e,f) of contra-directional modulated PPGW with stronger modulation depth of M = 0.2 .Dashed lines show the results obtained using our analytical model and solid lines illustrate COMSOL results.

Figure 5 .
Figure 5. Real part of normalized magnetic field, H y of PPGW-case2 at modulated graphene sheet boundary.(a-c) are for codirectional and (d-f) are for contra-directional modulation.(a,d) fundamental harmonic, (b,e) n = 1 , (c,f) n = −1 .Dashed lines show the results obtained using our analytical model and solid lines illustrate COMSOL results.

Figure 6 .
Figure 6.Real part of normalized electric field E x of modulated SMGW.(a-c) codirectional case and (d-f) contra-directional case of fundamental harmonic at f 0 and first harmonics at f 0 ± f m .The excitation frequency is f 0 .

Figure 7 .
Figure 7. Real part of normalized magnetic field H y of SMGW at modulated graphene sheet boundary. (a-c) are for codirectional and (d-f) are for contra-directional modulation.(a,d) fundamental harmonic, (b,e) n = 1 , (c,f) n = −1 .Dashed lines are obtained from our analytical model and solid lines are from COMSOL results.

Table 2 .
A n Coefficients of nth harmonics normalized to A 0 and n eff for structures with dimensions shown in Table1.

Table 3 .
A n Coefficients of nth harmonics normalized to A 0 and n eff for parameters of Table1.