Third harmonic generation from graphene lying on different substrates: Optical-phonon resonances and interference effects

Graphene is a nonlinear material which can be used as a saturable absorber, frequency mixer and frequency multiplier. We investigate the third harmonic generation from graphene lying on different substrates, consisting of a dielectric (dispersionless or polar), metalized or non-metalized on the back side. We show that the third harmonic intensity emitted from graphene lying on a substrate, can be increased by orders of magnitude as compared to the isolated graphene, due the LO-phonon resonances in a polar dielectric or due to the interference effects in the substrates metalized on the back side. In some frequency intervals, the presence of the polar dielectric substrate compensates the strongly decreasing with $\omega$ frequency dependence of the third-order conductivity of graphene making the response almost frequency independent.


I. INTRODUCTION
It was theoretically predicted 1 and then experimentally confirmed 2,3 that, due to the linear energy dispersion of quasiparticles in graphene -electrons and holes, this material should demonstrate strongly nonlinear electrodynamic properties. In recent years a great interest to the nonlinear electrodynamic phenomena in graphene arose, which was stimulated not only by purely academic interest but also by great perspectives which are opened up for using graphene in microwave-, terahertz-and opto-electronic devices. A variety of different nonlinear phenomena, such as harmonics generation 2,4-6 , frequency mixing 3,7,8 , saturable absorption [9][10][11] and so on, have been already experimentally observed; even more have been theoretically predicted  , for recent reviews see 40,41 .
The higher harmonics generation under a monochromatic irradiation of a nonlinear medium is one of the fundamental nonlinear phenomena. In uniform graphene the second harmonic cannot be observed due to its central symmetry; therefore the lowest higher harmonic which can be emitted from uniform graphene is the third one. An analytical quantum theory of the third-order response of graphene has been recently developed in Refs. [31][32][33] . It was shown that the largest up-conversion efficiency (ω → 3ω) is achieved in the low-frequency (microwave, terahertz), quasi-classical regimehω < ∼ E F where the inter-band electronic transitions between the valence and conduction bands of graphene can be neglected. At larger (infrared) frequencieshω ≃ E F the effect is quantitatively smaller but a number of sharp and narrow resonances related to the three-photon (3hω = 2E F ), two-photon (2hω = 2E F ) and one-photon (hω = 2E F ) transitions at the absorption edge have been predicted. The positions of these resonances depend on the electron/hole density in graphene (E F ∝ √ n s ) and hence on the gate voltage, therefore the infrared resonances in the third-order response function of graphene are very interesting for a potential fine electric control of the third-harmonic generation.
The theory [31][32][33] was developed for a single, isolated graphene layer. In most experiments, however, graphene lies on a dielectric substrate which is or can be covered by metal on the back side. Recently we have investigated how the presence of a substrate may influence the third harmonic intensity I 3ω emitted from such a structure 30 . In particular, we have shown that the third harmonic intensity, emitted from graphene lying on a dielectric metalized from the back side, can be larger than that from an isolated graphene layer by more than two orders of magnitude. The physical reason of such a huge enhancement of the up-conversion efficiency is the interference of both the incident wave (with the frequency ω) and the emitted wave (3ω) in the dielectric, see a discussion in Ref. 30 .
The results of Ref. 30 suppose that the dielectric constant ǫ (or the refractive index n = √ ǫ) of the dielectric substrate does not depend on the frequency, in particular, that n ω = n 3ω . In real materials the frequency dispersion of the dielectric permittivity ǫ(ω) is not negligible. On the other hand, the dielectric function ǫ(ω) itself may have strong resonances, related, e.g., to the excitation of the transverse and longitudinal phonons in polar dielectrics (for example, in SiO 2 , Al 2 O 3 ). These resonances may also influence the intensity of the third harmonics emitted from graphene placed on (polar) dielectric and possibly metalized substrates.
In this paper we investigate in details the influence of different types of substrates on the third-harmonic generation from graphene. In Section II we overview the applied theoretical approach. The results obtained are presented in Section III: we consider and analyze the structures AGA, AGDA, AGDM, AGPA, and AGPM, where the letters stand for air (A), graphene (G), dispersionless dielectric (D), polar dielectric (P), and metal (M). We show that the third harmonic intensity can be dramatically increased not only due to the interference of ω-and 3ω-waves in the back-metalized structures (the effect partly considered in Ref. 30 ) but also due to the LO-phonon resonances at the frequencies lying in the Reststrahlen-Band ω TO < ∼ ω < ∼ ω LO as well as above the LO-phonon frequency ω LO . We also study the response of graphene on a substrate at low (terahertz) frequencies. In Section IV we summarize our results and draw conclusions.

II. THEORETICAL APPROACH
We consider a structure shown in Figure 1. A linearly polarized (in the x-direction) electromagnetic wave with the frequency ω and intensity I ω propagates in the positive z-direction and is incident on graphene (G) lying on a substrate. The substrate consists of two layers, L1 and L2, which can be dielectric or metallic. Graphene produces the third harmonic which is emitted in the positive (forward) and negative (backward) z-direction.
To find the intensity of the third harmonic we solve the system of Maxwell equations with conventional boundary conditions: the continuity of E x and H y fields at all boundaries except the boundary z = 0 where the magnetic field has a jump The dielectric permittivity ǫ(ω, z) is constant (does not depend on z) inside each layer but may be a function of frequency. In the P-layer (polar dielectric) we assume that it has a single TO-phonon resonance so that where ω TO and ω LO are the frequencies of the transverse and longitudinal phonons and γ TO is a phenomenological relaxation rate of TO phonons; ǫ ∞ is the dielectric constant at infinite frequencies.
In the M-layer (metal) we describe the dielectric response by the Drude model where ω p is the plasma frequency in the metal and γ m = 1/τ m (τ m ) is the corresponding relaxation rate (time). The dielectric constants ǫ A and ǫ D of the A-and D-layers are frequency independent. The current j in Eq. (2) contains the first and the third harmonics, Here E x ω (0) is the complex amplitude of the electric field at the plane z = 0, c.c. means the complex conjugate, and the functions σ can be found, e.g. in Refs. 33,42 ; here Ω =hω/E F , Γ =hγ/E F =h/τ E F , γ (τ ) is the phenomenological relaxation rate (time) of electrons in graphene. The explicit analytical expressions for σ (3) αβγδ (ω 1 , ω 2 , ω 3 ) can be found in Refs. [31][32][33] ; they are very long [see, e.g. Eqs. (59)-(78) in Ref. 33 ] and we do not reproduce them here.
The problem is solved in two steps. First, we solve the linear response problem and calculate the transmission, reflection and absorption coefficients of the fundamental-harmonic wave (frequency ω), as well as the electric field E x ω (0) at the plane z = 0. Then we solve Maxwell equations for the 3ω-harmonic and calculate the intensities of the emitted third-harmonic radiation in the backward and forward directions. We present our results in terms of the parameter η (3) , which does not depend on the intensity of the incident radiation and describes the efficiency of the first-to thirdharmonic transformation. The quantity η (3) is measured in units (cm 2 /W ) 2 ; the value of η (3) = 10 −18 (cm 2 /W ) 2 (see Figures below) means that 1 MW/cm 2 of the input (ω) radiation produces 1 W/cm 2 of the output (3ω) signal. For specific results presented in the next Section III we use the following numerical parameters. In the dielectric layer we assume ǫ D = 4 which is close, e.g., to the dielectric constant (3.9) of SiO 2 . In polar dielectrics the typical values of the TO-and LO-phonon frequencies lie between ∼ 10 and 30 THz, see, e.g., Ref. 43 . We assume therefore ǫ ∞ = 1, ω TO /2π = 15 THz, ω LO /2π = 30 THz, and γ TO /2π = 0.2 THz. For the metallic layer we use parameters of gold 44,45 :hω p = 8.45 eV and τ m = 14 fs. The dielectric constant of air is ǫ A = 1.
In the next Section we present results obtained this way for isolated graphene and different layered structures.

A. Structure AGA
We begin our analysis from the isolated graphene layer, i.e. the structure AGA (air -graphene -air). In this case the third-harmonic intensity depends on two dimensionless parameters,hω/E F andh/τ E F , Ref. 33 . Figure 2(a) shows the parameter η (3) as a function of the electron density at two different values of the input-wave frequency and two values of the relaxation time τ , compare with Fig. 7(a) from Ref. 33 . One sees that, first, the up-conversion coefficient η (3) strongly depends on the input-wave frequency: the increase of the input-wave wavelength λ ω = 2πc/ω by a factor of 3 leads to the increase of η (3) by four orders of magnitude. Second, the density-dependence of η (3) contains three sharp resonances, athω = 2E F /3, E F , and 2E F , which corresponds to the three-photon (the largest resonance), two-photon and one-photon (the smallest one) inter-band transitions at the absorption edge. Third, the overall behavior of the function η (3) (n s ) is not very sensitive to the relaxation time τ , except the near-resonance regions. The resonances become very pronounced and sharp when the (inter-band) relaxation time gets bigger than ∼ 1 ps: one sees, for example, that the increase of τ by one order of magnitude leads to a more than (almost) two orders of magnitude increase of η (3) at λ ω = 10 µm (30 µm). Figure 2(b) shows the up-conversion coefficient as a function of frequency at a few values of n s and at τ = 1 ps. One sees a strong reduction of the effect with ω and sharp resonances related to the inter-band multi-photon transitions. The coefficient η (3) can be of order 10 −13 − 10 −12 (cm 2 /W ) 2 at the input-wave frequency f ≃ 1 THz and of order 10 −18 − 10 −17 (cm 2 /W ) 2 (at low densities and in out-of-resonance regions) at the input-wave frequency f ≃ 10 THz. Using the 2E F /3 resonances one can increase the last values by about two orders of magnitude if τ = 1 ps.

B. Structure AGDA
Now consider graphene lying on the dielectric substrate with the refractive index n ω = n 3ω = 2 and a dielectric thickness d (the structure AGDA). This situation was preliminary considered in Ref. 30 ; here we analyze it in more details, in particular, as a function of the input-wave frequency. We also need the results for the AGDA structure for the sake of comparison with the polar dielectric case (Section III D). Figure 3(a) shows the up-conversion efficiency η (3) as a function of the input-wave frequency at the electron density n s = 10 11 cm −2 and d = 12.5 µm. First, one sees that the third-harmonic intensity from the AGDA system is always smaller than from the isolated graphene layer (AGA). At several values of the frequency (≃ 6, 12, 18 THz) the value of η (3) for the AGDA structure coincides with that for the AGA structure. This is a consequence of the interference of both the first and third harmonics in the dielectric slab and these frequencies satisfy the conditions and i.e. the slab thickness d is a multiple of the half-wavelength of the first and third harmonics. Since we assume for the dielectric (D) that n ω = n 3ω = const, the second set of conditions (10) is a subset of the first one (9), so that the points where η AGA are determined by the condition (9). Second, the intensity of the third harmonic emitted in the forward direction is always higher than in the backward direction. The factor η AGDA for the third harmonic emitted back has additional oscillations with the local maxima determined by Eq. (10) (the interference of the third harmonic wave). These interference features are also seen in Figure 3(b) where we show the factor η (3) as a function of d at fixed values of τ , n s and f .
In Figure 3(a) we have chosen the electron density (n s ≈ 10 11 cm −2 ) and the dielectric thickness (d = 12.5 µm) so that one of the interference maxima (9), namely, the one with m 1 = 1, coincided with the largest inter-band resonance at otherwise, the large η (3) -value (≃ 7 × 10 −15 cm 4 /W 2 ) would be strongly suppressed by the destructive interference of waves reflected from the dielectric surfaces. The optimal condition for the observation of a resonant enhancement of the third harmonic in the AGDA structure is thus where n = n ω = n 3ω is assumed to be constant. The condition (12) relates the frequency, the electron density and the thickness of the dielectric slab. The number η (3) ≃ 7×10 −15 cm 4 /W 2 in our example corresponds to the emission of 7 W/cm 2 of the third-harmonic (18 THz) at the input wave (6 THz) intensity of 100 kW/cm 2 .

C. Structure AGDMA
Now consider what happens if the same structure that was analyzed in Section III B is covered by a thin metallic (Au) layer on the back side. Figure 4(a) shows the dependence of η (3) on the input-wave frequency at the same parameters as in Figure 3(a), in particular, at d = 12.5 µm; the only difference is that the backside of the dielectric is covered by a thin (0.2 µm) layer of gold. One sees that the spectrum of the up-conversion coefficient η (3) dramatically changes. First, since the metal reflects both the first and third harmonics, the emission of the 3ω-wave in the forward direction is suppressed by many orders of magnitude and can be completely neglected. Second, the emission of the 3ω-wave in the backward direction is strongly modified, and in certain frequency ranges, e.g. around the input-wave frequency f = 3, 9, and 15 THz the coefficient η (3) is substantially (by orders of magnitude) larger than in the AGA structure. For example, at f = 3 THz the value of η On the other hand, the resonance athω = 2E F /3 (≃ 6 THz), which is the largest inter-band resonance in the AGA structure, is completely suppressed (by more than ten orders of magnitude). This is explained, again, by the interference of waves in the dielectric, see 30 . In the presence of metal on its back side, however, the boundary condition at z = d (the boundary dielectric -metal) changes: the tangential electric field E x should vanish at this point and the interference maxima in the AGDMA structure are expected at i.e. the dielectric thickness is a quarter of wavelength plus an integer times a half-wavelength. In Figure 4(b) we plot the dependence of η AGDMA on the input-wave frequency f at two times smaller dielectric thickness which corresponds to the condition (13) with m 1 = 0. The resonances at 3 and 9 THz are suppressed while those at 6 THz and 18 THz (correspond to thehω = 2E F /3 and 2E F resonances) become substantially stronger and exceed the corresponding values of η AGA ≃ 1.9 × 10 −20 cm 4 /W 2 (the growth by a factor of ∼ 210). The value of η (3) ≃ 1.7 × 10 −12 cm 4 /W 2 , which can be achieved under the chosen parameters at the incident wave frequency of 6 THz, corresponds to the emission of 1.7 kW/cm 2 of the third-harmonic (18 THz) at the input wave intensity of 100 kW/cm 2 , compare with the last sentence in Section III B. Or, this means that the third-harmonic intensity of 1.7 W/cm 2 can be achieved already at the quite low intensity of the incident wave of only 10 kW/cm 2 .

D. Structure AGPA
The dielectric substrate with a dispersionless dielectric permittivity considered in the previous Section is not very common. Usually the dielectric function is frequency-dependent in the mid-IR range which is related with the opticalphonons poles. Real polar dielectrics typically have several TO-phonon poles, but in order to clarify the physics of the discussed phenomena we restrict ourselves here by the model (4) with a single TO-phonon pole, Figure 5(a).
The third-order response of the AGPA structure is much more complicated, as compared to the cases considered in the previous Sections, therefore it makes sense to consider first the linear response spectra. Figure 5(b)-(d) shows the transmission, reflection and absorption coefficients of the AGPA structure at d = 12.5 µm (the same dielectric thickness as in Figures 3(a) and 4(a)). The TO-phonon frequency is chosen to lie at 15 THz and the LO-phonon frequency -at 30 THz. Since ω LO /ω TO = 2 and ǫ ∞ = 1, the static dielectric constant equals 4 which coincides with ǫ D in the dielectric layer considered in Sections III B-III C. As seen from Figures 5(b)-(d), in the frequency window from ω TO to ω LO the AGPA structure reflects almost all the radiation and the transmission coefficient is close to zero (the Reststrahlen-Band effect). In the areas ω < ∼ ω TO and ω > ∼ ω LO the TRA-coefficients strongly oscillate which is related to the interference of the incident (ω) radiation in the dielectric slab, according to Eq. (9). In this equation now, the refractive index n ω = ǫ P (ω) depends on the frequency (strongly increases when ω → ω TO and tends to zero when ω LO ← ω), therefore the oscillation periods are not constant but vary approaching the TO-phonon frequency from the left and the LO-phonon frequency from the right. The absorption coefficient has two main maxima at the frequency slightly lower than ω TO and slightly higher than ω LO , with interference related oscillations. Now consider the third-order response. Having the goal to maximize the intensity of the third harmonic signal we can try to take advantage of combining the graphene resonance (11) with the optical phonon resonances at ω TO or ω LO . Therefore we consider four special cases.
1. The case 3ωres ≃ ωTO Figures 6(a,b) show the frequency dependence of η (3) for the AGA structure at τ = 1 ps, n s = 7.06858 × 10 10 cm −2 and the AGPA structure at the same values of τ and n s , and at d = 14.3 µm. The density n s is chosen so that the graphene resonance (11) lies at f res = ω res /2π ≃ 5 THz; then 3f res coincides with the TO-phonon frequency 15 THz. The chosen thickness d = 14.3 µm corresponds to one of the interference maxima, Figure 7(a), of the coefficient η (3) at the incident wave frequency f = 5.02 THz (the exact position of the resonance maximum in Figures 6(a,b)). Several interesting features are seen in Figures 6(a,b). First, the third-harmonic intensities, emitted in the forward and backward directions, are quite close to each other at f > ∼ 10 THz, but are substantially different at the lower frequencies, especially at 5 < ∼ f < ∼ 10 THz where the forward-emitted 3ω-intensity falls down by many orders of magnitude. This is explained by the fact that at 5 < ∼ f < ∼ 10 THz the frequency 3f lies in the Reststrahlen-Band, 15 < ∼ 3f < ∼ 30 THz, therefore the 3ω-wave produced by graphene does not penetrate into the polar dielectric and can therefore be emitted only in the backward direction. Second, the oscillations of η (3) at 10 < ∼ f < ∼ 15 THz are due to the interference of the input-wave frequency harmonic in the dielectric: the oscillation period becomes shorter when approaching the TO-phonon frequency from the left, similar to the behavior of TRA coefficients discussed above, Figure 5(b)-(d). It is interesting that a similar oscillating behavior of η (3) is also seen at the frequencies below 5 THz, for a detailed picture see Figure 6(b). This is due to the same interference effect but of the third harmonic: when f is approaching 5 THz, its third harmonic tends to 15 THz which is just the TO-phonon pole where the refractive index n 3ω diverges.
The influence of the polar dielectric on the third-harmonic intensity at the frequency f = f res ≃ 5 THz is seen in the detailed picture Figure 6(b) as well as in Figure 7(a). Although the thickness d = 14.3 µm was chosen to maximize the third order response, Figure 7(a), the coefficient η AGPA is by a factor of ∼ 50 smaller than η response at f ≃ f res .
However, at frequencies far from the resonance (11), at f > ∼ 20 THz in Figure 6(a), the polar-dielectric substrate does help to substantially increase the third-harmonic generation. For example, at the incident-wave frequency f ≃ 30 THz the coefficient η AGPA with the dielectric thickness d is illustrated in Figure 7(b). One sees that η (3) AGPA substantially grows when d varies from zero up to ≈ 5 µm; then it essentially saturates. The weak oscillations of the backward emitted third-harmonic intensity (the red curve in Figure 7(b)) is the consequence of the interference of the third-harmonic wave (90 THz in our example) in the dielectric slab: at 3f = 90 THz the polar dielectric has the refractive index slightly smaller than one, and the oscillation period in Figure 7(b) is well described by the interference formula (10).
Even more interesting feature seen in Figure 6(a) is that the third-harmonic intensity very weakly depends on the input-wave frequency f if it lies in the Reststrahlen-Band, 15 < ∼ f < ∼ 30 THz (the third-harmonic frequency 3f then lies between ∼ 45 and 90 THz). This fact can be very interesting for application. Indeed, the third-order conductivity of graphene σ xxxx (ω, ω, ω), responsible for the third-harmonic generation effect, falls down very quickly with ω. Therefore the coefficient η AGA of the isolated graphene quickly decreases with the frequency (by four orders of magnitude when f varies from 15 THz to 30 THz, see Figure 6(a)). The factor η (3) AGPA is almost frequency independent in this frequency range which means that the decrease of σ xxxx (ω, ω, ω) should be compensated by something else. To understand the reason of this effect we plot the frequency dependence of the electric field E x (z = 0) at the plane z = 0, as well as of the cube of this field, see Figure 8. One sees that the field grows in the interval from f TO to f LO , even stronger grows the cube of the field; therefore, the frequency dependence of the third-order conductivity is to a large extent compensated. This explains the almost flat frequency dependence of the third-harmonic intensity in the range 3f TO < ∼ 3f < ∼ 3f LO when the input-wave frequency lies in the Reststrahlen-Band f TO < ∼ f < ∼ f LO . When the incident wave frequency exceeds f LO , the coefficient η

The case ωres ≃ ωTO
Now consider the case when the resonance frequency ω res , Eq. (11), is close to the TO-phonon frequency. Figure  9 shows the dependence of the parameter η (3) of the AGPA structure on the input-wave frequency f at d = 10.54 µm, τ = 1 ps and n s = 0.636 × 10 12 cm −2 . The density of electrons is chosen so that the resonance frequency (11) corresponding to the conditionhω = 2E F /3 is close to the TO-phonon frequency f TO = 15 THz. The general behavior of η (3) is similar to that shown in Figure 6. Near the graphene resonance at 15 THz the third harmonic intensity is strongly (by many orders of magnitude) suppressed, but at the larger frequencies ( > ∼ 20 THz) the presence of the polar-dielectric substrate increases the isolated graphene response (by a factor ≃ 50 at 30 THz and by a factor of ≃ 10 at 35 THz).
As seen from Figures 6 and 9, in the cases considered in Sections III D 1 and III D 2, the third harmonic intensity is strongly suppressed by the polar dielectric substrate, if the incident wave frequency is close to the TO-phonon frequency, and substantially enhanced, if it is close to the LO-phonon frequency. Let us now consider what happens if the triple resonance frequency 3ω res is close to the LO-phonon frequency. 3. The case 3ωres ≃ ωLO Figure 10 demonstrates the input-frequency dependence of the parameter η (3) in the AGPA structure at the electron density n s = 0.2827 × 10 12 cm −2 . The resonance frequency (11) equals 10 THz in this case, then 3f res = f LO = 30 THz. One sees that at frequencies higher than ≃ 20 THz the behavior of η (3) is similar to the one observed in two previous cases: the presence on the substrate quite substantially increases the third harmonic intensity and flattens the frequency dependence of η (3) in the range ω TO < ∼ ω < ∼ ω LO . The strong suppression of η AGPA at the frequencies between 5 and 10 THz is due to the fact that the triple frequency lies in the Reststrahlen-Band, and the oscillations of η AGPA between 10 and 15 THz have the same origin as the oscillations of the TRA coefficients shown in Figures  5(b)-(d).
What is new in Figure 10 as compared to the previous cases is that the third-harmonic intensity (at the frequency 3f ≃ f LO ≃ 30 THz) is increased, when the input-wave frequency is close to the graphene resonance (11), f ≃ f res ≃ 10 THz. Figure 10(b) shows details of the frequency dependence of η (3) around this point. For the backward emitted third harmonic the factor η (3) is increased by a factor of ≃ 2.86 under the chosen conditions. This is due to the interference of the incident (ω) wave in the dielectric, as seen from Figure 10(c) where the dielectric thickness dependence is shown. The oscillation maxima of η (3) in this Figure correspond to the interference condition (9) where n ω ≈ 2.518 is the refractive index of the substrate at 10 THz. As seen from Figure 10(c) the substrates of certain thicknesses increase the effect, however, if the optimal conditions are not satisfied, the effect can be strongly suppressed as compared to the case d = 0.

The case ωres ≃ ωLO
Finally we consider the case when the graphene resonance (11) coincides with the LO-phonon resonance, Figure 11. In this case the density of electrons is about 2.5447 × 10 12 cm −2 which gives the resonance position at ≃ 30 THz. As seen from Figure 11(a) and especially Figure 11(b), the combination of two resonances helps to substantially increase the third-harmonic (90 THz) intensity if the input-wave frequency is close to the LO-phonon frequency and to the graphene resonance frequency f res (30 THz in our example). The amplification factor, under the chosen conditions, is about 52. The d-dependence of the up-conversion factor η (3) at f = 30 THz is shown in Figure 11(c). One sees that η (3) monotonously grows with d when d < ∼ 10 µm and then practically saturates. The polar dielectric substrate thus strongly increases the effect if d < ∼ 10 µm; at larger d the choice of the dielectric thickness is not very crucial. The intensity of the backward emitted third harmonic slightly oscillates with d while that of the forward emitted wave does not (see the discussion of Figure 7(b)). Notice that at the parameters of Figure 11(a) one of the interference maxima also coincides with the second graphene resonancehω = E F at f ≃ 45 THz. The third-harmonic response is also resonantly enhanced near this frequency, from η If the graphene resonance frequency (11) slightly deviates from ω LO the intensity of the third harmonic falls down drastically. This is illustrated in Figure 11(d) where the dependence of η AGPA on the electron density, and hence on the resonance frequency f res , is shown at the input-wave frequency f = 30 THz (which is equal to f LO in this case) and the polar dielectric thickness d = 14.82 µm. One sees that the resonance around the density n s = 2.54469 × 10 12 cm −2 is very sharp; small deviations from the resonance number leads to a strong suppression of the effect. Two other resonances at n s = 1.13 × 10 12 cm −2 and n s = 0.28 × 10 12 cm −2 correspond to the graphene resonanceshω = E F andhω = 2E F respectively. The sharpness of the resonances also depends on the relaxation time τ ; if it is smaller than 1 ps, the height of the resonances will be smaller, compare with Figure 2(a).

E. Structure AGPMA
Now we study how the thin metallic layer on the backside of the polar dielectric (the structure AGPMA) may influence the up-conversion efficiency η (3) as compared to the AGA structure. Figure 12(a) shows the frequency dependence of η (3) under the condition ω res = ω LO , which corresponds to the electron density n s = 2.54469 × 10 12 cm −2 . The polar dielectric thickness is chosen to be d = 11.31 µm which corresponds to the highest interferometric resonance at 30 THz, see Figure 12(b). One sees that, due to many intereferometric maxima, the spectrum of η AGPMA has many frequency bands in which η AGA ≈ 3.47 × 10 −19 (cm 2 /W) 2 . Like in the case of the AGDMA structures, the metalization of the back-side of the polar dielectric also leads to a strong enhancement, under certain conditions, of the third-harmonic generation from the AGPMA structure. The correct choice of the dielectric thickness is thus crucial for a successful experiment. Notice that at small d ( < ∼ 1 µm) the effect practically disappears due to the screening of the electric fields by metal.

F. Terahertz response of a AGDMA structure
Finally we show results obtained for the low-frequency (one to a few THz) response of the AGDMA structures. At such low frequencies the third-order response of graphene is not resonant, since at realistic electron densities the Fermi energy is typically larger thanhω. On the other hand, the absolute value of η (3) grows when the frequency decreases, therefore it makes sense to quantitatively investigate the role of the dielectric substrate at f ≃ 1 THz. At such low frequencies the dielectric permittivity of the polar dielectric can be considered to be frequency independent, therefore we ignore the difference between the AGPMA and AGDMA structures. Figure 13 shows the frequency dependence of the up-conversion parameter η AGDMA at different values of the dielectric thickness and two values of the graphene relaxation time, τ = 1 ps (Figure 13(a)) and τ = 0.1 ps ( Figure  13(b)). The back side of the dielectric is assumed to be metalized by gold of the thickness 0.4 µm, the refractive index n ω = n 3ω = 2 is considered to be frequency independent. One sees that the third-harmonic intensity can be more than two orders of magnitude larger than in the AGA structure and this enhancement can be achieved in at the input-wave frequency f ≃ 1 THz can be as large as ∼ 5 × 10 −11 (cm 2 /W) 2 in structures with the graphene relaxation time τ = 1 ps and ∼ 3 × 10 −12 (cm 2 /W) 2 at τ = 0.1 ps. This corresponds to the output signal intensity (at 3f ≃ 3 THz) of 50 and 3 mW/cm 2 respectively, at the input-wave intensity of only 1 kW/cm 2 .

IV. SUMMARY AND CONCLUSIONS
We have investigated the third harmonic generation effect from structures consisting of the nonlinear material graphene lying on different substrates. We have shown that, dependent on the type and thickness of the substrate and/or on a frequency interval, the third harmonic intensity can be both drastically suppressed and substantially enhanced as compared to the freely hanging graphene. The most essential growth of the third harmonic is shown to exist in the structures with the metalized back side and in structures with the polar dielectric serving as a substrate. In the first case the strong enhancement of the third harmonic is caused by the interference of the input (ω) and output (3ω) waves in the dielectric substrate, due to the reflection from the back-side metallic mirror. In the second case the growth of the third harmonic is caused by the interaction of either the input (ω) and output (3ω) waves with the LO phonons in the polar dielectric substrate. A correct choice of the substrate thickness is thus crucially  important for the optimization of the emitted third-harmonic intensity.
Another interesting effect that we have found is a flattening, in certain frequency intervals, of the frequency dependence of the graphene-on-a-substrate system, as compared to the very strong frequency dependence of the puregraphene response. We have shown that, if graphene lies on the polar-dielectric substrate, this strong frequency dependence can be compensated by the dispersion of the dielectric permittivity of the substrate, so that in the quite broad frequency interval between the TO and LO phonon frequencies, ω TO < ∼ ω < ∼ ω LO , the resulting intensity of the third harmonic, at 3ω TO < ∼ 3ω < ∼ 3ω LO , turns out almost frequency independent. To summarize, our work shows that the strongly nonlinear electrodynamic properties of graphene can be further increased by a proper design of the underlying substrate. Results of this work can be used in different terahertz and optoelectronic applications.