Optimized nonlinear terahertz response of graphene in a parallel-plate waveguide

Third harmonic generation of terahertz radiation is expected to occur in monolayer graphene due to the nonlinear relationship between the crystal momentum and the current density. In this work, we calculate the terahertz nonlinear response of graphene inside a parallel-plate waveguide including pump depletion, self-phase, and cross-phase modulation. To overcome the phase mismatching between the pump field and third-harmonic field at high input fields due to self-phase and cross-phase modulation, we design a waveguide with two dielectric layers with different indices of refraction. We find that, by tuning the relative thicknesses of the two layers, we are able to improve phase matching, and thereby increase the power efficiency of the system by more than a factor of two at high powers. With this approach, we find that dispite the loss in this system, for an incident frequency of $2$ THz, we are able to achieve power efficiencies of $75 \%$ for graphene with low Fermi energies of $20$ meV and up to $35\%$ when the Fermi energy is $100$ meV.

FIG. 1. The metallic parallel-plate waveguide with graphene inside which forms the system being modelled. The inner material of the waveguide is polyolefin and the graphene is placed at the center of the waveguide at y = b/2. The pump field propagates in the +z-direction and is polarized in the x-direction.

I. INTRODUCTION
Graphene, as a zero-bandgap two-dimensional semiconductor with a linear electron band dispersion near the Dirac points has the potential to exhibit very interesting nonlinear optical properties [1][2][3][4] . The linear dispersion relation of the electrons near the Dirac points leads to a constant electron speed 5,6 . Thus, the intraband current induced in graphene by terahertz (THz) fields displays clipping as the amplitude of the incident field increases, which generates odd harmonics in the current and transmitted electric field [7][8][9][10] . Exploiting the nonlinear response of graphene enables one to produce higher frequency THz radiation through the generation of harmonics. Several experimental and theoretical groups have examined third-harmonic generation from graphene at terahertz frequencies. Almost all have employed a configuration where the field is normally incident on the graphene [11][12][13] .
However, here we consider a configuration where the radiation propagates in a metallic parallel-plate waveguide (PPW), with the graphene sheet lying at the midpoint between the two plates as shown in Fig. 1 14 . With this configuration, we increase the interaction time between the radiation and graphene, and thereby generate a larger-amplitude harmonic field.
We have shown in previous work that this configuration can increase the power efficiency of the system by more than a factor of 100, relative to the results for the normal-incidence configuration and that the power efficiency is relatively insensitive to the plate separation, but depends strongly on the Fermi energy 14 . However, in that work, we did not include the effects of pump depletion, self-phase modulation (SFM), and cross-phase modulation (XFM) in our calculations. In this work, we develop a coupled-mode theory including all propagating lossy modes to calculate the power efficiency for third-harmonic generation in a PPW and use this model to examine the impact of these effects on the power conversion efficiency.
For weak input fields, there is generally good phase matching in the waveguide between the pump field in the TE 1 mode at ω and third harmonic field in the TE 3 mode at 3ω.
However, as we shall show, when the pump field amplitude increases, the phase matching degrades due to SFM and XFM. To overcome this, we propose a new configuration in which the waveguide contains two layers of dielectric materials: cyclic polyolefin (n 1 = 1.53) and phenol formaldehyde resin (n 2 = 1.70) 17 . One goal in this work is to optimize the thickness of the dielectric layers and the Fermi energy of the graphene to obtain phase matching and thereby maximize the generated third-harmonic electric field.
The paper is organized as follows. In Sec. II we expand the electric field at the fundamental and third harmonic in terms of the lossy modes of the PPW and use the slowly-varying envelope approximation to derive the differential equations for the mode amplitudes. In Sec. III we compare the results obtained for the generated third-harmonic field, using our coupled-mode theory in the undepleted pump approximation, with pump depletion, and using a full calculation, which includes pump depletion and self-and cross-phase modulation.
In Sec. IV we propose new configuration PPW and demonstrate that this configuration allows us to essentially eliminate phase mismatch over a wide range of Fermi energies and input field amplitudes. Finally, in Sec. V we summarize our results.

II. THEORY
In this section, we first solve for the lossy linear modes of the waveguide with graphene present. We then expand the fields at ω and 3ω in terms of these linear modes, to derive the nonlinear coupled mode equations.
Our parallel-plate waveguide consists of two metallic plates placed at y = 0 and y = b, with the graphene midway between the plates at y = b/2 as shown in Fig. 1. The inner material of the waveguide is chosen to be cyclic polyolefin, with a refractive index of n 1 = 1.53, due to its compatibility with graphene, ease of fabrication and low loss at THz frequencies 17 . The THz wave propagates in the +z direction and is polarized in the x direction. For simplicity we take the plates to be perfect conductors that are infinite in the x direction.
From Maxwell's equations, we obtain the inhomogeneous wave equation, where E(r, t) is the total electric field and J (r, t) is the total current density in the graphene.
The incident electric field is taken to be harmonic with frequency of ω and a third-harmonic electric field is generated, so that the total electric field is and the total current density is Now, the current density can be broken up into its linear and nonlinear components as Here, J L is the current due to the linear conductivity of the graphene and is given by where σ (1) (ω ) is the linear conductivity of the graphene, with ω = {ω, 3ω}. The linear conductivity has both intraband and interband contributions 8,14,20 . At low Fermi energies (E F ≤ 20 meV), a photon with frequency of f 0 ≡ ω 0 /2π = 2 THz is not able to cause interband transitions. Thus, the dominant contribution to the linear conductivity is the intraband conductivity at ω 0 14 . However, for the third harmonic, the interband transition does contribute to some degree in the linear conductivity at low Fermi energies and so is included in our calculations 14 . The intraband contribution to the linear conductivity at zero temperature is proportional to the Fermi energy and is given by where E F is the Fermi energy and τ is phenomenological scattering time, which in this work is taken to be 50 f s. Thus, to limit linear loss, it is better to work at low Fermi energies.
The nonlinear current density of the graphene, J N L , arises from the third-order nonlinear conductivity, σ (3) . At 3ω it is given by The first term in Eq. (6) is the most important, as it is the source of the third-harmonic electric field. The next two terms are respectively related to SFM and XFM.
The nonlinear current density at ω is given by, The first term in Eq. (7) gives the nonlinear current at the graphene that results in to pump depletion and the next two terms represent SPM and XFM, respectively. The factors of 3 and 6 in Eqs. (6) and (7) arise from the number of ways of generating the nonlinear polarization in those cases.
In Eqs. (6) and (7), σ (3) is the nonlinear conductivity of the graphene. In this work, we use the theoretical expression of Cheng et al. 19 , which is derived from perturbative calculations at zero temperature, for electrons close to the Dirac points, with the neglect of the effects of the scattering (τ → ∞). Under these assumptions, where σ 0 = e 2 4 is the universal conductivity of the graphene, v F is the Fermi velocity of the electrons in the graphene, taken to be 1.1 × 10 6 m/s and in which, where Θ(x) is the Heaviside step function. The other nonlinear conductivities in Eqs. (6) and (7) are related to σ (3) (3ω; ω, ω, ω) by

A. Linear Modes
For harmonic waves travelling in the +z direction with angular frequency ω, the linear electric field for the n th transverse electric (TE) mode is given by where n = 1, 2, 3, ..., E n is the amplitude of the n th mode, andk n is the complex wave number for the field's y dependence. This wavenumber depends on the linear conductivity of the graphene and is obtained by enforcing the boundary conditions at the graphene, which leads to the following transcendental equation 14 whereφ n ≡k n (ω)b 2 and µ is the permeability of the dielectric. The complex propagation constant,β n , of the TE n mode is given bỹ where c 0 is the speed of light in vacuum and n 1 is the refractive index of the dielectric material. If there is no graphene, i.e., for a bare waveguide,k n ≡ k 0 n = nπ/b, where n is an integer.

B. Coupled Mode Equations
In this work, we take the input field at ω to be in the TE 1 mode at z = 0. We expand the field at ω in terms of the lossy TE n modes as and expand the generated third-harmonic electric field as where the summation is over all of the lossy modes propagating in the waveguide and the A n (z; ω) are slowly varying envelopes. Although this is not an exact expansion (as the lossy modes are not complete), we showed in our previous work 14 that using this expansion in the undepleted pump approximation, we obtain almost identical results for the generated third-harmonic field as were obtained using an exact Green function approach as long as the frequency is not close to the cut-off frequency.
The initial conditions are: A n (0; ω) = δ n,1 E input /sin(k n (ω) b 2 ) and A n (0; 3ω) = 0 where E input is the amplitude of the incident field at the graphene. We now employ our mode expansions to solve Eq. (1) for the fundamental and third harmonic fields including pump-depletion, SFM, and XFM. Using Eqs. (6), (7) and (14) in Eq. (1) along with the facts that ∇ · E = 0 and the modes are essentially orthogonal gives Now, using Eqs. (15) and (16) and employing the slowly-varying envelope approximation (i.e neglecting d 2 A n (z; 3ω) dz 2 ) we obtain the following differential equation for the amplitude of the electric field at 3ω for the m th mode (See appendix for details): Similarly, for the electric field modes at ω we obtain:

III. RESULTS
In this section, we solve the coupled dynamic equations for the amplitudes of the electric fields at 3ω and ω, given by Eqs. (18) and (19) respectively. In all that follows, we take the incident (pump) field to have frequency f 0 ≡ ω 0 /2π = 2 THz and take the plate separation to be 70 µm. We choose this plate separation because it is the largest separation for which there are only two propagating modes in the waveguide at 3ω. Note also that for this plate separation, only the TE 1 mode is a propagating mode at ω. In our previous work we showed that there is a perfect phase matching between the first mode at ω and third mode at 3ω when there is no graphene. Thus, we only include first and third modes in our calculations.
To solve the coupled equations of Eq. (18) and (19), we employ a Runge-Kutta algorithm; solving these coupled equations numerically takes less than one minute on an i7 processor.
In Fig. 2, we plot the generated third-harmonic electric field at the graphene as a function  18) and (19) and taking . Note however that we still include linear loss in all modes. Initially the third-harmonic grows rapidly, while the field at ω (see insets) decays exponentially until it is essentially gone after a propagation distance of about 2 mm. Oscillations in the third-harmonic arise due to the phase mismatch between the third-harmonic field in the n = 1 and n = 3 modes, which is why they persist even after the fundamental field is essentially gone.
Almost identical results are obtained for the generated third-harmonic field at low input fields when pump depletion is included, as seen in Fig. 2(a). However, as we increase the input field, we see a significant reduction in the generated third-harmonic field with pump depletion relative to the undepleted pump approximation. This is due to the increased transfer of power from the incident (pump) field to the third harmonic field. This can also be seen in the insets of Fig. 2(b) and (c), for strong input fields, where the pump field decays faster when the pump field is increased.
We now turn to the effects of self-and cross-phase modulation. As we increase the input field, we see in Fig. 2(b) and (c) that including SFM and XFM results in a decrease in the generated third-harmonic electric field. This is due to a degradation in the phase matching between the third mode at 3ω and first mode at ω due to SFM and XFM. For an input field of 15 kV/cm, this results in a 6% reduction in the peak field. As we shall see, the effect is much more significant at lower Fermi energies and/or higher input fields. Let us now examine the effects of SFM and XFM on phase matching in the PPW in more detail.
It is easy to show, using Eq. (14) that for a bare waveguide there is a perfect phase matching between the n = 3 mode at 3ω and the first mode at ω 14 . To generate a strong third-harmonic electric field we need to have a very small effective refractive index difference between these two modes in the presence of graphene loss and SFM and XFM. This effective refractive index difference for a lossy waveguide with SFM and XFM is approximately given by where and The first terms in each of Eqs. (21) and (22) are the input-field-independent effective refractive indices for the first mode at ω and third mode at 3ω, respectively. The second terms are added to approximately account for the change in the effective index due to the SFM and XFM. Using Eqs. (18) and (19) we obtain In deriving these expressions, we have only included the terms proportional to the square of the electric field at ω and have taken the pump field to be given by its value at z = 0.
We find that the terms proportional to the square of the third-harmonic electric field are negligible relative to the linear electric field. However, in our full numerical calculations, all terms are retained, as is the z-dependence of A 1 (ω).
The effective index difference between the third mode at 3ω and the first mode at ω as a function of Fermi energy is shown in Fig. 3 for different input pump field amplitudes.
When E input = 0, ∆n ef f linearly increases with Fermi energy, due to the dependence of the propagation constants on the doping level of the graphene. As the input field increases, We now consider the power efficiency S ef f of the device, which is defined as the ratio of the power in the third harmonic at the end of the waveguide to the power in the fundamental at the beginning of the waveguide. In Fig. 4 we plot the maximum power efficiency as a function of Fermi energy in three different schemes: undepleted-pump approximation, with pump depletion, and full calculation. In all cases, the length of the waveguide is chosen to be the distance at which the power in the third harmonic field is a maximum, as seen in In this section we examine the effect of the layers thickness on the linear phase mismatch.
We begin by examining the linear modes. In Fig. 6 we plot the normalized TE 1 mode at ω and TE 3 mode at 3ω for a 2 THz field in a 70 µm waveguide in our original configuration and in our new configuration. As an example we choose d 1 = 25 µm. Note that the TE 1 mode peaks at the centre of the waveguide, inside of the low-index (n 1 ) material, while the TE 3 mode also has peaks inside of the high-index (n 2 ) material. As a result, it is expected that by increasing the width of the high-index material, we can raise the effective index of the TE 1 mode more than that of the TE 3 mode and thereby modify the index mismatch.
This new configuration can be used to not only help to decrease the phase mismatching induced by the SFM and XFM but to overcome the linear phase mismatch introduced by the graphene. Note also that the new configuration leads to an increase in the amplitude of the field at the graphene in the TE 1 mode. This will yield a slight increase in the generated field for a given input power.
2 (3ω) to zero. It is seen that perfect phase matching (∆n ef f ≡ 0) occurs at two points for each value of E F as we increase d 1 . More importantly, we see that we can reduce ∆n ef f by up to 0.1, which is more than enough to compensate for the effective index difference shown in Fig. 3

B. Power Efficiency
In Fig. 8 we plot the power efficiency of the waveguide in our new configuration as a function of d 1 for an incident field of 5 kV/cm at Fermi energies of E F = 20 meV, E F = 50 meV, and E F = 100 meV. Decreasing the Fermi energy leads to higher power efficiency due to the increase in the nonlinear conductivity 20 and reduced loss due to a decrease in the linear conductivity 8,19 . Note that when d 1 = 35 µm, our new configuration PPW is identical to our original PPW, and so we obtain the same efficiency as is given in Fig. 4(a). At this low input power, we know from Fig. 3 that the effects of XFM and SFM are almost negligible except for E F = 20 meV. Therefore, we expect to obtain the peaks in the efficiency close to the values of d 1 where ∆n ef f is zero in Fig. 7. This is indeed what we find, but with small shifts that arise from the need to also compensate for XFM and SFM. As expected from Fig. 3, this shift is largest for the Fermi energy of 20 meV. We also note that for this small input field, the increase in the efficiency over the initial-configuration PPW is rather modest (< 30%). In Fig. 9 we plot the power efficiency of the waveguide in the new configuration as a function of d 1 for a Fermi energy of 20 meV for input fields of 5, 10, and 15 kV/cm. As in Fig. 8, the efficiency peaks at two different values for each input field. However, for the higher input fields, we see that these peaks are much larger and have shifted to values of d 1 that are closer to the minimum in ∆n ef f seen in Fig. 7. Both of these effects are the result of the need to compensate for the much larger phase mismatch that arises from XFM and SFM at high input fields.
In Fig. 10, we compare the power efficiency of the waveguide in the original configuration with and without SFM and XFM to the efficiency of the waveguide in the new configuration with the full calculation. In the new configuration, the optimized power efficiency is improved such that it equals or improves upon the results we obtained for the original configuration when SFM and XFM is neglected. This is because our new configuration not only overcomes the phase mismatch induced by SFM and XFM but also that induced by the linear response of the graphene. For example, for E F = 20 meV, and E input = 10 kV/cm the power efficiency increases from 28% to 48% and for E input = 15 kV/cm it increases from 35% to 75%.
Due to the experimental difficulties in achieving a uniform doping over graphene sheets that are millimetres in length, achieving low Fermi energies is very challenging. Therefore, we now examine what efficiencies can be achieved for higher Fermi energies if we move to higher input fields. In Fig. 11(a), we plot the optimized power efficiency of the waveguide in the new configuration as a function of input field for different Fermi energies. Note that the optimized length and d 1 are different for each input field and Fermi energy (See Fig. 11(b) and (c)). We note that for all three Fermi energies, the efficiency initially rises with input power, but it then reaches a peak and settles to a value of around 30% at high input fields.
To understand this, consider Fig. 3, where we see that for E F = 20 meV, the index mismatch due to XFM and SFM is already 0.06 at an input field of 15 kV/cm. It is easy to see therefore that for the fields considered in Fig. 10, we will quickly reach index differences greater than 0.1 which is the maximum that can be compensated for using our new configuration PPW. Therefore, the input field at which the efficiency peaks for a given Fermi energy is the field at which the nonlinear effective index difference reaches about 0.1. Due to the strong dependence of XFM and SFM on the Fermi energy, this field amplitude is different for the different Fermi energies. The highest efficiencies are obtained for the lowest Fermi energy of 20 meV, largely because the loss is lower in this system and because the higher nonlinearity means that the structure length is less. Note that for all three Fermi energies, the peak efficiencies occur for devices with a length of only a few hundred microns. The reason that the high-field efficiency is essentially independent of Fermi energy is because in all cases, the pump is depleted over a distance that is less than the linear loss distance (1.2 mm, 0.48 mm and 0.24 mm for Fermi energies of 20 meV, 50 meV, and 100 meV, respectively) and so the different losses at the different Fermi energies do not play a significant role. Therefore, we find that very good efficiencies can be obtained for higher Fermi energies if one can attain the higher input fields. For example, for a Fermi energy of 100 meV, we are able to obtain a power efficiency of 30% at an input field of 50 kV/cm. Because at high-fields, our structure is not able to compensate for SFM and XFM, we find (not shown) that efficiencies of up to 40% can be obtained at high fields and large Fermi energies even in our original configuration. This is a very promising configuration that we believe should be achievable in the lab using current graphene samples and THz sources 11,21 .

V. SUMMARY
We have developed a coupled-mode theory for the propagating lossy modes of the pump and third harmonic fields in a PPW to calculate third harmonic generation, including pump depletion, SFM, and XFM. We find that SFM and XFM degrades the phase matching between the TE 1 mode at ω and TE 3 mode at 3ω and thereby decreases the generated third-harmonic electric field. We have shown that one can overcome the phase mismatch due to SFM and XFM by designing a new configuration PPW. We found that by optimizing the dielectric layer thickness, the power efficiency can be increased by more than a factor of two relative to the original configuration. We have also shown that even for graphene with Fermi energy of 100 meV, where the nonlinearity is relatively modest, efficiencies of up to 30% can be achieved for input field amplitudes of 50 kV/cm. We therefore believe that our PPW system is an excellent platform to produce and examine harmonic generation in graphene.
If we define S (1) n y = sin(k n (ω)y),  To determine C, we take integrate S (1) n y from y = b/2 − to y = b/2 + , for b.
Then, we obtain The electric field below and above the graphene is defined as below (y, z; 3ω) =E n e iβn(3ω)z sin(k n (3ω)y) x 0 < y < b/2 (A.6) above (y, z; 3ω) =E n e iβn(3ω)z sin(k n (3ω)(y − b)) x b/2 < y < b using Maxwell's equations we have, ∇ × E below (y, z; 3ω) =β n (3ω) 3ωµ 0 E n e iβn(3ω)z sin(k n y) y (A.8) The current density at the graphene is given by where the current density at the graphene for n th mode is related to the field by We also have that, using the slowly-varying envelop approximation, where we neglect the second derivative of the envelope, that −β 2 n (3ω)A n (z; 3ω)}e iβn(3ω)z S (1) n (y). From this we obtain Eq. (18) for the differential equation of the amplitude of the electric field at 3ω for mth mode. A similar calculation yields Eq. (19) for the differential equation for the electric field at ω. Note that we have confirmed numerically that the slowly-varying envelope approximation, where we neglect the second derivatives of the envelopes is an excellent approximation for all fields and Fermi energies.