Nonlinear terahertz frequency conversion via graphene microribbon array

While it is believed that the universe is bathed in terahertz (THz) energy, and myriad applications in chemical and biological imaging, security technology, and astronomy are envisioned for the THz frequency range (roughly from 100 GHz − 10 THz), until relatively recently, it has remained one of the least tapped regions of the electromagnetic spectrum compared to its neighbors (the infrared and microwave bands) due to lack of efficient laboratory facilities such as emitters, detectors and so on. Since a large number of materials are transparent for THz waves, the development of THz technology and the implementation of THz components require materials with proper response to THz frequencies [1, 2]. Although metamaterials have been proven to offer a promising approach for this dilemma, the complexity in design and difficulty in the fabrication process have remained obstacles. Graphene, a recently isolated two-dimensional crystal of ${{\rm{sp}}}^{2} \mbox{-} {\rm bonded}$ carbon atoms, arranged in a honeycomb lattice, is a material of interest in this context. Over the last few years, whole ranges of experiments have revealed the fascinating electrical, optical and physical properties of graphene, mainly originated from its exceptional Dirac cone-type band structure [3–5].

While in the conventional plasmonics of noble metals there is no opportunity for efficient excitation of surface plasmon polaritons (SPPs) at the THz band and it is limited to the study of transverse magnetic (TM) SPPs at the near infrared and visible frequencies, graphene can support either TM or transverse electric (TE) SPPs depending on the sign of the imaginary part of its conductivity. TM plasmons in graphene offer a higher level of confinement, a lower degree of loss and therefore greater propagation length than their metal counterparts. TE plasmons are also specified by low loss even in room temperature but weak localization at the surface [6–9]. The strong coupling of light to SPPs and the consequent creation of strong localized electromagnetic fields and the high sensitivity of SPPs to the properties of the dielectric bounding medium, besides the fast response time of plasmonic excitations, which can strengthen the traditionally observed nonlinear phenomena in bulk semiconductors, make plasmonics a good platform for nonlinear studies [10]. The nonlinear plasmonics field includes both the nonlinear response of the dielectric environment and the active medium such as noble metals. Although the nonlinear response of metallic films and nanostructures is thoroughly investigated in the visible and near infrared frequencies, it cannot be extended to the THz regime since noble metals cannot effectively support SPPs at these frequencies. The limitations of noble metal based plasmonics can be effectively addressed by graphene, a material of atomic thickness. The capability of graphene to support SPPs at THz frequencies combined with its dynamically tunable, strong nonlinear response due to the particular linear energy dispersion of its charge carriers nominates graphene as a promising candidate for the relatively unexplored field of THz nonlinear plasmonics. The particular nonlinearity of graphene is theoretically predicted in [11–14] and experimentally confirmed in several cases including direct measurement of the third-order refractive index with the Z-scan technique [15], observation of four-wave mixing in graphene flakes [16], and third-harmonic generation (THG) [17].

Over the last few years, some aspects of graphene nonlinearity have been reported in the literature including its bistability [18], Kerr nonlinearity [19], self-induced regenerative oscillations [20], four-wave mixing [16, 20], second-harmonic generation (SHG) [21, 22], saturable absorption [23], THG [17] and soliton formation [24] with applications in SPP wave excitation [25, 26], phase shifting [27], filtering [28], switching [29], ultrafast pulsed lasers [30] and so on. In this contribution, we aim to study the enhancement of THG from graphene sheets under normal illumination of THz waves due to plasmon resonance. While in our recently published paper, we have investigated the effect of the quasi phase matching technique to enhance the efficiency of THG for propagating SPPs on graphene [31], here we show how patterning graphene as a microribbon array can considerably promote the efficiency of THG by several orders of magnitude due to plasmonic interactions, compared to infinite graphene sheets under normal illumination of THz waves. Although resonance-enhanced THG via azo-polymer resonant waveguide gratings [32], optical fibers [33], Ag nanoparticles [34] and film-coupled plasmonic strip resonators [35] has been studied in the visible and near infrared frequencies, the use of graphene due to its capability of supporting plasmonic modes in the THz band gap (providing opportunity to work below the diffraction limit), its large third-order nonlinearity (reducing power consumption for a given nonlinear response) and the gate voltage tunability of its first- and third-order conductivities introduces fascinating aspects to this topic, worthy of exploring. This is besides the atomic thickness and planar nature of graphene, which make it a promising candidate for miniaturized and highly integrated nonlinear THz plasmonic devices and circuits.

The current study presents the extension of the nonlinear plasmonic field to the THz regime. Here, for the first time to our knowledge, we reveal a noticeable enhancement in the conversion efficiency (CE) of dynamically tunable THG over the THz gap due to plasmonic resonance from a graphene microribbon array of nanometer thickness, which cannot be achieved by conventional noble metal based plasmonics.

Since graphene is a two-dimensional centrosymmetric crystal, SHG in it is prohibited [36]. However, in some particular cases, symmetry breaking achieved by fabricating nanoislands of graphene with finite size [37] or by an oblique incidence of radiation on infinite graphene sheets [21] can lead to SHG from graphene which is outside of the focus of this study.

Within random phase approximation and in the absence of a magneto-static bias field and spatial dispersion (where the conductivity components become wavenumber dependent), the graphene isotropic complex surface conductivity can be expressed as (assuming $\mathrm{exp}(j\omega t)$ time harmonic dependency) [38]:

$$\begin{eqnarray}&&\sigma (\omega ,{\mu }_{c},{\rm{\Gamma }},T)=-j\displaystyle \frac{{e}^{2}{K}_{B}T}{\pi {\hslash }^{2}(\omega -j{\tau }^{-1})}\\ &&\quad \times ({\mu }_{c}/{K}_{B}T+2\;\mathrm{ln}(\mathrm{exp}(-{\mu }_{c}/{K}_{B}T)+1))\\ &&\quad -\displaystyle \frac{j{e}^{2}}{4\pi \hslash }\;\mathrm{ln}\left[\displaystyle \frac{2|{\mu }_{c}|-(\omega -j{\tau }^{-1})\hslash }{2|{\mu }_{c}|+(\omega -j{\tau }^{-1})\hslash }\right]\end{eqnarray} \tag{ 1 }$$

where the first and second terms refer to the intraband and interband graphene conductivities originated from the intraband and interband transitions of electrons respectively. In this relation -e is the electron charge, $\hslash$ is the reduced Planck's constant, $\omega$ is the angular frequency, $\tau$ is the electron–phonon relaxation time, T is the temperature, ${K}_{B}$ is the Boltzmann constant, ${\mu }_{c}$ is the chemical potential representing the graphene Fermi level, and ${\epsilon }$ is the energy. The interband term in this relation is calculated approximately by considering the constraint of ${K}_{B}\;T\ll |{\mu }_{c}|,\;\;\hslash \omega .$ For frequencies $\hslash \omega \gtrsim 2{\mu }_{c},$ the graphene conductivity is dominated by the interband transition of electrons which is known as a major mechanism of plasmon loss (via excitation of electron–hole pairs) [6, 38]. Furthermore, in this frequency range the sensitivity of graphene conductivity to the Fermi level and frequency is negligible and has a universal value of ${\sigma }_{0}\approx {e}^{2}/4\hslash$ for $\hslash \omega \gg 2{\mu }_{c}$ corresponding to a noticeable absorption constant of 2.3% per graphene layer with atomic thickness [39–41]. On the other hand, for frequencies $\hslash \omega \lt 2{\mu }_{c},$ the Pauli principle will block interband transitions (in p-doped graphene there is no electron available for the interband transition and in n-doped graphene electron states in resonance in the conduction band are occupied) and the intraband transition of electrons (free carrier absorption) plays a significant role. In this frequency range, which is considered as a low-loss regime (the real part of conductivity is much smaller than the imaginary part) for graphene operation, the graphene conductivity reveals great tunability with the Fermi level [38, 42].

The nonlinear conductivity of graphene is also defined as [11, 12, 43]:

$$\begin{eqnarray}&&{\sigma }_{3}(\omega )=j\displaystyle \frac{3}{32}\displaystyle \frac{{e}^{2}}{\pi {\hslash }^{2}}\;\displaystyle \frac{{(e{v}_{F})}^{2}}{\mu {\omega }^{3}}(1-j{\alpha }_{T})\end{eqnarray} \tag{ 2 }$$

In the low-frequency limit, where $\hslash \omega \lt 2{\mu }_{c},$ the interband transition of electrons occurs at a negligible rate compared to the intraband one and the third-order conductivity of graphene is proportional to $1/{\omega }^{3}.$ In contrast, for higher frequencies where $\hslash \omega \geqslant 2{\mu }_{c},$ the interband transition of electrons plays the dominant role and the frequency dependency of third-order conductivity would be proportional to $1/{\omega }^{4},$ according to quantum mechanical analysis [16, 28]. The coefficient ${\alpha }_{T}$ in (2) represents the effect of two-photon absorption (TPA) in graphene which in recent experiments is shown to be $\approx 0.1$ [20, 43]. Since TPA is essentially an interband mechanism and we work in low THz frequency on moderately doped graphene samples in which the intraband transition of electrons has the dominant role besides the smallness of the real part of the graphene nonlinear conductivity compared to its imaginary part and also for numerical simulation simplicity, we ignore the effect of TPA in graphene in our study.

Our simulations are performed in room temperature $T=300\;{\rm{K}}$ and the electron–phonon relaxation time is calculated as $\tau ={\mu }_{c}\mu /e{v}_{F}^{2}$ [44] where ${v}_{F}={10}^{6}\;{\rm{m}}\;{{\rm{s}}}^{-1}$ is the Fermi velocity and $\mu$ is the mobility of charge carriers taken in the order of $100\;000\;{{\rm{cm}}}^{2}/{\rm{V}}.{\rm{s}}$ in our simulation corresponding to $\tau =3\;{\rm{ps}}$ for ${\mu }_{c}=0.3\;{\rm{eV}}$ which is a rather large value for graphene samples on a substrate in room temperature (while an ordinary one for suspended graphene sheets in low temperatures and also reported in several works for fabrication under special conditions for the lowest level of impurities, with multilayer structures or with special types of dielectric environments) [45–47]. This is only by considering the graphene operation in the nonlinear regime in our study. Actually for lower typical values of relaxation time in the order of 0.1 ps graphene loss would be noticeable and to work in the nonlinear regime, high and possibly unachievable light intensities may be required. In this regard, for nonlinear purposes, low-loss and high-quality samples are preferred.

To investigate the THG capability of graphene, schematically illustrated in figure 1(a), and calculate the power transmission, we employ a subcell finite-difference time-domain (FDTD) numerical technique which can precisely model the atomic thickness of graphene. Taking the assumption that the graphene sheet is infinite and invariant in the y direction, the simulation can be performed in a two-dimensional case, according to figure 1(b), without loss of accuracy. Certainly, more accurate results by accounting the effect of finite length in this direction can be accomplished by three-dimensional subcell FDTD but at the cost of a large amount of computer memory and computing time. The periodic boundary condition is enforced to the x direction to simulate the periodicity and infinity and the convolutional perfectly matched layer (CPML) is employed in the z direction to dissipate the outgoing waves [48–51]. The number of CPML cells is set to be 80, to ensure the complete absorption of power at the boundaries. The FDTD grid size is chosen as ${\rm{\Delta }}x={\rm{\Delta }}z=0.1\;\mu {\rm{m}}$ [52] (the results for a smaller grid converge to the analytic ones). The time step is also set following the Courant stability condition ${\boldsymbol{\Delta }}{\boldsymbol{t}}={\bf{0}}.85/({\boldsymbol{c}}\sqrt{{\rm{\Delta }}{{\boldsymbol{x}}}^{-2}+{\rm{\Delta }}{{\boldsymbol{z}}}^{-2}})$ [51], where c is the speed of light in the free space (according to our numerous simulations, our developed code is stable). We recall that by working in the frequencies $\omega \ll 2{\mu }_{c}/\hslash$ only the intraband transitions of electrons would be relevant and only this term of graphene linear and nonlinear conductivities should be included in the FDTD code.

Zoom In Zoom Out Reset image size

By defining a volume conductivity for a ${\rm{\Delta }}-$thick graphene layer as ${\sigma }_{v}\equiv \sigma /{\rm{\Delta }},$ the volume current density and thus Ampere's law can be written as ${\bf{J}}={\sigma }_{v}{\bf{E}}$ and ${\boldsymbol{\nabla }}\times {\bf{H}}={\bf{J}}+j\omega {\varepsilon }_{0}{\bf{E}}=({\sigma }_{v}+j\omega {\varepsilon }_{0}){\bf{E}}$ respectively. So we can ascribe an equivalent complex permittivity of ${\varepsilon }_{\mathrm{eq}}\equiv {\varepsilon }_{0}+{\sigma }_{i}/\omega {\rm{\Delta }}-j{\sigma }_{r}/\omega {\rm{\Delta }}\;\;$ and thus an equivalent effective linear susceptibility of ${\chi }_{{\rm{e}}{\rm{f}}{\rm{f}},g}^{(1)}=-\tfrac{j\sigma }{\omega {\varepsilon }_{0}{\rm{\Delta }}}\;$ to the ${\rm{\Delta }}-$thick graphene layer where ${\sigma }_{r}$ and ${\sigma }_{i}$ denote the real and imaginary parts of graphene conductivity [53, 54]. In the same way, an equivalent nonlinear susceptibility can be designated to graphene defined as ${\chi }_{{\rm{e}}{\rm{f}}{\rm{f}},g}^{(3)}=-\tfrac{j{\sigma }_{3}}{\omega {\varepsilon }_{0}{\rm{\Delta }}}$ [55, 56].

Although it seems that the more appropriate implementation of graphene as a two-dimensional material in the FDTD code would be by its linear and nonlinear conductivities (and so by its linear and nonlinear currents), the bulk representation of graphene (certainly of finite thickness) by its linear and nonlinear susceptibilities (and so by its linear and nonlinear polarizations) in this study is only due to the presence of a Kerr-nonlinear medium in the graphene environment which must be modeled by polarization.

By working in the frequency range where only the intraband term of graphene conductivity is relevant, the effective linear susceptibility can be written as [50, 51]:

$$\begin{eqnarray}&&{\chi }_{{\rm{e}}{\rm{f}}{\rm{f}},g}^{(1)}=F/(j\omega Z-{\omega }^{2}G)\end{eqnarray} \tag{ 3 }$$

$$\begin{eqnarray}&&F=\displaystyle \frac{{e}^{2}{K}_{B}T}{a{\varepsilon }_{0}{\rm{\Delta }}}({\mu }_{c}/{K}_{B}T+2\;\mathrm{ln}(\mathrm{exp}(-{\mu }_{c}/{K}_{B}T)+1))\end{eqnarray} \tag{ 4 }$$

$$\begin{eqnarray}&&G=\pi {\hslash }^{2}\end{eqnarray} \tag{ 5 }$$

$$\begin{eqnarray}&&Z=\pi {\hslash }^{2}{\tau }^{-1}\end{eqnarray} \tag{ 6 }$$

where a is the portion of the cell occupied by graphene which in our simulation is taken to be 0.01 (the graphene thickness is taken to be 1 nm). However it is shown that inasmuch as the graphene thickness is sufficiently small compared to the working wavelength, this special choice is not of critical importance [53].

To simulate the spectral response of the system and evaluate its THG capability, Maxwell's equations are solved by using an FDTD algorithm. For this purpose, at first electric flux density D is calculated using a discrete form of Ampere's law (in the time and space domains) $\;{\boldsymbol{\nabla }}\times {\bf{H}}=\tfrac{\partial }{\partial {\boldsymbol{t}}}{\bf{D}}$ where H is the magnetic field vector. The key point to model graphene in our FDTD code is in the next step by use of the following constitutive relation [48]:

$$\begin{eqnarray}&&{\bf{D}}={\varepsilon }_{0}{\varepsilon }_{\infty }{\bf{E}}+{{\bf{P}}}_{{\rm{L}}}+{{\bf{P}}}_{{\rm{NL}}}^{{\rm{D}}}+{{\bf{P}}}_{{\rm{NL}}}^{{\rm{G}}}\end{eqnarray} \tag{ 7 }$$

in which we have assigned the linear polarization of ${{\bf{P}}}_{{\rm{L}}}$ and the nonlinear polarization of ${{\bf{P}}}_{{\rm{NL}}}^{{\rm{G}}}\;$ to characterize the linear and nonlinear dispersive nature of graphene and the nonlinear polarizations of $\;{{\bf{P}}}_{{\rm{NL}}}^{{\rm{D}}}\;$ to describe the Kerr-effect of the graphene bounding medium. The relations used to incorporate ${{\bf{P}}}_{{\rm{L}}}$ and ${{\bf{P}}}_{{\rm{NL}}}^{{\rm{D}}}$ in our FDTD code can be found in [50, 51] and ${{\bf{P}}}_{{\rm{NL}}}^{{\rm{G}}}$ is implemented as a fourth-order time-derivative (considering the $1/{\omega }^{4}\;$ frequency dependency of ${\chi }_{eff,g}^{(3)})$ following the procedure reported in [48, 57]. A Newton iteration technique is applied to find the electric field from equation (7). Finally, the magnetic field can be computed by a discretized expression of Faraday's law ${\boldsymbol{\nabla }}\times {\bf{E}}=-{\boldsymbol{\mu }}\tfrac{{\boldsymbol{\partial }}}{{\boldsymbol{\partial }}{\boldsymbol{t}}}{\bf{H}}$ where ${\boldsymbol{\mu }}$ is the magnetic permeability.

Throughout this work, we employ light of TM polarization in normal incidence (the electric field is in the x direction and along the width of the ribbon). This condition is necessary to excite TM SPPs in the following explored resonant structures (if the incident electric field is oriented along the length of the ribbon instead of the width, no resonance is observed since the ribbon is too long to support any resonance in our working frequency range besides that by considering the sign of the imaginary part of graphene conductivity; for the chosen graphene Fermi levels in this study and also our working frequency range, the graphene can only support TM SPPs). To study THG, a cosine-modulated Gaussian wave packet centered at the frequency of the fundamental wave (FW), narrow in the frequency domain and thus long-standing in the time domain, is launched to the structure (along a line in the x direction on one side of the structure) which requires a rather large computational box. Certainly to calculate the transmittance, reflectance and absorbance spectra over a wide spectral range, the illuminated Gaussian wave should be sufficiently broadband (wide in the frequency domain and narrow in the time domain). Our excitation source can be expressed mathematically as ${\boldsymbol{Exc}}({\boldsymbol{t}})={\boldsymbol{A}}\;\mathrm{cos}(2{\boldsymbol{\pi }}{{\boldsymbol{f}}}_{{\boldsymbol{c}}}{\boldsymbol{t}})\mathrm{exp}(-{({\boldsymbol{t}}-{{\boldsymbol{t}}}_{0})}^{2}/{{\boldsymbol{\tau }}}^{2})$ where ${\boldsymbol{t}}$ is the time in the FDTD loop, ${\boldsymbol{A}}$ is the amplitude, ${{\boldsymbol{f}}}_{{\boldsymbol{c}}}$ is the center frequency, ${{\boldsymbol{t}}}_{0}$ is the time delay and ${\boldsymbol{\tau }}$ determines the pulse width in the time and frequency domains.

To compute the transmitted power $({{\bf{P}}}_{{\bf{t}}}),$ the structure (infinite graphene sheet or ribbon array) is excited from one side (e.g. the top of the structure) and the FDTD calculated fields are recorded over a line (along the x direction) on the other side (e.g. the bottom of the structure). By taking the fast Fourier transform from the recorded fields and integrating the pointing vector over the sampling line, ${{\bf{P}}}_{{\bf{t}}}$ can be calculated. Removal of the structure (a situation in which the medium in the simulation region is homogenous and no reflection occurs) is done to calculate the incident power $({{\bf{P}}}_{{\bf{i}}{\bf{n}}}).$ The normalized transmitted power is defined as ${\bf{T}}={{\bf{P}}}_{{\bf{t}}}/{{\bf{P}}}_{{\bf{i}}{\bf{n}}}.$ Sampling for the reflected power $({{\bf{P}}}_{{\bf{r}}})$ is performed on the same side of the excitation and the difference in calculated power in the presence and absence of the structure gives ${{\bf{P}}}_{{\bf{r}}}$.

All of our simulations are performed by exploiting the C++ programming language and MATLAB software on a P4 Core(TM)2 Quad CPU 3 GHz with 8 GB RAM. In the THG study, more than 300 000 time steps (corresponding to 260 min running time) are considered, and in the calculation of transmittance, reflectance and absorbance spectra, around 80 000 time steps (corresponding to 70 min running time) are taken into account.

To validate the FDTD code we developed, the normalized transmitted (T), reflected (R) and absorbed (A) powers for the normally incident light wave upon an extended graphene sheet are computed and the results are compared with the analytic solutions defined as [49]:

$$\begin{eqnarray}&&T=\displaystyle \frac{{\eta }_{1}}{{\eta }_{2}}| 2{\eta }_{2}/({\eta }_{2}+{\eta }_{1}+\sigma {\eta }_{1}{\eta }_{2}){| }^{2}\end{eqnarray} \tag{ 8 }$$

$$\begin{eqnarray}&&R=\;| \displaystyle \frac{{\eta }_{2}-{\eta }_{1}-\sigma {\eta }_{1}{\eta }_{2}}{{\eta }_{2}+{\eta }_{1}+\sigma {\eta }_{1}{\eta }_{2}}{| }^{2}\end{eqnarray} \tag{ 9 }$$

$$\begin{eqnarray}&&A=1-T-R\end{eqnarray} \tag{ 10 }$$

where ${\eta }_{n}=\sqrt{{\mu }_{n}/{\varepsilon }_{n}\;}$ and $n=1,\;2$ refer to the regions above and below the graphene sheet. The result of this comparison, which is depicted in figure 2, confirms the accuracy of our code. (The effect of the finite thickness of the substrate and possibly related Fabry–Perot (FP) resonances are ignored here).

Zoom In Zoom Out Reset image size

We have also examined our code in the case of SPP wave propagation on a graphene sheet. For this purpose the frequency and Fermi level dependency of the wavelength of the propagating TM SPP wave are computed by the dispersion relation according to [6]

$$\begin{eqnarray}&&\displaystyle \frac{{\varepsilon }_{r1}}{\sqrt{{\beta }^{2}-{\varepsilon }_{r1}{\omega }^{2}/{c}^{2}}}+\displaystyle \frac{{\varepsilon }_{r2}}{\sqrt{{\beta }^{2}-{\varepsilon }_{r2}{\omega }^{2}/{c}^{2}}}=\displaystyle \frac{j\;\sigma (\omega )}{\omega {\varepsilon }_{0}}\end{eqnarray} \tag{ 11 }$$

and are compared with the ones obtained from our FDTD code where ${\varepsilon }_{r1}$ and ${\varepsilon }_{r2}$ are the relative dielectric constant of the surrounding medium, $\beta$ is the plasmon wave vector and $c$ is the free space light speed. The results show very good agreement as depicted in figures 3(a) and (b) respectively. The inset in figure 3(b) shows the spatial distribution of ${E}_{z}$ in the x–z plane in the steady-state case for a suspended graphene sheet at the frequency of 5 THz and ${\mu }_{c}=0.2\;\;{\rm{eV}}$.

Zoom In Zoom Out Reset image size

Figure 4 shows the normalized output power (the output power is normalized to the incident power at the FW) for the case depicted in figure 1. Here the intensity of the FW is ${\rm{0.1}}\;{\rm{MW}}\;{{\rm{cm}}}^{-2}$ and the graphene Fermi level is set as ${\mu }_{c}=0.1\;\;{\rm{eV}}.$ The first peak is attributed to the excitation of the structure with a narrowband, cosine-modulated Gaussian wave packet for which the center frequency is taken to be the frequency of the FW and here is set at 0.5 THz. While for low-intensity illumination there is no second peak in the output spectrum, by increasing the illumination intensity a peak appears at the place of the third harmonic (TH) (1.5 THz) due to third-order nonlinearity in the graphene and its substrate (in fact there is no direct excitation at TH). The substrate chosen is ${{\rm{As}}}_{2}{{\rm{Se}}}_{3}$ glass as a Kerr-nonlinear medium, with a linear refractive index of 2.86 and a nonlinear refractive index of $\;\;{n}_{2}=11.4\times {10}^{-14}\;{{\rm{cm}}}^{2}/{\rm{W}}$ in the THz regime [58]. This Kerr-nonlinear substrate is taken to be isotropic, homogenous and dispersionless over our frequency range. In what follows, to have a concept of THG strength, we evaluate the CE defined as ${\rm{CE}}=10\;{\mathrm{log}}_{10}({{\rm{P}}}_{{\rm{out}}\_{\rm{TH}}}/{{\rm{P}}}_{{\rm{in}}\_{\rm{FW}}})$ where ${{\rm{P}}}_{{\rm{out}}\_{\rm{TH}}}$ and ${P}_{{\rm{in}}\_{\rm{FW}}}\;$ stand for the produced TH power by graphene and the incident pump power at the FW respectively. The CE for the result presented in figure 4 is obtained around −80 dB.

Zoom In Zoom Out Reset image size

THG from the graphene sheet under normal illumination is mainly affected by three factors. According to (2) stronger nonlinear response is expected for graphene layers of lower Fermi level. This is mainly due to the fact that for operations in the nonlinear regime, the total energy gained by electrons in graphene under high intensity illumination should be high enough compared to their average Fermi energy. The intensity of the FW $({I}_{{\rm{FW}}})\;$ is also a key factor that largely contributes to THG efficiency. While for lower FW intensities the response of the system is linear, for higher ones it enters the nonlinear regime and higher values of THG efficiency are obtained. Finally, the frequency of the FW $({f}_{FW})$ is considered as an influencing factor. According to (2), the graphene third-order conductivity is higher at lower frequencies so a higher TH level is expected for operation in this region. Figure 5 shows how CE varies with these three factors.

Zoom In Zoom Out Reset image size

While for lightly doped graphene samples under illumination of low-frequency THz waves the contribution of graphene in THG is dominant and THG from a Kerr-nonlinear substrate is negligible, for heavily doped graphene samples under high-frequency THz wave illumination, the graphene nonlinear response becomes weak and comparable with that of the substrate. Certainly, this fact cannot be overlooked that to get a nonlinear response from graphene samples of a high Fermi level with higher frequency of the FW, the intensity of illuminating THz waves should be increased which in turn enhances the contribution of the substrate in THG. This phenomenon is explored in figure 6.

Zoom In Zoom Out Reset image size

In what follows, we investigate how plasmonic interactions by creating a strong localized electromagnetic field may help to improve the efficiency of THG from graphene under normal illumination of THz waves. For this purpose, we investigate a periodic array of graphene microribbons, as depicted in figure 7.

Zoom In Zoom Out Reset image size

The resonance condition in the normal illumination of the graphene microribbon array, depicted in figure 7 (as far as the coupling strength between the ribbons would be negligible, especially for arrays with periodicity sufficiently greater than the ribbon widths), is achieved following the excitation of SPPs in each ribbon where a constructive interference occurs for an SPP round trip [59]:

$$\begin{eqnarray}&&2{\rm{Re}}\{\beta \}{\rm{w}}+2{\varphi }_{r}=2m\pi \end{eqnarray} \tag{ 12 }$$

where $\beta$ is the propagation constant of the SPPs, m is an integer (the first resonance can be determined by placing $m=0)$, Re indicates the real part, ${\rm{w}}$ is the ribbon width and ${\varphi }_{r}$ is the phase induced upon reflection at the ends of the ribbon, commonly taken to be $-\pi$ based on the assumption that the field vanishes at the ribbon edge. However, considering the complex excitation of deeply evanescent modes near the ribbon edges and imposing the continuity of electric and magnetic fields leads to a nontrivial phase of about $-3\pi /4$ which is essentially independent of the intrinsic properties of graphene (Fermi level, temperature....), the wavelength of operation, and the dielectric environment, and strongly affects the resonant frequencies [50, 59].

Figure 8 shows the transmittance (T), reflectance (R) and absorbance (A) spectra for the schematically depicted case in figure 7 under normal illumination of low-intensity THz waves with ${\mu }_{c}=0.4\;{\rm{eV}},$ $w=\;4.2\;\mu {\rm{m}}\;$ and $p=15\;\mu {\rm{m}},$ and a silicon substrate of ${\varepsilon }_{r}=11.9$ at THz frequencies [60]. Consistent with our expectation, a notch in the transmittance spectrum and a peak in the reflectance and absorbance spectra at the resonant frequency are created due to the excitation of SPPs. Charge accumulation along the edge of the graphene ribbon may affect its local conductivity which for the sake of numerical simulation simplicity is not taken into account in our study; however, including this effect may lead to more realistic results [61].

Zoom In Zoom Out Reset image size

Figure 9(a) shows the FDTD calculated transmittance spectra for different values of the graphene Fermi level with $w=5.4\;\mu {\rm{m}},$ $p=15\;\mu {\rm{m}}$ and a silicon substrate of ${\varepsilon }_{r}=11.9.$ Figure 9(b) illustrates that the resonant frequency undergoes a blue shift with an increase in Fermi level due to a reduction in the $Re(\beta )$ of the excited SPPs following equation (12) [50, 51]. The resonant frequencies obtained from our FDTD code reveal very good agreement with those obtained from equation (12). The quality factor $(Q)$ of resonance shows an increase with the graphene Fermi level as depicted in figure 9(c). This is mainly attributed to the higher resonant frequency and also an increase in the figure of merit (FOM) of propagating SPPs on graphene defined as $\;{\rm{FOM}}=Re(\beta )/Im(\beta )$ [31]. The reduction in minimum transmittance at the resonant frequency with an increase in graphene Fermi level is also attributed to the higher number of charge carriers contributed in the plasmonic oscillation [62]. Besides that, the minimum transmittance is strongly affected by the value taken for the mobility of the charge carriers. In fact, for lower-quality graphene samples (with lower mobility of charge carriers), the loss is more noticeable and the minimum transmittance (the quality factor of resonance) increases (decreases).

Zoom In Zoom Out Reset image size

It is well known that a change in graphene Fermi level can be achieved by chemical doping or applying a gate voltage between a Au electrode on graphene and a highly doped silicon substrate which can act as a bottom electrode while an insulating layer, typically silicon dioxide ${{\rm{SiO}}}_{2},$ is positioned between the graphene and silicon substrate. However, in this study, we have ignored the finite thickness of the substrate and related FP resonances. Certainly, the coincidence of the frequency of FP and plasmon resonances and the proper positioning of the graphene microribbon array at the place of maximum field enhancement in the FP cavity (or even the coincidence of FP resonance with the TH) may help to increase the TH level which is not the subject of the current study.

Following equation (12) a red shift in the resonant frequency by increasing the ribbon width is expected. Figure 10(a) shows the transmittance spectra for different values of the ribbon width while the graphene Fermi level and period of the array are kept constant at ${\mu }_{c}=0.4\;{\rm{eV}}$ and $p=9\;\mu {\rm{m}}$ respectively. The relation reported in equation (12) is valid for single ribbons or arrays of periodicity much greater than the ribbon width. By matching the ribbon width to the period of the array due to an increase in coupling strength between adjacent ribbons, the resonant frequency undergoes a red shift compared to a single ribbon and the quality factor of resonance also shows a reduction. The reduction in quality factor with an increase in ribbon width may also be attributed to the lower resonant frequency besides the lower FOM of the SPP mode at these frequencies. Figures 10(b) and (c) respectively explore the variation of resonant frequency and the quality factor with ribbon width.

Zoom In Zoom Out Reset image size

We also investigate how the period of the array can affect the resonant frequency and quality factor of resonance without any change in the graphene Fermi level and ribbon width. Figure 11(a) depicts the transmittance spectra for different values of the period at ${\mu }_{c}=0.3\;{\rm{eV}},$ $w=1.8\;\mu {\rm{m}}$ and ${\varepsilon }_{r}=11.9.$ The variation of resonant frequency and the quality factor with the array period is also illustrated in figures 11(b) and (c) respectively. Consistent with our discussion above, by increasing the period (and thus the inter-ribbon space) the resonant frequency approaches the value of a single ribbon and the quality factor of resonance also increases due to the higher resonant frequency and greater confinement of the SPP mode.

Zoom In Zoom Out Reset image size

Now let us investigate how the resonance in the above investigated graphene microribbon array can enhance TH strength. We expect in the case of the coincidence of the resonant frequency with the FW frequency, the localization of the electric field around the graphene ribbon, which boosts the nonlinear response of the structure arising from both the graphene and the nonlinear substrate, leads to considerable improvement in CE. To show the effectiveness of patterning graphene as a microribbon array and the influence of different parameters on the improvement of the THG CE, we define ${\rm{CEI}}({\rm{dB}})={\rm{CERA}}({\rm{dB}})-{\rm{CEIN}}({\rm{dB}})$ which compares the CE obtained by the graphene microribbon array denoted by CERA with a reference case of an infinite graphene sheet denoted by CEIN (schematically depicted in the insets of figures 12(a) and (b) respectively). To study each parameter, other parameters involved in both cases (including the intensity and frequency of the illuminating wave and other structural and graphene parameters) are taken to be the same.

Zoom In Zoom Out Reset image size

Figure 12 explores this fact. The transmittance spectrum for a graphene microribbon array of ${\mu }_{c}=0.42\;{\rm{eV}},$ $w=2\;\mu {\rm{m}}$ and $p=10\;\mu {\rm{m}}$ on a substrate made of ${{\rm{As}}}_{2}{{\rm{Se}}}_{3}$ glass is depicted in figure 12(a) which shows a resonance around ${f}_{0}=4.3\;{\rm{THz}}.$ This resonant frequency is chosen as ${f}_{{\rm{FW}}}.$ For the infinite graphene sheet a very low TH level with noiselike behavior for ${I}_{{\rm{FW}}}=0.1\;\;{\rm{MW}}\;{{\rm{cm}}}^{-2}$ is observed according to figure 12(b), by patterning it as a microribbon array; interestingly, the TH level at the same ${{\boldsymbol{I}}}_{{\bf{FW}}}\;$ and ${{\boldsymbol{f}}}_{{\bf{FW}}}\;$ is increased by around six orders of magnitude.

According to our discussion above, increasing the graphene Fermi level and period of the array and reducing the ribbon width lead to resonances of a higher quality factor and hence a field enhancement on the graphene which is expected to improve the efficiency of THG. This expectation is confirmed by our FDTD simulation results depicted in figures 13(a)–(c) respectively.

Zoom In Zoom Out Reset image size

Certainly, it cannot be overlooked that by increasing the graphene Fermi level which is accompanied by a corresponding increase in the resonant frequency of the structure, a reduction in the nonlinear susceptibility of the graphene is observed. Therefore, even a higher quality factor (and hence a greater field enhancement) of the resonance obtained from the structures of graphene samples of higher Fermi levels may not lead to a continuous and linear increase in CEI as our results in figure 13(a) indicate.

Very clean and high-quality graphene samples which possess carriers of high mobility (corresponding to long relaxation time) are preferred for our study as they lead to stronger resonances with higher quality factors and reduce the required energy for nonlinear response. Figure 14(a) presents the transmittance spectra for different values of carrier mobility and figure 14(b) explores how patterning graphene samples with higher carrier mobility as a microribbon array can more noticeably improve the CE of THG. In the case of using conventionally produced graphene samples of lower mobility (lower than 20 000 cm² V^–1 s^–1) [63–65] depending on the power consumption criteria and experimental facilities, one should consider the optimization of the structure for resonances of higher quality factors to improve the nonlinear response of the system. Besides that although the commonly used chemical vapor deposition (CVD) grown graphene samples have lower mobility than exfoliated graphene, mainly attributed to the transfer process in the CVD technique, this challenge can be addressed by advanced transfer methods [66].

Zoom In Zoom Out Reset image size

Figure 15 shows CEI as a function of FW frequency. As expected, when the frequency of the FW is tuned around the plasmon resonant frequency, a noticeable enhancement in CEI is observed. The frequency of the maximum CEI shows a small red shift in the order of 35 GHz with respect to the plasmon resonant frequency (obtained from our linear simulation) which is attributed to the modification of graphene conductivity by high-intensity illuminating waves [28].

Zoom In Zoom Out Reset image size

Regarding the asymmetric behavior of CEI versus the frequency of the FW around the resonant condition in figure 15, one should notice the contribution of the concentration of the electromagnetic field at the ribbon edge in THG. This effect would be much more important in the case of high-intensity illumination at higher frequencies where the difference between an infinite graphene sheet and a ribbon array (sensitivity to the inter-ribbon space) is more distinguishable for the illuminating wave than at lower frequencies (corresponding to longer wavelengths).

In all of the above simulations, to calculate the amount of enhancement in CE, we have set the intensity of the illuminating wave to be as low as possible (the TH level from the infinite graphene sheet is around −140 dB) such that the red shift in the resonant frequency of the system is negligible and the resonant frequency in the linear case provides a good approximation of the one in the nonlinear case. However, certainly by increasing the illumination intensity, to achieve higher levels of TH, one should be aware that the resonant frequency of the system undergoes a red shift and by working in the resonant frequency obtained from linear simulation, a reduction in THG improvement due to working out of the resonant region may be observed which is demonstrated in figure 16.

Zoom In Zoom Out Reset image size

All of the results reported in this study are calculated for the transmission mode (according to figures 1 and 7). However, according to our simulations, the results obtained by working in the reflection mode follow the same trend of variation with different parameters consistent with those obtained from the transmission mode and a negligible difference in the order of 2–3 dB in the TH level is observed which may be attributed to the transmission and reflection coefficients of the incident wave at the graphene boundary. Therefore, depending on the experimental facilities or convenience one of these two modes can be employed.

It is worth mentioning that the mechanisms of laser-induced damage to graphene, when illuminated by femtosecond (fs) pulsed lasers compared to CW lasers, are different. With fs pulsed lasers, both thermal and non-thermal effects have a contribution, whereas with CW lasers thermal effects play a dominant role. In terms of peak intensity, graphene can tolerate higher intensities for a shorter pulse at lower average intensity. In total, the reported damage threshold for fs pulsed lasers is ${I}_{th}\sim 2.7\;{\rm{TW}}\;{{\rm{cm}}}^{-2}$ which is much higher than that observed for CW lasers, $\sim 1\;{\rm{MW}}\;{{\rm{cm}}}^{-2}$ [50, 51, 67–69]. The intensities in our study are much lower than the threshold intensity of graphene damage. Since no loss is taken for the graphene bounding dielectric medium in our simulations, obtaining the same results in an experimental case may be achieved by increasing the intensity of the illuminating wave.

The intensities required for our predicted nonlinear responses can be provided by currently available or recently developed THz lasers [70] such as CO₂ lasers or molecular THz lasers (up to $10\;{\rm{MW}}\;{{\rm{cm}}}^{-2})$ [71], THz lasers with peak electric fields of $\sim \;4\;{\rm{MV}}\;{{\rm{m}}}^{-1}$ [72], single-cycle THz pulses with electric fields greater than $\sim 100\;{\rm{MV}}\;{{\rm{m}}}^{-1}$ [73] and coherent transition radiation sources from relativistic electron bunches with field strength in the order of $\mathrm{1-100}\;{\rm{MV}}\;{{\rm{cm}}}^{-1}\;$ [74].

In this study, we have demonstrated how electromagnetic field enhancement in a graphene microribbon array, under plasmon resonance conditions, can noticeably contribute to the improvement of THG strength, by more than five orders of magnitude based on our FDTD numerical simulations, compared to an infinite graphene sheet. By considering the electric tunability of the graphene first- and third-order conductivities besides the ability of graphene to support the SPP mode at THz frequencies, valuable potential applications for the THz gap region can be developed. We have shown that the resonant frequency of the ribbon array undergoes a blue shift with an increase in graphene Fermi level and period length and a decrease in ribbon width. We have revealed that for a graphene microribbon array of a higher Fermi level and quality (higher carrier mobility and longer relaxation time), a longer period and a narrower ribbon width, a more noticeable improvement in TH level is expected due to the formation of plasmon resonances of higher quality factors. We have explored how the modification of graphene conductivity under intense illumination and the corresponding red shift in the plasmon resonant frequency may lead to a reduction in THG CE in the case of excitation of the system at the resonant frequency obtained from linear simulation, due to working out of the resonance state.