Terahertz rectangular waveguides with inserted graphene films biased by light and their quasi-linear electromagnetic modeling

Novel rectangular waveguides with graphene inserts biased by light are proposed herein. The graphene films short the conductor plates of waveguides and support the localized transverse-electric modes. Their electric fields are parallel to the wide walls of these waveguides, and the eigenmodes have decreased conductor loss. The designs do not involve the conductor and graphene strips with their sharp edges, and the loss associated with the current crowding effect is excluded. The waveguides are treated in the quasi-linear regime using a rigorous field matching method, and the complex dispersion eigenmodal equation is solved using a validated iteration algorithm. At the terahertz frequencies of amplification, where the real part of graphene conductivity is negative, a gain increase is found with the eigenmodal number. This gain can be tuned by the waveguide geometry, dielectric filling, and the level of quasi-Fermi energy. The ideal waveguide theory is corrected using a perturbation approach and the Drude model of surface resistance of waveguide plates.


Introduction
Today, many efforts are being devoted to developing passive and active graphene devices [1][2][3]. These include interconnects, antennas, attenuators, transistors [4,5], and nonlinear devices [6]. Many techniques are used to realize these and other prospective components, such as chemical doping [7], integration of graphene with other materials [3], and graphene nanoribbon or nanodot elements [8]. Some problems of active and passive components are solved by staking the graphene nanosheets or using the quantum-mechanical interaction of graphene layers [9]. An analysis of the current state of graphene and other one-atom materials is given in [3], where the industrial mass production of 2-D materials and components is predicted for the next 30 years of this century.
The development of new components, with their analytical and semi-analytical models, will necessarily make these efforts easier. For instance, to avoid the lengthy ab initio computation of graphene components, the semiclassical scalar and dyadic models of conductivity of graphene are known based on the approach described by Kubo [10][11][12][13]. A theory based on a non-Kubo technique is discussed in Ref. [14]. In recent years, much attention has been paid to the experimental verification of conductivity models [15][16][17] and theories of guided plasmons in a multi-strip environment [18]. For instance, in Ref. [15], a difference between the results provided by a Kubo formula and the measurements is explained with the random microstructural properties of graphene films deposited on a dielectric substrate. A statistical analysis of defects or irregularities of graphene and dielectric layers is performed in Ref. [19]. To reduce this influence, researchers have recommended a high-level polishing of silicon dioxide or the use of complementaryto-graphene boron nitride [20] to avoid the morphological obstacles of graphene covering these materials. The influence of impurities on the atomistic level is considered in Ref. [16], where a modified Drude conductivity formula is given for the few-terahertz range.
Special interest is focused on models of the nonlinear properties of graphene when its chemical potential depends on the applied AC voltage [6]. Here, mixing of signals or their parametric amplification is possible [21]. An analytical model of dynamic graphene conductivity was developed in Ref. [10], where it was shown that this material biased by a pumping light demonstrated negative conductivity, and terahertz signal amplification was possible. This effect is discovered and described numerically and in experiments in Refs. [12,[22][23][24][25]. Most papers from this field are on the planar designs when graphene films and conductor strips are placed on a surface of a dielectric substrate, and diffraction of the EM waves is studied in numerical simulation or by measurements. An increase in amplification has been demonstrated with the use of multiple layers of graphene, and formulas for calculating the effective conductivity are given in Ref. [26]. Traveling wave amplification is studied numerically in a conductor-backed slot-waveguide in Ref. [27]. A detailed review of the research on optically pumped terahertz amplifiers and lasers is provided in Ref. [9].
In this paper, we propose a novel design for graphene loaded waveguides illuminated by light. In these waveguides, there are no sharp-edged conductors and graphene strips causing current crowding effect. The wall loss is decreased using the transverse-electric modes with the electric field parallel to the conductor plates. The modal fields are tightly pressed to the graphene film due to its surface nature, and the evanescent effect is realized by properly chosen crosssection geometry. In addition, the multiple-element effect is excluded, which usually leads to a rather dirty spectrum response. A new result found in this paper is the increased gain with the modal number as a result of the designed waveguide structures. In Sect. 2, the waveguides and their semianalytical theory are given using the quasi-linear approach. The simulation results for the waveguides are presented in Sect. 3. The obtained results are discussed and conclusions drawn in Sect. 3.

Proposed waveguides and their electromagnetic quasi-linear theory
In contrast to planar components proposed previously by other authors, our waveguides have a rectangular enclosure design made of an entirely ideally conducting electric wall (x = ±a, y = 0, b) = 0 (Fig. 1a) or combined electric (y = 0, b) and magnetic (x = ±a) = 0 walls with the cross-section size comparable to several hundred or tens of micrometers depending on the dielectric filling.
The graphene film shorts the wide-walls, and the lowloss TE y modes can be employed similarly to the known parallel-plate and the non-radiative dielectric (NRD) waveguides used in the millimeter-wave and terahertz frequency range [28,29]. The biasing or pumping infrared light can be launched along the waveguide and illuminates uniformly the graphene film due to reflections from the sidewalls (Fig. 1a). In the case of a magnetic wall (Fig. 1b), the light is directly incident to the graphene film from outside.
The complexity and accuracy of the electromagnetic (EM) theory of these waveguides are highly dependent on the model of graphene conductivity ↔ g used here, known as a Kubo formula. For instance, if anisotropy of conductivity should be taken into account at increased terahertz frequencies, then the hybrid modes having all six components of the EM field appear in consideration [11]. It is known that at strong signal fields, the conductivity of graphene shows moderate nonlinearity, and again, this should be treated. Some additional factors involve the substrate influence, impurity of graphene and its doping, quantum discretization of energy for nano-size strips, and edge effect.
Our theory supposes the perfect electric and magnetic walls initially, scalar linear conductivity of graphene σ g , and weak signals only. It allows a separate description of the TE y and TM y modes in a quasi-linear manner in the considered waveguides. The subject of our study is the TE y modes, because other modes TM y are damped due to the shortening effect of graphene film, although they can propagate if the real part of the conductivity of graphene is negative.
To describe these TE y modes, the theory presented by Balanis [30] is pertinent, allowing us to write all required field components using one vector potential function F m (1,2) in each subdomain 1, 2 divided from each other by an infinitely thin conducting graphene film (Fig. 1). For the waveguide of the perfect electric wall (Fig. 1a), they are where 0 is the vector unit, A m (1,2) are the unknown modal a m p l i t u d e s , , f is the driving frequency, ɛ r and μ r are the relative permittivity and permeability, correspondingly, c is the velocity of light, j = √ −1 , and k z (m) is the unknown modal longitudinal propagation constant which should be defined from our treatment.
In the case of electric/magnetic wall enclosure (Fig. 1b), these vector potential functions are where B m (1,2) are the unknown modal amplitudes. The fields are obtained from Ref. [30], p. 275, and they are in the case of Fig. 1a where ɛ = ɛ 0 ɛ r and μ = μ 0 μ r are the absolute permittivity and permeability of dielectrics inside waveguides, with ɛ 0 and μ 0 as the corresponding vacuum constants.
The electric/magnetic wall geometry fields are described in a similar manner using the functions Φ m (1,2) in Ref.
(3) instead of F m (1,2) . The chosen functions F m (1,2) and Φ m (1,2) satisfy the boundary conditions for the EM fields everywhere except the graphene film plane x = 0. At this boundary, the following equations are written, valid for both waveguide geometries (Fig. 1a, b) with graphene conductivity according to Ref. [10,25] (1) (1) where e is the electron charge, ℏ is the reduced Planck's constant, k B is the Boltzmann constant, T is the temperature in Kelvin, τ is the graphene relaxation time, E F is the quasi-Fermi energy, and G E, (1,2) and Φ m (1,2) into Eqs. (3) and (4), two sets of linear algebraic homogeneous equations are obtained regarding the unknown amplitudes A m (1,2) and B m (1,2) for the electric and electric/magnetic wall enclosures, correspondingly. Calculating the determinants of these equations at their zero points gives the equations to obtain complex k z (m) . The formulas convenient for an iteration method of the solution are given below.
The propagation constant k z (m) for the electric wall waveguide (Fig. 1a) is found from Eq. (6) For the electric/magnetic wall waveguides (Fig. 1b), a similar Eq. (7) is Additionally, Eqs. (6) and (7) can be solved using any complex root search algorithm. Still, it has an approximate it can be used as an initial value. On the other hand, a = ± ∞ leads to an H-waveguide parallel-plate configuration with an analytical formula (8) for propagation constant

Results of simulation and discussion
Initially, let us analyze the waveguides without graphene inserts. Analogously to the non-radiative dielectric waveguides [28], their vertical size b is chosen to provide the condition of non-propagation for the TE (m) y eigenmodes in the hollow rectangular waveguide. It allows us to study the waves tightly coupled to graphene and reduce the influence of electric or magnetic walls at x = ± a.
Considering the chosen modal propagation constant k z (m) along the z-axis in expressions (1), (2) and satisfying the energy conservation law, then the obtained formulas (6) and (7) describe the backward modes with the phase velocity oriented against the group one k z (m) = −(β (m) + jα (m) ). It follows that if α (m) > 0, then a mode propagates exponentially decreasing along the z-axis according to formulas (1) and (2). The negative sign of α (m) is with the amplification occurring in graphene-based waveguides when Re g < 0.
Before applying formulas (6) and (7) to study these waveguides, the proposed iteration method to calculate the propagation constant should be validated. For this purpose, the parameters of graphene were chosen according to Ref. [27], where these values are properly discussed. It allows some qualitative comparisons with the mentioned paper and other works of the authors from this field.
It has been found that for the geometries a/b > 0.02 from Fig. 1a, b, the propagation constant k z (1) is close to the value provided by the parallel-plate waveguide model (8), and no additional iteration steps are required except one.
At very narrow waveguides, when a/b < 0.02, several iterations are needed to obtain accurate values for the propagation constant. Figures 2 and 3 show the convergence curves for one of the worst-case scenarios when (a/b) = 0.01 (waveguide geometry Fig. 1a). A notable feature of this algorithm is its numerical stability, as seen from Figs. 2 and 3.
Consider the waveguides (Fig. 1a) filled by quartz or boron nitride (ɛ r = 3.84) and calculate the frequency dependence of the loss α (m) and phase β (m) constants of the first two modes at two different values of quasi-Fermi energy E F . As an initial value, the propagation constant (8) of an infinite parallel-plate waveguide is chosen. The waveguide sidewalls are far enough from the graphene layers (a/b = 1), and only one iteration is needed to reach a converged value (Figs. 4, 5).
It is seen that after a specific frequency where Re g = 0 , the loss constant α (m) < 0 (Fig. 4), and its absolute value increases with E F . This effect and frequency is in accord with the results of [27] where the scattering of terahertz waves on multi-strip structures was studied.
At the difference to |α (m) |, the phase constant β (m) /k 0 decreases with the value of the quasi-Fermi energy E F (Fig. 5).
The increased gain with the modal number is explained by the increase in the sub-wavelength localization of the field near the graphene film (Fig. 6). It is caused by a Due to the abovementioned effect, the loss α (m) and phase β (m) constants show strong dependence at small values of the waveguide height b (Figs. 7 and 8), and it can be used to obtain a gain of a desired level without tuning E F . On the other hand, a geometry can be chosen, providing stable modal parameters towards manufacturing intolerances at larger b/a.
The abovementioned analysis is confirmed by a study of the field distribution near the graphene film (Fig. 9).  At smaller b waveguides, the field is more closely concentrated near graphene, leading to significant amplification. Again, this localization is caused by the substantial reflection of waves from the evanescent waveguide parts in ± x-directions. Thus, varying the height b of the waveguide allows for field localization control in the orthogonal direction.
In addition to the influence of the waveguide height b, the dependence of the modal parameters on the waveguide width a is simulated. Figure 10 shows an example of these calculations. It is seen that the influence of the sidewalls is strong enough when a/b ≪ 1. It is interesting to note that magnetic walls increase |α (1) /k 0 |, in contrast to the electric ones. For most practical applications, the sidewalls do not affect the propagation constant of the TE (1) y mode, and this waveguide would be rather stable towards manufacturing inaccuracies.
In addition to the simulation of these ideal waveguides, the dielectric and conductor loss must be considered. The dielectric loss can be calculated using formulas (6), (7), or (8) and complex representation of dielectric permittivity. The use in this paper of quartz (tan δ = 0.007) as waveguide filling shows a negligible loss in comparison to the graphene effect.
The conductor loss constant α c (m) is calculated here with the known perturbation approach where the fields and propagation constant are taken from the ideal waveguide model (3) where R s is the conductor surface resistance, and integrations are along the waveguide perimeter L and cross-section S, correspondingly.  normalized |E x | field component dependence in the x-direction in the right part of the rectangular waveguide (Fig. 1a). Waveguide parameters: a = 0.009 mm, ɛ r = 3.84, μ r = 1, and f = 4 THz. Graphene parameters used for calculations of conductivity (5) are from Ref. [27]: E F = 40 meV and τ = 1 ps. Number of iterations N = 1 As the surface resistance R s , for comparison, two approximations are used. The first is just the conventional skineffect formulas used to calculate, for instance, micro-strips [31,32] where Z 0 = (1 + j) √ 0 2 , k = (1 − j) √ 0 2 , σ is the specific conductivity at DC, and t is the thickness of the waveguide walls.
At terahertz frequencies, one can consider the displacement current in conductors [34], and the surface resistance of a thick sheet of metal is where the Drude conductivity σ D is where τ c is the mean free path time of the electrons in the conductors [33]. Figure 11 shows a comparison of calculations α (1) /k 0 obtained for an ideal conductor rectangular waveguide and those calculated with the help of formulas (10) and (11). It (10) R s = Re −jZ 0 cot (kt) , is seen that loss essentially corrects the calculation results of loss/gain constant, as expected. Nevertheless, the proposed waveguide shows still amplification, but the frequencies where the amplification appears are shifted to higher values. Another conclusion from Fig. 11 is that the Drude and conventional skin-effect formulas show comparable results. Still, the use of the first formula in the terahertz range is preferable due to higher accuracy according to the analysis of Refs. [34,35].
In full confirmation of the previous theoretical and experimental works cited above, our simulation shows that our proposed optically biased graphene-based waveguides can provide an essential gain for modes propagating in them.
The developed lossy model allowed us to plot more realistic E x (x, z) field and current J z (y, z) = -σ g E z (y, z) distributions Fig. 11 Loss constants of TE (1) y mode of studied rectangular waveguide (Fig. 1a) versus frequency. Waveguide data: a = b = 0.009 mm, ɛ r = 3.84, μ r = 1, t = 20 μm, = 5.96 ⋅ 10 4 Sim∕mm , and c = 3.6 ⋅ 10 14 s are for copper [33]. Graphene parameters used for calculations of conductivity (5) are from Ref. [27]: E F = 40 meV and τ = 1 ps. Number of iterations N = 1  (Figs. 12, 13). The first one is normalized by the value of |E x (x = a, z = a)|, and the second by |J z (z = h)|. The figures take into account the corrected values of k x and k z . In Fig. 12, the exponential dependence of the field along the x-and z-directions can be observed.
The longitudinal current (Fig. 13) takes into account the Drude model corrections as well, but the y-dependence is taken from the ideal model according to the perturbation technique, and then J z (y = 0, h) = 0 . Similarly to the electric field, the current increases exponentially along the z-axis.

Conclusions
Novel terahertz graphene-based waveguides and their quasi-linear analytical models have been proposed for lightpumped applications. In these components, at the difference to planar integrated circuits, the following design ideas have been offered: 1. The rectangular waveguides have been proposed evanescent for propagation of the TE (m) y and TM (m) y (m > 0) modes for housing the graphene films illuminated by light inducing negative graphene conductivity in the terahertz range 2. These graphene films short two parallel plates of hollow rectangular waveguides while other walls are made of conductors or left open 3. The designs do not employ conductor and graphene strips or micro-strips having sharp edges that cause a current crowding effect It enables the following to be achieved: a. To decrease the conductor-plate-associated loss in the case of TE (m) y modes, similarly to the non-radiative dielectric rectangular waveguides b. To exclude the loss caused by sharp edges of conductor and graphene strips, unlike known planar designs c. To provide a better concentration of the field close to the graphene surfaces to gain increase d. To provide enhanced conditions for the integration of many of these waveguides of the increased terahertz power e. To be open using the multilayer graphene films or heterostructures for more effective employment of light power f. To be designed for the current injection lasing configurations in the case of conductor/open wall housing configuration These proposed waveguides have been modeled using the EM field matching method valid for low-level signal applications when the terahertz field is weak enough so as not to influence the graphene conductivity. Otherwise, the solution of nonlinear Maxwell equations should be employed [21].
The formulas for the iterative calculation of complex propagation constants have been obtained and validated. In these waveguides, similarly to the known planar designs, the negativity of the real part of graphene conductivity caused by infra-red pumping light provides the exponential increase in signals along the longitudinal propagation direction, as shown in the simulation. It revealed a gain increase with the mode number caused by the increasing concentration of high-order modal field to the graphene surface in the considered frequency band. The parameters of modes can be regulated by the level of quasi-Fermi energy E F , dielectric filling ɛ r , waveguide height b, and width 2a. The developed model considers the conductor loss according to the Drude theory of metals and dielectric loss at terahertz frequencies.
Funding Open Access funding provided by NTNU Norwegian University of Science and Technology (incl St. Olavs Hospital -Trondheim University Hospital).
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creat iveco mmons .org/licen ses/by/4.0/.