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

###### Abstract

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 THz, we are able to achieve power efficiencies of for graphene with low Fermi energies of meV and up to when the Fermi energy is meV.

## 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 mikailov (); Mikhailov2010 (); Gu (); Glazov (). The linear dispersion relation of the electrons near the Dirac points leads to a constant electron speed castro (); sarma (). 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 Bowlan2014 (); Alnaib2015 (); Ibra2015 (); ibraheem2014 (). 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 Hafez (); Paul (); PBowlan (). 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 parvin ().

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 , 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 parvin (). 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 mode at and third harmonic field in the TE mode at . 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 () and phenol formaldehyde resin () Nestor (). 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 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 as shown in Fig. 1. The inner material
of the waveguide is chosen to be cyclic polyolefin, with a refractive index of , due to its compatibility with graphene, ease of fabrication and low loss at THz frequencies Nestor (). The THz wave propagates in the direction and is polarized in the direction. For simplicity we take the plates to be perfect conductors that are infinite in the direction.

From Maxwell’s equations, we obtain the inhomogeneous wave equation,

(1) |

where is the total electric field and 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

(2) | ||||

and the total current density is

(3) | ||||

Now, the current density can be broken up into its linear and nonlinear components as . Here, is the current due to the linear conductivity of the graphene and is given by

(4) |

where is the linear conductivity of the graphene, with . The linear conductivity has both intraband and interband contributions Alnaib2015 (); parvin (); horng (). At low Fermi energies ( meV), a photon with frequency of THz is not able to cause interband transitions. Thus, the dominant contribution to the linear conductivity is the intraband conductivity at parvin (). 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 parvin (). The intraband contribution to the linear conductivity at zero temperature is proportional to the Fermi energy and is given by

(5) |

where is the Fermi energy and is phenomenological scattering time, which in this work is taken to be . Thus, to limit linear loss, it is better to work at low Fermi energies.

The nonlinear current density of the graphene, , arises from the third-order nonlinear conductivity, . At it is given by

(6) | ||||

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,

(7) | ||||

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 and in Eqs. (6) and (7) arise from the number of ways of generating the nonlinear polarization in those cases.

In Eqs. (6) and (7), is the nonlinear conductivity of the graphene. In this work, we use the theoretical expression of Cheng et al. cheng (), 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,

(8) |

where is the universal conductivity of the graphene, is the Fermi velocity of the electrons in the graphene, taken to be m/s and

(9) |

in which,

(10) |

where is the Heaviside step function. The other nonlinear conductivities in Eqs. (6) and (7) are related to by

(11) | ||||

### ii.1 Linear Modes

For harmonic waves travelling in the direction with angular frequency , the linear electric field for the transverse electric (TE) mode is given by

(12) |

where , is the amplitude of the n mode, and is the complex wave number for the field’s 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 parvin ()

(13) |

where and is the permeability of the dielectric. The complex propagation constant, , of the TE mode is given by

(14) |

where is the speed of light in vacuum and is the refractive index of the dielectric material. If there is no graphene, i.e., for a bare waveguide, , where is
an integer.

### ii.2 Coupled Mode Equations

In this work, we take the input field at to be in the TE mode at . We expand the field at in terms of the lossy TE modes as

(15) |

and expand the generated third-harmonic electric field as

(16) |

where the summation is over all of the lossy modes propagating in the waveguide and the are slowly varying envelopes. Although this is not an exact expansion (as the lossy modes are not complete), we showed in our previous work parvin () 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: and where 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 and the modes are essentially orthogonal gives

(17) | ||||

Now, using Eqs. (15) and (16) and employing the slowly-varying envelope approximation (i.e neglecting ) we obtain the following differential equation for the amplitude of the electric field at for the mode (See appendix for details):

(18) | ||||

where, and . Similarly, for the electric field modes at we obtain:

(19) | ||||

## Iii results

In this section, we solve the coupled dynamic equations for the amplitudes of the electric fields at and , given by Eqs. (18) and (19) respectively. In all that follows, we take the incident (pump) field to have frequency THz and take the plate separation to be . We choose this plate separation because it is the largest separation for which there are only two propagating modes in the waveguide at . Note also that for this plate separation, only the TE 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 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 of for a Fermi energy of meV for three different input field amplitudes ( , and kV/cm). We present the results of three different calculations: undepleted pump approximation (red solid curves), calculation with pump depletion (green dotted curves), and the full calculation that includes the SFM and XFM and pump depletion (blue dashed curves). To demonstrate the importance of pump depletion, we
first present the results of the calculation when the SFM, XFM and pump depletion are neglected. This is accomplished by keeping only the first term on the right hand side of Eqs. (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 mm. Oscillations in the third-harmonic arise due to the phase mismatch between the third-harmonic field in the and 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 and first mode at due to SFM and XFM. For an input field of kV/cm, this results in a 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 mode at and the first mode at parvin (). 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

(20) |

where

(21) |

and

(22) |

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 , 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

(23) | ||||

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 . 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 -dependence of .

The effective index difference between the third mode at and the first mode at as a function of Fermi energy is shown in Fig. 3 for different input pump field amplitudes. When , 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, increases for all Fermi energies, but increases the most for low Fermi energies. At the lowest Fermi energy of meV and input field of kV/cm, is very large, reaching a value of approximately . This strong dependence on Fermi energy has its origin in the strong dependence of the nonlinear conductivity on Fermi energy, as seen in Eqs. (8) to (10).

We now consider the power efficiency 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 Fig. 2. Decreasing the Fermi energy leads to higher nonlinear conductivity and lower linear conductivity by the graphene. Therefore, we obtain a higher power efficiency as the Fermi energy is decreased.

We see in Fig. 4(a) that for a weak input field of kV/cm, the power efficiency is almost the same in all three calculation schemes. It is only at low Fermi energies ( meV) that there is a noticeable difference in the power efficiency. As we increase the input field, we see in Figs. 4(b) and (c) that the effects of pump depletion, SFM, and XFM become very significant, particularly at low Fermi energies. For strong enough fields and low Fermi energies, our calculations for the undepleted pump approximation yield a non-physical power efficiency that is greater than parvin (). However, including pump depletion results in power efficiencies less than , as required. It is seen that the power efficiency decreases when SFM and XFM are included. For example, for Fermi energy of meV, when the input field is kV/cm, the power efficiency decreases from in the PPW with pump depletion to in the PPW with full calculation, while for an input field of kV/cm, the power efficiency decreases from to . It is therefore worth examining if we can modify the structure to decrease the phase mismatch introduced by SFM and XFM.

## Iv New Configuration: Mitigating self- and cross-phase modulation

In this section we define a new configuration of the PPW in order to deal with phase mismatching due to the SFM and XFM. We consider the waveguide shown in Fig. 5, where there are two different dielectric materials in the waveguide: cyclic polyolefin with refractive index of and phenol-formaldehyde resin with refractive index of . The graphene layer is located at midway between two metallic plates. The material is in the regions to and to , while the material in the region to . This new configuration allows us to control to some degree the phase matching between the third mode at and the first mode at . In the following, we optimize to obtain the best phase matching between the TE mode at and TE mode at , and thereby maximize the third harmonic generation and the power efficiency. Note that in all cases, the total plate separation is fixed at .

### iv.1 Phase Matching

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 mode at and TE mode at for a THz field in a waveguide in our original configuration and in our new configuration. As an example we choose . Note that the TE mode peaks at the centre of the waveguide, inside of the low-index () material, while the TE mode also has peaks inside of the high-index () 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 mode more than that of the TE 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 mode. This will yield a slight increase in the generated field for a given input power.

We now examine how the effective refractive index changes with the thickness of the region with lower refractive index (). In Fig. 7, we plot the linear effective refractive index difference, for a waveguide with plate separation of for different Fermi energies as a function of . In this calculation, we set and to zero. It is seen that perfect phase matching () occurs at two points for each value of as we increase . More importantly, we see that we can reduce by up to , which is more than enough to compensate for the effective index difference shown in Fig. 3 over the full range of Fermi energies, up to incident fields of at least kV/cm.

### iv.2 Power Efficiency

In Fig. 8 we plot the power efficiency of the waveguide in our new configuration as a function of for an incident field of kV/cm at Fermi energies of meV, meV, and meV. Decreasing the Fermi energy leads to higher power efficiency due to the increase in the nonlinear conductivity horng () and reduced loss due to a decrease in the linear conductivity Alnaib2015 (); cheng (). Note that when , 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 meV. Therefore, we expect to obtain the peaks in the efficiency close to the values of where 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 meV. We also note that for this small input field, the increase in the efficiency over the initial-configuration PPW is rather modest ().

In Fig. 9 we plot the power efficiency of the waveguide in the new configuration as a function of for a Fermi energy of meV for input fields of and 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 that are closer to the minimum in 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 meV, and kV/cm the power efficiency increases from to and for kV/cm it increases from to .

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 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 at high input fields. To understand this, consider Fig. 3, where we see that for meV, the index mismatch due to XFM and SFM is already at an input field of kV/cm. It is easy to see therefore that for the fields considered in Fig. 10, we will quickly reach index differences greater than 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 . 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 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 ( mm, mm and mm for Fermi energies of meV, meV, and 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 meV, we are able to obtain a power efficiency of at an input field of 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 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 Hafez (); Hadi ().

## 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 mode at and TE mode at 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 meV, where the nonlinearity is relatively modest, efficiencies of up to can be achieved for input field amplitudes of kV/cm. We therefore believe that our PPW system is an excellent platform to produce and examine harmonic generation in graphene.

Acknowledgements We thank the Natural Sciences
and Engineering Research Council of Canada and Queen’s University for financial support. The authors would like to thank Lukas Helt for useful discussions.

## Appendix A Dynamic equations of the linear and third-harmonic electric field

In this Appendix, we give the details of the derivation of our coupled-mode equations: Eqs. (18) and (19).

The LHS of Eq. (17) can be written as

(A.1) |

If we define ,

(A.2) |

where

(A.3) |

However, Eq. (A.3) is not valid at . According to the results obtained for the first derivative of the electric field at the graphene, shown in Fig. A.1, we must add a term to our calculations for the electric field at the graphene.

This term is given by , such that,

(A.4) |

To determine , we take integrate from to , for . Then, we obtain

(A.5) |

The electric field below and above the graphene is defined as

(A.6) | ||||

using Maxwell’s equations we have,

(A.7) | ||||

Thus, Eq. (A.7) can be written as

(A.8) | ||||

The current density at the graphene is given by

(A.9) |

where the current density at the graphene for mode is related to the field by

(A.10) |

where

(A.11) |

Thus, Eq. (A.10) can be written as

(A.12) |

The current density at the graphene is related to the linear conductivity by

(A.13) | ||||

where

(A.14) |

Equality of Eq. (A.12) and Eq. (A.13) gives

(A.15) |

Using Eq. (A.15) in Eq. (A.5), we obtain

(A.16) |

Thus, Eq. (A.4) becomes

(A.17) |

and Eq. (A.2) can be written as

(A.18) | ||||

We also have that, using the slowly-varying envelop approximation, where we neglect the second derivative of the envelope, that

(A.19) | ||||

If we use Eqs. (A.18) and (A.19) in Eq. (17) of the main text and retain only those terms that are maximally phase marched, then we obtain

(A.20) | ||||

For a lossy waveguide, the propagation constant is defined as . Using this removes the second term in Eq. (A.20). We now multiply Eq. (A.20) by and integrate over , using

(A.21) | ||||