## Abstract

We propose a hybrid structure where graphene is inserted to the interface of two one-dimensional photonic crystals (1D PCs). The two PCs are designed to have opposite topological properties, and at the interface, topological edge modes exist. The edge modes exist at both the fundamental frequency (FF) and the third harmonic (TH) frequency. This double resonant structure will enhance the nonlinear responses of graphene greatly, including Kerr nonlinearity and TH generation. We discuss these two kinds of nonlinearities both at terahertz (THz) and near-infrared (NIR) frequencies. The influence of Kerr nonlinearity on the resonant frequencies is considered, when we calculate the TH generation. At THz frequency, low-threshold bistability (about 8MW/cm^{2}) is obtained and the TH generation efficiency of 2.5% is achieved with incident intensity of 10MW/cm^{2}. At NIR frequency, the nonlinear conductivities of graphene are about 7 orders lower. Bistability is unlikely to happen with incident intensity below 1GW/cm^{2}. The TH generation efficiency is only about 5×10^{−6} with incident intensity of 25MW/cm^{2}. The proposed structure is more suitable to work as a low-threshold saturable absorber at NIR frequency. These results may be helpful both for a better understanding of graphene’s nonlinear responses in a double resonant structure and for potential applications in THz nonlinear devices and NIR nanophotonics.

© 2019 Optical Society of America under the terms of the OSA Open Access Publishing Agreement

## 1. Introduction

Nonlinear optical effects, which occur with strong light illumination, are indispensable in modern physics. Various phenomena and important applications arise with sufficiently high field strength ranging from frequency conversion, all-optical signal processing to ultrafast switching [1]. However, optical nonlinearities of common materials are inherently weak. To allow significant nonlinear responses to build up, usually bulky materials with enlarged interaction length and intense light are used. Moreover, for frequency conversion, precise and cumbersome phase matching techniques are always required to compensate for the momentum mismatch between pump and generated signals. The inherently weak nonlinearities hinder the efficient light manipulation at subwavelength scale and the applications in ultrafast compact optical devices [2–4]. Resonant structures including metasurfaces [2], metal plasmonics [3] and dielectric resonators [4] have been utilized to boost the nonlinearities at subwavelength scale. However, materials with inherent large nonlinearity are always pursued in modern photonics.

Graphene, a two-dimensional crystal consisting of a single layer of carbon atoms, has attracted great interests since it was discovered in 2004 [5]. Graphene is first known to have remarkable electronic properties due to its massless Dirac fermions [6]. Soon after, its unique linear optical responses are demonstrated both theoretically and experimentally [7,8], such as the universal linear absorption of ∼2.3% at visible wavelengths [9], the highly confined and low loss plasmonics at infrared and terahertz frequencies [10], and the highly tunable linear surface conductivity [11]. Meanwhile, it was theoretically predicted that the nonlinear electromagnetic response of graphene is inherently strong due to its Dirac cone band structure [12]. Since then, a lot of efforts have been paid to demonstrate its ultra-strong nonlinearities, including four-wave mixing [13–15], second/third-harmonic generation [16–19], Kerr nonlinearity and nonlinear optical absorption [20–22].

Though, graphene’s atomic thickness facilitates its applications in nanophotonics and no phase matching technique is needed if light propagates perpendicularly to the graphene, such short interaction length between graphene and light hinders its real applications. Covering graphene on dielectric waveguides can enlarge the interaction length between light and graphene [15,22], while phase matching is still required. Other resonant structures or mechanisms are also studied where no phase matching is required, including graphene plasmonics [23–25], photonic crystals [26,27], and simple Fabry-Perot (FP) cavity [28]. The enhanced nonlinear response originates from the giant local field enhancement when the incident light is at the resonant frequency. So far, the resonant enhancement of graphene’s nonlinearity has been demonstrated experimentally in a FP cavity at NIR wavelength [28] and graphene plasmonics at THz frequency [24], over two orders of enhancement was achieved. More delicate resonant structures were studied theoretically to show much larger nonlinear enhancement, including double resonant graphene plasmonics [23] and photonic crystals [26], where both the FF and TH fields are resonant with the structure. However, in most theoretical studies, the nonlinearity of graphene is modeled simply by a third order susceptibility *χ ^{(3)}* measured experimentally at a specific wavelengths, or a much simplified nonlinear surface conductivity. According to the recent experimental results [18,19,22] and theoretical predictions [29–31], the nonlinear responses of graphene are very complicated and are influenced by many conditions, such as temperature, Fermi level, electron relaxation time, nonlinear saturation, etc. In addition, in resonant structures, the large Kerr nonlinearity of graphene may shift away the resonant frequency and lower the local field enhancement at the desired frequency. All these effects should be carefully analyzed in theoretical studies to model graphene’s nonlinearities correctly.

Recently, topological edge modes in one-dimensional (1D) PCs have attracted much attention, due to its analogous to the 1D topological insulators in solid-state physics [32–34]. Apart from its significance in fundamental physics, i.e. the topological description of interface states, it also provides a systematic method to design 1D PCs with interface states appearing in the desired photonic band gaps [32]. When the edge mode is excited, large field enhancement can be expected at the interface. Especially, it has been demonstrated that the topological edge modes can be designed to be resonant with both the FF and TH frequency [33], while the traditional defect state in a 1D PC or equivalently described as a FP cavity state, can only be resonant with either the FF or TH frequency [27]. In addition, the topological properties of the bulk 1D PC guarantees the creation of edge modes which is robust for enhancing light-matter interaction.

In this work, we incorporate graphene into the interface of two 1D PCs with opposite topological properties. By carefully designing the two 1D PCs, topological edge modes will appear at both the FF and TH frequency, and very large enhancement of graphene’s nonlinearities can be achieved. To model the graphene’s nonlinearities precisely, nonlinear surface conductivities derived by a full quantum theory method are used [31]. However, this method is based on the perturbative theory of density matrix. One should always carefully estimate the saturation intensity before using these conductivities [30]. In addition, we also consider the influence of finite temperatures according to [29], and realistic relaxation rates are used. To calculate the linear 1D PCs with nonlinear graphene, a quite simple transfer-matrix method (TMM) is proposed. The previous nonlinear TMM treat graphene as a traditional material with *χ ^{(3)}* [35] and is complicated. A scattering-matrix method is also proposed in [25], which is also quite complicated, and treat graphene as a surface current boundary. In addition, these methods can only work for frequency conversion, the Kerr nonlinearity and two-photon absorption can not be included. The TMM in this work treat graphene as a nonlinear surface current boundary, which simplified the method greatly. The phenomenon of Kerr nonlinearity, two-photon absorption and TH generation can be solved precisely simultaneously. To guarantee the accuracy of our TMM, all the results in this work are also calculated by the finite element method (FEM, Comsol Multiphysics).

In this study, we discuss the performance of the proposed structure both at THz and NIR frequencies. From mid-infrared to THz frequencies, the nonlinear responses of highly doped graphene are mainly attribute to the intra-band transition and scale with frequency by 1/*ω*^{3}. Thus, the nonlinear responses are inherently strong at THz frequency. We demonstrate that the Kerr nonlinearity of graphene shifts the resonant frequency of the edge modes so much that low-threshold bistability is achieved. Meanwhile, the tuning effect of Femi energy on the linear response of graphene also influences the resonant frequency greatly. Combining these two effects, the resonant frequencies of FF and TH modes are matched. Maximum TH generation efficiency of more than 2% is achieved with incident intensities of several MW/cm^{2}. At NIR wavelengths, several resonant peaks of graphene’s nonlinear conductivities may appear at zero temperature [31], corresponding to resonant transitions induced by one, two and three photon absorption in graphene. However, at room temperature these resonant peaks are smeared out, and the nonlinear conductivities of graphene are 7 orders lower than that at THz frequency. Nevertheless, using the double resonant structure, maximum TH generation efficiency of 5×10^{−6} is achieved. We believe the results will not only have potential applications in nonlinear THz devices and NIR nanophotonics, but also deepen our understanding of graphene nonlinearities in double resonant structures.

## 2. Theoretical models and methods

Firstly we will discuss the strategy of designing the topological edge modes shortly, and one can find more details in Refs. [32–34]. As shown in Fig. 1(a), the red dashed rectangles mark the unit cells of PC1 and PC2. According to [33], if the topological edge modes are expected to exist at the odd photonic band gaps (in this work, the first and third band gap are chosen), we should first design PC1 with opposite sign of surface impedance at the first and third band gap. Following [33], we define $\alpha \Lambda = {n_A}{d_A} + {n_B}{d_B}$, where *n _{A}* and

*n*are the refractive indexes of the dielectric layers,

_{B}*d*and

_{A}*d*are the thicknesses, $\Lambda = {d_A} + {d_B}$ is the thickness of a unit cell and α is the ratio of optical path length of a unit cell to its thickness. α is chosen so that $2/\alpha = 1/{n_\alpha } + 1/{n_B}$, which is derived from ${n_\alpha }{d_\alpha } = {n_B}{d_B}$, and ${n_A} < {n_\alpha }$ should be fulfilled [33]. Thus, we set

_{B}*n*= 1.46,

_{B}*n*= 1.9, α = 1.665 and the thicknesses

_{A}*d*and

_{A}*d*can be derived by the above equations. For consistency, we use the same

_{B}*n*and

_{A}*n*at both THz and NIR frequencies. At THz frequency, layer B can be polymethylpentene (TPX), layer A can be SiO

_{B}_{2}, and the substrate is high-resistivity silicon with n

_{s}= 3.415. At NIR wavelengths, layer B can be SiO

_{2}, layer A can be ZnO and the substrate is silicon with n

_{s}= 3.49. Note that using other material parameters will not change our results qualitatively, only α and Λ should be changed for a different design. In Fig. 1(b), we show the band structure (black line) of PC1 defined as:

*i*= A or B). In this figure, the wave vector and frequency are all normalized by Λ. According to the method used in [32], we derive the Zak phase of every band of PC1, marked by red ‘0’ and ‘π’. The sign of surface impedance of the nth band gap can be derived as:

*l*is the number of crossing points under the

*n*th gap, $\theta _m^{\textrm{Zak}}$ is the Zak phase of the

*m*th band. For the case in Fig. 1(b), we mark every band gap by red rectangle if ${\zeta ^{(n )}} > 0$, and by blue rectangle if ${\zeta ^{(n )}} < 0$. For PC1, ${\zeta ^{(1 )}} > 0$ and ${\zeta ^{(3 )}} < 0$. To obtain topological edge modes in both the first and third band gaps, we should design PC2 to have opposite sign of the surface impedance. Fortunately, the strategy is quite simple [33], all we need is to exchange layer A and layer B in PC2, as shown in Fig. 1(a). In Fig. 1(c), we show the band structure (black lines) of PC2, which is exactly the same as that of PC1. However, due to the opposite formation of a unit cell, PC2 have opposite Zak phases for every band comparing to PC1. The sign of surface impedance of PC2 at the first and third band gap is just opposite to that of PC1. We can expect topological edge modes to occur in these two gaps simultaneously. In Fig. 1(d), we show the reflectance spectrum of the PC1-PC2 structure without graphene, we set

*N*=

_{1 }*N*

_{2}_{ }= 7. Just as we expected, in the first and third photonic band gaps, resonant edge modes exist appearing as reflectance dips. Then if we incorporate graphene to the interface of PC1 and PC2, marked by black dashed line in Fig. 1(a), resonant enhanced nonlinear responses can be expected.

Next, we introduce our TMM for completeness. This TMM is based on the well known linear TMM of [36] where Gaussian units are used. We rewrite the matrix using SI units as follows:

*E*and

*H*,

*ɛ*and –

*μ*,

*Z*and 1/

_{0}*Z*. In Eq. (3), a dielectric layer with thickness

_{0}*d*and dielectric constant

*ɛ*is considered.

*k*is the

_{z}*z*component of the wave vector in this dielectric,

*k*=

_{0 }*ω*/c is the vacuum wave vector,

*Z*is the vacuum impedance,

_{0}*H*and

_{y}*E*are the transverse magnetic and electric fields. This kind of TMM is based on the boundary condition of continuous transverse fields, i.e. $H_y^{{0^ + }} = H_y^{{0^ - }}$ and $E_x^{{0^ + }} = E_x^{{0^ - }}$ if the interface locates at

_{x}*z*= 0. However if graphene exist at

*z*= 0, the boundary condition should be modified by the linear surface conductivity of graphene

*σ*, i.e. $H_y^{{0^ - }} - H_y^{{0^ + }} = {\sigma _l}E_x^{0 + }$ and $E_x^{{0^ + }} = E_x^{{0^ - }}$. We can define the transfer matrix of graphene as:

_{l}*H*and

_{y}*E*. One can derive them easily using Maxwell equations, whenever it is necessary. Then the total response of the layered structure in Fig. 1(a) can be described as:

_{x}In this work, TM light is considered to normally incident on the left surface of PC1, as shown in Fig. 1(a). Our method is also valid for oblique incidence. We assume the transmitted fields to be ${E_{xs}} = T$, $ {H_{ys}} = ({\omega {\varepsilon_0}{\varepsilon_s}/{k_{zs}}} )T$ according to the Maxwell equation ${E_x} = 1/({i\omega {\varepsilon_0}\varepsilon } )\partial {H_y}/\partial z$. We set the location of graphene to be at *z* = 0. Then the field on the right of graphene is:

*n*is the refractive index.

After we get the electric field at the graphene $E_x^{{0^ - }} = E_x^{{0^ + }}$, we can calculate the TH fields as follows. The TH generation associates with graphene’s nonlinear conductivity $\sigma _{xxxx}^{(3 )}({\omega ,\omega ,\omega } )\; $[31], rewrite as $\sigma _{TH}^{(3 )}$. All the variables at the TH frequency (3*ω*) will be marked by tilde. For the TH fields, no incident light exists, we can define the reflected and transmitted TH fields as ${\tilde{H}_{y0}} = \tilde{R}$, ${\tilde{E}_{x0}} ={-} {\tilde{k}_{z0}}/({\tilde{\omega }{\varepsilon_0}} )\tilde{R}$ and ${\tilde{H}_{ys}} = \tilde{T}$, ${\tilde{E}_{xs}} = {\tilde{k}_{zs}}/({\tilde{\omega }{\varepsilon_0}{{\tilde{\varepsilon }}_s}} )\tilde{T}$. The following equations are fulfilled defining $\tilde{N} = {({{{\tilde{M}}_{{d_B}/2}}{{\tilde{M}}_{{d_A}}}{{\tilde{M}}_{{d_B}/2}}} )^{ - {N_1}}}$, $\tilde{M} = {({{{\tilde{M}}_{{d_A}/2}}{{\tilde{M}}_{{d_B}}}{{\tilde{M}}_{{d_A}/2}}} )^{{N_2}}}$:

*ω*and finally the fields $\tilde{R}$ and $\tilde{T}$ can be calculated by Eq. (10). To conclude, our TMM can calculate both the Kerr nonlinearity and TH generation of nonlinear graphene embed in linear layered structures, and the appearance is much simpler than that of [35]. Our method is also valid for other 2D materials like the edge state of topological insulator, MoS

_{2}, etc.

All the TMM results in this work are confirmed by the FEM, and we discuss the FEM shortly. In the FEM, two frequency domain calculations are implemented corresponding to the FF and TH frequency. The nonlinear graphene is modeled as surface current boundaries at both the FF and TH frequency. At the FF, the surface current is ${J_x} = {\sigma _l}{E_x} + 3\sigma _K^{(3 )}{|{{E_x}} |^2}{E_x}$. At the TH frequency, no incident light exists and the surface current is ${\tilde{J}_x} = {\tilde{\sigma }_l}{\tilde{E}_x} + \sigma _{TH}^{(3 )}E_x^3$, note that the variables with tilde are evaluated at 3*ω* and $\sigma _{TH}^{(3 )}E_x^3$ without tilde is the source term evaluated at the FF. It is vital important that one should use small enough mesh size comparing both to the fundamental and TH wavelengths, for example 1/100 of fundamental wavelength. Using a mesh size 1/30 of fundamental wavelength can achieve correct results only for FF. However, the TH generation may be incorrect because the mesh size is only 1/10 of the TH wavelength.

At the end of this section, we discuss the linear and nonlinear conductivities of graphene. For the linear conductivity, we use the local-RPA model with finite temperature *T* [11,38]:

*E*is the Fermi energy, $E_\textrm{F}^T = {E_\textrm{F}} + 2{k_\textrm{B}}T\textrm{log}({1 + {\textrm{e}^{ - {E_\textrm{F}}/{k_\textrm{B}}T}}} )$, and ${f_\textrm{E}} = 1/[{1 + {\textrm{e}^{({E - {E_\textrm{F}}} )/{k_\textrm{B}}T}}} ]$ is the Fermi-Dirac distribution function.

_{F}*τ*is the phenomenological relaxation time, ${v_F} \approx 1 \times {10^6}\textrm{m}/\textrm{s}$ is the Fermi velocity. From mid-infrared to THz frequency, $\tau = \mu {E_\textrm{F}}/ev_\textrm{F}^2$ where

*μ*is the impurity-limited dc mobility, indicating that the linear loss of graphene is mainly attributed to impurities. The reported maximum

*τ*estimated from

*μ*is about 0.5ps for graphene on a clean SiO

_{2}substrate [39], while for low quality graphene grown by CVD method

*τ*is only about 0.1ps [39]. Recently, very high relaxation times over 1ps have been demonstrated using graphene/hBN hybrid structures [40]. For photon energy higher than 0.2eV ($\lambda < 6.2$µm), the relaxation-time reduces significantly due to the electron-phonon coupling, usually around 0.1ps [41]. The nonlinear conductivities of graphene are calculated according to the full quantum theory [31]. The expressions are very complicated, so we do not list them in this work. We guarantee the correctness of the conductivities by repeating all the figures in [31]. We will discuss these nonlinear conductivities whenever it is necessary. However, the results in [31] were calculated at zero temperature. We should extend them to finite temperature by [29]:

*T*= 300K is assumed whenever one-photon absorption is prohibited.

## 3. Results and discussions

The reflectance in Fig. 1(d) is calculated with normalized frequency $\omega \Lambda /({2\pi c} )$ or equivalently $\Lambda /\lambda $, where λ is the wavelength. We can achieve the desired resonant wavelength of edge modes according to Fig. 1(d). Firstly, we will show the performance of our structure in THz range. The normalized frequency of the first edge mode is about $ \Lambda /\lambda \approx 0.3$, if $\lambda = 100$µm (3THz) we should set $\Lambda = 0.3\lambda = 30$µm. The thicknesses of A and B layers can be calculated according to $\alpha \Lambda = {n_A}{d_A} + {n_B}{d_B}$ and $\Lambda = {d_A} + {d_B}$. The linear responses of the edge modes without and with graphene are shown in Fig. 2.

In Figs. 2(a) and 2(b), the reflectance around 3THz and 9THz are shown respectively. The blue solid line is the result without graphene, and a large mismatch exists between the resonant frequencies of two edge modes, for convenience these two resonant frequencies are denoted as *f _{FF}* and

*f*. The mismatch between

_{TH}*f*and

_{FF}*f*/3 is due to the asymmetry of the dielectric environment on the two sides of the interface. However, as we will show below this mismatch is beneficial for the hybrid structure with graphene at THz frequency. In Figs. 2(a) and 2(b), the reflectance where graphene has different Fermi levels is also shown and

_{TH}*τ*= 0.5ps is fixed. At THz frequency, the interaction between graphene and light is much stronger than that of infrared or visible frequency. The influence of graphene’s Fermi level on the resonant frequency of various structures has been demonstrated extensively [38,39]. Here, with increased

*E*the resonant frequencies of both FF and TH blue shift. At THz frequency, the linear conductivity of highly doped graphene can be simplified to the traditional Drude type expression [38], ${\sigma _l}(\omega )= i{e^2}{E_F}/({\pi {\hbar^2}\omega + \textrm{i}\pi {\hbar^2}{\tau^{ - 1}}} )$, and finite temperature has no influence on it.

_{F}*σ*is proportional to 1/

_{l}*ω*, so the linear response of graphene is expected to be much smaller around

*f*. As expected, the blue shift of

_{TH}*f*in Fig. 2(b) is much smaller than

_{TH}*f*. For

_{FF}*E*larger than 0.4eV, 3

_{F}*f*becomes larger than

_{FF}*f*. This condition is necessary if the Kerr nonlinearity of graphene is considered, which will red shift

_{TH}*f*and match

_{FF}*f*and

_{FF}*f*/3. In Figs. 2(c) and 2(d), we show the normalized electric field intensity at

_{TH}*f*and

_{FF}*f*, where

_{TH}*E*

_{F}_{ }= 0.5eV is fixed. In Fig. 2(c), the results with

*f*= 3.05THz,

_{FF}*τ*= 0.5ps or 0.25ps are shown. As expected, field intensity enhancement is achieved near the interface. When

*τ*is decreased by a half, the field intensity enhancement also reduces by about a half, showing the significant influence of graphene’s impurities on the resonant enhanced fields. In Fig. 2(d), the results with

*f*= 9.07THz,

_{TH}*τ*= 0.5ps or 0.25ps are shown. The field intensity enhancement is much larger than that in Fig. 2(c), and reducing the quality of graphene influences less on the TH resonant field enhancement. We can expect that if

*f*is matched with

_{FF}*f*/3, field intensity enhancement will occur simultaneously at the FF and TH frequency. According to [42,43], the TH generation can be estimated qualitatively by these linear fields:

_{TH}*E*= 0.5 eV and

_{F}*τ*= 0.5ps, the linear mode field ${E_x}(\omega )$ is enhance by a factor of about $\sqrt 6 $, the linear mode field ${\tilde{E}_x}({3\omega } )$ is enhance by a factor of about $\sqrt {30} $. Thus, the TH intensity will be enhanced by a factor of ${6^3} \times 30 = 6480$. Note that Eq. (13) is only for estimation, TMM or FEM is always used to calculate the TH quantitatively.

In Fig. 3, the influence of graphene’s Kerr nonlinearity on the resonance at the FF is shown. For the frequency interval considered here, $\sigma _K^{(3 )}$ is nearly constant. For convenience, we transform $\sigma _K^{(3 )}$ to the bulk nonlinear third order susceptibility $\chi _K^{(3 )}$ by $\chi _K^{(3 )} = \chi _\textrm{K}^{(3 );xxxx}({\omega ,\omega , - \omega } )= 3i\sigma _K^{(3 )}/({{\varepsilon_0}\omega {d_{gr}}} )$ according to Eq. (67) in [29], where *d _{gr }*= 0.34 nm is the thickness of graphene. Note that the factor of 3 in the numerator is due to the frequency permutation symmetry, i.e $\sigma _{xxxx}^{(3 )}({\omega ,\omega , - \omega } )= \sigma _{xxxx}^{(3 )}({\omega , - \omega ,\omega } )= \sigma _{xxxx}^{(3 )}({ - \omega ,\omega ,\omega } )= \sigma _K^{(3 )}$. At

*f*= 3THz,

*E*= 0.5 eV and

_{F}*τ*= 0.5ps we get $\chi _K^{(3 )} = ({6.84 \times {{10}^{ - 10}} - 3.62 \times {{10}^{ - 11}}i} ){\textrm{m}^\textrm{2}}/{\textrm{V}^\textrm{2}}$, the influence of temperature on $\chi _K^{(3 )}$ is negligible if

*E*is much larger than ${k_B}T \approx 26\textrm{meV}$. The positive real part of $\chi _K^{(3 )}$ indicates self-focusing type of Kerr nonlinearity, and the negative imaginary part indicates perturbative saturable absorption. Strictly speaking, the saturable absorption is a nonperturbative nonlinear phenomenon [44]. Nevertheless, we can estimate the real saturable field using the perturbative $\sigma _K^{(3 )}$ though ${E_{sat}} = \sqrt { - \textrm{Re}({{\sigma_l}} )/3\textrm{Re}({\sigma_K^{(3 )}} )} $ [30]. The perturbative $\sigma _K^{(3 )}$ can be used if the enhanced electric field is smaller than

_{F}*E*. At

_{sat}*f*= 3THz,

*E*= 0.5eV,

_{F}*τ*= 0.5ps, the estimated saturation intensity $ {I_{sat}} \approx 84\textrm{MW}/\textrm{c}{\textrm{m}^\textrm{2}}$. The influence on Kerr nonlinearity and saturation intensity is negligible even if we decrease

*τ*to 0.1ps. In Fig. 3(a), we show the influence of incident light intensity

*I*on the edge mode of FF. Obvious red shifts of

_{0}*f*can be observed. At around

_{FF}*I*= 8MW/cm

_{0}^{2}, the reflectance at 3THz drops suddenly, which is the indication of optical bistability [45]. In Fig. 3(b), the reflectance at 3THz and

*E*= 0.5eV with varying incident intensity and different

_{F}*τ*is shown. For

*τ*= 0.5ps low-threshold bistability appears as blue solid line, the blue dashed line is calculated by FEM with increasing or decreasing

*I*. The upper threshold where reflectance drops suddenly with increasing

_{0}*I*is about 8MW/cm

_{0}^{2}, the lower threshold where reflectance increase suddenly with decreasing

*I*is about 5MW/cm

_{0}^{2}. When

*τ*is decreased to 0.35ps, these thresholds both are increased. If

*τ*is around 0.25ps, the bistability disappears. Figure 3(b) shows the significant influence of the graphene’s quality on its Kerr nonlinear applications.

In Fig. 4, the TH generation is discussed. According to [31], $\sigma _K^{(3 )}$ at THz frequency is mainly attribute to the intra-band transition, and scales with 1/ω^{3}. Thus, the effect of graphene’s Kerr nonlinearity on *f _{TH}* can be ignored. In Fig. 4(a), the reflectance spectra with different

*I*are shown, with

_{0}*E*= 0.5eV and

_{F}*τ*= 0.5ps. The vertical red dash-dotted line marked the 1/3 of

*f*= 9.07THz, which is constant with

_{TH}*I*considered in this figure. Similar to $\sigma _K^{(3 )}$, $\sigma _{TH}^{(3 )}$ is nearly constant in this frequency interval. And bulk nonlinear third order susceptibility is $\chi _{TH}^{(3 )} = \chi _{\textrm{TH}}^{(3 );xxxx}({\omega ,\omega ,\omega } )= i\sigma _{TH}^{(3 )}/({3{\varepsilon_0}\omega {d_{gr}}} )$ [29]. At

_{0}*f*= 3.023THz, $|{\chi_{TH}^{(3 )}} |= 2.55 \times {10^{ - 11}}{\textrm{m}^\textrm{2}}/\textrm{V}$ which is about 7 orders larger than that at NIR frequency [18,19]. In Fig. 4(a), with increased

*I*,

_{0}*f*red shifts due to graphene’s Kerr nonlinearity. At around

_{FF}*I*= 4MW/cm

_{0}^{2},

*f*approximates

_{FF}*f*/3, and double resonant enhancement of TH generation can be expected. In Fig. 4(b), the conversion efficient of the reflected TH field defined by

_{TH}*η*=

*I*/

_{3ω}*I*is shown. As expected, when

_{0}*f*and

_{FF}*f*/3 are matched, large

_{TH}*η*is achieved. To be more specific, for

*I*= 2MW/cm

_{0}^{2}

*f*and

_{FF}*f*/3 are not matched,

_{TH}*η*= 0.05% at

*f*, and

_{FF}*η*is maximized to about 0.1% when frequency is between

*f*and

_{FF}*f*/3. For frequency far away from

_{TH}*f*or

_{FF}*f*/3,

_{TH}*η*is only about 10

^{−7}%. For non-resonant frequencies, TH field is proportional to $ I_0^{3/2}$, and

*η*will be enhanced proportionally to $ I_0^2$, we call this kind of enhancement as intrinsic. The enhancement of

*η*due to the double resonant is called extrinsic enhancement. When

*I*is increased to 4MW/cm

_{0}^{2}, the intrinsic enhancement is about 4 times. However, maximized

*η*= 1.4% appears around

*f*and

_{FF}*f*/3, which is 14 times larger. When

_{TH}*I*= 6MW/cm

_{0}^{2}, though

*f*seems to be away from

_{FF}*f*/3, combining the intrinsic and extrinsic enhancement, at around

_{TH}*f*and

_{FF}*f*/3, maximum

_{TH}*η*= 2.4% is still achieved. When

*I*is further increased to 10MW/cm

_{0}^{2}, the efficiency reduction due to the mismatch of

*f*and

_{FF}*f*/3 and the intrinsic enhancement are balanced, maximum

_{TH}*η*= 2.5% is achieved. In Fig. 4(c), the process discussed above is clearly shown by showing

*η*with increasing

*I*, and cases with different

_{0}*τ*are shown. We fix

*f*=

*f*/3 = 3.023THz, and

_{TH}*E*= 0.5eV. The blue solid line show the result with

_{F}*τ*= 0.5ps. Clearly,

*η*shows a saturation behavior, which is not due to the nonlinear saturation of graphene (all the intensities here are well below the saturation field even considering the local field enhancement at

*f*). Below

_{FF}*I*= 6MW/cm

_{0 }^{2}, as

*I*increases

_{0}*η*increases quickly due to the combined intrinsic and extrinsic enhancement. From 6 to 10MW/cm

^{2},

*η*is nearly constant due to the balance between intrinsic enhancement and the reduction of extrinsic enhancement. Above 10MW/cm

^{2},

*η*reduces very slowly due to the large mismatch between

*f*and

_{FF}*f*/3. When

_{TH}*τ*is reduced to 0.25ps,

*η*is reduced by more than 2 times, this is mainly attributed to the reduction of local field enhancement at

*f*as shown in Fig. 2(c). When the quality of graphene is very low and

_{FF}*τ*= 0.1ps,

*η*is reduced significantly as shown by the black dash-dotted line, the value of

*η*is multiplied by 10 here. As can be seen,

*η*is increased with increasing

*I*, and no saturation effect appears. This indicates that the intrinsic enhancement of TH generation plays the leading role when the graphene quality is low.

_{0}*η*= 0.09% when

*I*= 15MW/cm

_{0}^{2}, however this is also 2 orders larger than that at non-resonant frequencies (not shown here).

To conclude, when graphene’s quality is high, the limitation of even larger *η* is the mismatch of *f _{FF}* and

*f*/3. If graphene’s quality is low, larger

_{TH}*I*will results in larger

_{0}*η*, the limitation then will be the saturation effects of graphene’s nonlinear conductivity $\sigma _{TH}^{(3 )}$. The effects of the nonlinear saturation on $\sigma _{TH}^{(3 )}$ is complicated and non-perturbative method should be used [30,44]. This is out of the scope of this work, and

*I*here is chosen to be smaller than the estimated saturation intensity.

_{0}The above results are vailed for finite temperatures below *T *= 300 K, where the influences of *T* on both the linear and nonlinear responses of graphene are negligible. However, recent experiments clarify strong effect of hot electrons on the graphene’s THz response [46,47]. Thus, it is necessary to discuss the influence of higher electron temperatures on the graphene’s linear and nonlinear THz responses. According to [18,48–50], higher incident intensities and higher absorptance will result in higher electron temperature *T _{e}* of the graphene. We can estimate

*T*as follows [18]:

_{e}*T*= 300 K is the lattice temperature, $\textrm{F}\; ({\textrm{J} \cdot {\textrm{m}^{-2}}} )$ and $\Delta t$ are the fluence and duration of the incident light pulse, so that the incident intensity ${I_0} = \textrm{F}/\Delta t$. When the incident pulse is off,

_{0}*T*relaxes towards

_{e}*T*on a time scale ${\tau _r} = 100\textrm{fs}$.$\; \alpha $ is the absorptance of the structure and ${c_\textrm{v}}$ is the electronic heat capacity of the photoexcited graphene. According to [18], ${c_\textrm{v}}({\mu ,{T_\textrm{e}}} )= {\pi ^2}\nu ({{E_\textrm{F}}} )k_\textrm{B}^2{T_\textrm{e}}/3$, if ${k_B}{T_\textrm{e}} \ll |\mu |$ is fulfilled. Here $\nu ({{E_F}} )= 2|{{E_F}} |/[{\pi {{({\hbar {v_F}} )}^2}} ]$ is the density of electronic states per unit area at the Femi energy. Note that the Femi energy ${E_F}$ is the chemical potential $\mu $ at ${T_e} = 0K$. For finite ${T_e} \lesssim {T_F}$, $\mu = {E_F}[{1 - {\pi^2}T_\textrm{e}^2/({6T_\textrm{F}^2} )} ]$, where ${T_F} = {E_F}/{k_B}$. However, the difference between ${E_F}$ and $\mu $ below 300 K is negligible, we can safely replace $\mu $ by ${E_F}$ at low temperatures. When ${T_e}$ is larger, we will always use $\mu $ instead of ${E_F}$.

_{0}In Fig. 5, we show the influence of large ${T_e}$ on the results obtained before. Firstly, in Fig. 5(a), we show the estimated ${T_e}$ with increasing ${I_0}$ (blue solid line, left axis). The ratio of $\mu /{E_F}$ is also shown as a reference (red dashed line, right axis). The absorptance $ \alpha $ is about 0.33, which is obtained at *f _{FF}* in Fig. 2(a) with

*E*= 0.5 eV,

_{F}*τ*= 0.5ps. As we can see, with increased ${I_0}$ from 0.1 to 15 MW/cm

^{2},

*T*is increased from 300 K to 1500 K. Meanwhile, $\mu $ do not deviate ${E_F}$ much and the influence of increased $ {T_e}$ on the linear conductivity of graphene is small (not shown here). For example, at

_{e}*T*= 1500 K, the absolute value of graphene’s linear conductivity is reduced by about 10%, which is proportional to the reduction of $\mu $. In Fig. 5(b), we show the influence of increased $ {T_e}$ on the nonlinear susceptibilities of graphene. The $|{\chi_K^{(3 )}} |$ (blue solid line, left axis) and $|{\chi_{TH}^{(3 )}} |$ (red dashed line, right axis) are shown respectively. We can see that both $|{\chi_K^{(3 )}} |$ and $|{\chi_{TH}^{(3 )}} |$ are increased a little, about 30% at

_{e}*T*= 1500 K. In Fig. 5(c), we show the influence of

_{e}*T*on

_{e}*η*, the other parameters remain the same as that in Fig. 4(c). As can be seen, the influence of higher

*T*is not significant. When

_{e}*I*< 5MW/cm

_{0 }^{2},

*η*seems to be larger with higher

*T*, however, according to Fig. 5(a),

_{e}*T*may not be larger than 1000 K with

_{e}*I*< 5MW/cm

_{0 }^{2}. To summarize, the influence of hot electrons on the nonlinear responses of our structure in the THz range is not significant, and the results we obtained with low temperature are not violated qualitatively.

We note that in [46], a very large $|{\chi_{TH}^{(3 )}} |$ is predicted, which is ∼${10^{ - 9}}{\textrm{m}^\textrm{2}}/\textrm{V}$, considering the large thermal responses of graphene at THz frequencies. This is two orders larger than that predicted by the perturbative theory. The conversion efficiency is about $4 \times {10^{ - 6}}$ for the graphene on a SiO2 substrate, the incident intensities are about 1-20MW/cm^{2} [46]. In our structure, *η* is still 2∼3 orders larger even with 2 orders lower $|{\chi_{TH}^{(3 )}} |$, which suggests giant enhancement factor.

The above results are obtained at THz frequency. When frequency increases, the linear conductivity of graphene scales with 1/*ω* before inter-band transition plays a role. Thus, the linear response of graphene is expected to be much weaker at NIR frequency than that at THz frequency. And the shift of resonant frequency is much smaller when *E _{F}* is varied at NIR frequency [38,39]. According to Figs. 2(a) and 2(b), without graphene

*f*and

_{FF}*f*/3 of the structure in Fig. 1(a) have a large mismatch at THz frequency. Note that we can achieve totally the same linear response at NIR frequency by only scaling down Λ according to Λ/

_{TH}*λ*= 0.3. The large mismatch can be eliminated by the effects of

*E*and

_{F}*I*at THz frequency. However, as we will see, the influence of

_{0}*E*and

_{F}*I*on the

_{0}*f*and

_{FF}*f*is much smaller at NIR frequency. We need to reduce this mismatch by modifying the structure in Fig. 1(a). After our investigation, the structure in Fig. 6(a) fulfilled the demand. Comparing the structure in Fig. 6(a) to that in Fig. 1(a), the unit cell of PC1’ is shifted toward positive

_{TH}*z*by

*d*/2 and the unit cell of PC2’ is shifted toward positive

_{B}*z*by

*d*/2. Around the new interface between PC1’ and PC2’, the dielectric environment is symmetric. This kind of interface state has also been used as topological edge modes [34], because PC1 (PC2) and PC1’ (PC2) is equivalent indeed if we ignore how these PCs are truncated at the interface of

_{A}*z*= 0. In this new structure,

*f*and

_{FF}*f*/3 are matched before graphene is inserted.

_{TH}The reason why we did not use this symmetric structure at THz frequency is that, at THz frequency, the blue shift of *f _{FF}* induced by

*E*is significant. Using

_{F}*E*= 0.05eV in the PC1’-PC2’ structure at THz frequency will cause similar mismatch comparing to

_{F }*E*

_{F}_{ }= 0.5eV in PC1-PC2 structure. If we use

*E*= 0.05eV and the PC1’-PC2’ structure at THz frequency, the nonlinear saturation effect will be significant and the influence of finite temperature is also vital. The saturation intensity will reduce to 30kW/cm

_{F }^{2}due to the combined effect of low Fermi energy and finite temperature.

*I*of several kW/cm

_{0}^{2}should be used for the calculations, and the TH generation efficiency is much lower than that in Fig. 4. We note in [31], very large TH generation efficiency is predicted at 3THz even for suspended graphene with very low

*E*(about 0.012eV according to ${E_F} = \hbar {v_F}\sqrt {\pi {n_g}} $) and

_{F}*I*= 1MW/cm

_{0}^{2}. However, for such low

*E*, the saturation intensity will be much lower and the perturbative nonlinear conductivities will not be valid for incident light intensity of 1MW/cm

_{F}^{2}.

For the PC1’-PC2’ structure shown in Fig. 6(a), the *n _{A}*,

*n*and

_{B}*α*remain unchanged. To make the resonant wavelength of the fundamental edge mode (

*λ*) at 1.55µm, $\Lambda /{\lambda _{FF}} = 0.3028$ is set. To further enhance the local field enhancement at the interface, we set

_{FF}*N*= 7 and

_{1 }*N*= 14. In Fig. 6(b), we show the linear reflectance spectra at NIR wavelength. At NIR wavelength, the relaxation time

_{2}*τ*is mainly limited by the coupling of electrons and phonons [41],

*τ*= 0.1ps is fixed for the results below. As shown in Fig. 6(b), even for

*E*= 0.5eV,

_{F}*λ*is not shifted obviously. For

_{FF}*E*= 0.8eV,

_{F }*λ*is blue shifted by only about 2nm. Meanwhile, the resonant wavelength of the TH edge mode (

_{FF}*λ*) is not shifted by the

_{TH}*E*and we mark the wavelength of 3

_{F}*λ*by the vertical red dash-dotted line. In Figs. 6(c) and 6(d), the normalized electric fields are shown at

_{TH}*λ*and

_{FF}*λ*respectively. Comparing to Fig. 2, the field enhancement is much larger.

_{TH}In Fig. 7(a), we show the nonlinear susceptibility $ \chi _K^{(3 )}$ at 1.55µm (ℏ*ω = *0.8eV) with varying *E _{F}* and

*τ*= 0.1ps. If 2

*E*≤ 0.8eV, one-photon absorption occurs, the saturation intensity of graphene due to the inter-band transition is very low, about only 10MW/cm

_{F}^{2}[44]. On the contrary, the saturation intensity of graphene due to intra-band transition is very high, about 100GW/cm

^{2}is obtained at

*E*= 0.8eV [44]. In addition, according to [31], at

_{F}*E*= 0.8eV two-photon absorption occurs, a resonant peak of $\textrm{Re}({\sigma_K^{(3 )}} )$ exists. Thus, in Fig. 7(a) $\chi _K^{(3 )}$ with

_{F}*E*larger than 0.5eV is shown, the results with zero temperature and

_{F}*T*= 300K are shown. As expected, with

*T*= 0 a resonant peak exists at

*E*= 0.8eV. However, this peak is smeared out by the finite temperature effect. Nevertheless, several times larger $\chi _K^{(3 )} = 7.05 \times {10^{ - 17}} + 4.71 \times {10^{ - 17}}{\textrm{m}^\textrm{2}}/{\textrm{V}^\textrm{2}}$ can still be obtained at

_{F }*E*= 0.8eV. Note that at

_{F }*E*= 0.8eV, $\textrm{Im}({\chi_K^{(3 )}} )$ is positive indicating two-photon absorption process. Comparing to $ \chi _K^{(3 )}$ at THz range, here $\chi _K^{(3 )}$ is 7 orders smaller. Even considering much higher field intensity enhancement in Fig. 6(c), which is 2 orders larger, the Kerr nonlinearity response at 1.55µm can be expected to be about 5 orders lower than that in Fig. 3. In addition, the relaxation time

_{F}*τ*is low even for high quality graphene at 1.55µm due to the phonon effect. Kerr nonlinear applications such as bistability are not likely to happen at NIR wavelength. On the contrary, for low doping levels of graphene (2

*E*≤ 0.8eV), one-photon absorption is allowed at NIR wavelength, graphene has been demonstrated to be a fast saturable absorber with low saturation intensity in a wide frequency range [7,8,44]. In Fig. 7(b), the reflectance spectra with different

_{F}*I*are shown, and

_{0}*E*= 0.8eV. Considering the saturation intensity is about 100GW/cm

_{F}^{2}and the field enhancement at the interface is about 100 times, the maximum

*I*is chosen to be 1GW/cm

_{0}^{2}. Even for

*I*= 1GW/cm

_{0}^{2}, the red shift of

*λ*is not significant. Nevertheless, the small red shift due to the Kerr nonlinearity may compensate the small blue shift due to large

_{FF}*E*.

_{F}In Fig. 8, the TH generation at NIR wavelength is shown. We show the value of $|{\chi_{TH}^{(3 )}} |$ with varying *E _{F}* at 1.55µm in Fig. 8(a). Two resonant peaks can be seen for $|{\chi_{TH}^{(3 )}} |$ with zero temperature, corresponding to two and three photon absorption. However the two peaks are smeared out by finite temperature effect. With

*T*= 300, $|{\chi_{TH}^{(3 )}} |$ is around 7×10

^{−18}m

^{2}/V

^{2}for

*E*= 0.8eV which is on the same order with recent experiments [18,19]. In Fig. 8(b), the TH generation efficient

_{F }*η*is shown with varying wavelengths and different

*I*. The maximum

_{0}*η*= 1.17% is obtained with

*I*= 1GW/cm

_{0 }^{2}at 1.55µm. In Fig. 8(c),

*η*with varying

*I*at 1.55µm is shown.

_{0}*η*increases nearly linearly with increasing

*I*, according to Fig. 7(b) and Fig. 8(b), both the intrinsic and extrinsic enhancement play vital roles. The influence of reducing

_{0}*τ*to 0.05ps is negligible, the maximum

*η*= 1.07% is obtained with

*I*= 1GW/cm

_{0 }^{2}. For intensity larger than 1GW/cm

^{2}, the nonlinearity is probable influenced by saturation effect considering the large local intensity enhancement of about 100 times.

Similar to the THz case, the above results for NIR wavelengths are valid only for *T* = 300 K. The intense incident light used above can result in large electron temperature *T _{e}* and its effects on the graphene’s responses are shown in Fig. 9. In Fig. 9(a), we show the estimated ${T_e}$ with increasing ${I_0}$ (blue solid line, left axis), according to Eq. (14). The ratio of $\mu /{E_F}$ is shown as a reference (red dashed line, right axis). The absorptance $ \alpha $ is about 0.03, which is obtained at

*λ*with

_{FF}*E*= 0.8 eV,

_{F}*τ*= 0.1ps. As we can see, with ${I_0}$ increased from 1 to 500 MW/cm

^{2},

*T*is increased from 300 K to 1900K. Meanwhile, $\mu $ do not deviate ${E_F}$ much. The influence of increased $ {T_e}$ on the linear conductivity of graphene is a little larger than that of $\mu /{E_F}$. At

_{e}*T*= 1900K, the absolute value of graphene’s linear conductivity is reduced by about 17%, and $\mu $ is reduced by only about 7%. In Fig. 9(b), we show the influence of increased $ {T_e}$ on the nonlinear susceptibilities of graphene. The imaginary part of $\chi _K^{(3 )}$ (blue solid line, left axis) and $|{\chi_{TH}^{(3 )}} |$ (red dashed line, right axis) are shown respectively. The influence of larger ${T_e}$ on $\textrm{Im}({\chi_K^{(3 )}} )$ is significant, which is entirely different from the case at THz frequency. For

_{e}*T*< 500 K, $\textrm{Im}({\chi_K^{(3 )}} )> 0$ indicates the two-photon absorption, and the saturation intensity is very large [44]. However, for

_{e}*T*> 500 K, $\textrm{Im}({\chi_K^{(3 )}} )< 0$ which indicates low saturation intensities. In the inset of Fig. 9(b), we show the enlarged curve around

_{e}*T*= 500 K. The $- \textrm{Im}({\chi_K^{(3 )}} )$ is increased by about two orders from 550 K to 1000 K. We can estimate the saturation intensity as before. At

_{e}*T*= 550 K, the saturation intensity is about 2.5GW/cm

_{e}^{2}. As

*T*increased to 700 K, the saturation intensity is about 250MW/cm

_{e}^{2}. For

*T*= 1000 K, the saturation intensity is only about 50MW/cm

_{e}^{2}. Considering that the intensity enhancement in Fig. 6(c) is ∼80, we should limit

*I*< 30MW/cm

_{0}^{2}for

*T*= 550 K. Meanwhile, according to Fig. 9(a), for

_{e}*I*< 25MW/cm

_{0}^{2},

*T*is smaller than 550 K. Thus, we limit

_{e}*I*below 25MW/cm

_{0}^{2}. If larger

*I*is used as before,

_{0}*T*will be larger, the nonlinear responses will be saturate and the perturbative theory may be failed. In Fig. 9(b), we can see that $|{\chi_{TH}^{(3 )}} |$ decreases with increased

_{e}*T*, just as demonstrated by the experiments [48] for highly doped graphene. In Fig. 9(c) we show the influence of

_{e}*T*on

_{e}*η*, the other parameters remain the same as that in Fig. 8(c). For larger

*T*, $|{\chi_{TH}^{(3 )}} |$ is decreased and

_{e}*η*is decreased accordingly.

To summarize, the influence of hot electrons on the nonlinear responses of our structure at NIR wavelength is significant. Higher *T _{e}* will result in lower saturation intensities and the incident intensity is limited. The maximum

*η*is only about $5 \times {10^{ - 6}}$ with

*I*= 25MW/cm

_{0}^{2}and

*T*= 550 K. At NIR wavelength, the main limitation on the TH generation of graphene in the double resonant structure, is the inherently small third order conductivities of graphene, which is 7 orders lower than that at THz frequencies. Meanwhile, the proposed structure may work as a low-threshold saturable absorber for large

_{e}*T*.

_{e}To highlight the enhancement of TH intensity arising from the double resonant structure, we compare *η* at NIR wavelength with those experimental results. In [15], graphene was covered on SiN waveguides and *η* is about 10^{−6}, which is on the same order as our results. However, the length of the waveguide is 100µm, about 10 times of our structure. In [17], for a few-layer graphene, *η* is only about 10^{−10}, though very intense incident light (2GW/cm^{2}) was used. In [18], a single layer graphene on a SiO_{2} substrate was studied, *η* is only about 10^{−12} and *I _{0}* is about 240MW/cm

^{2}. In [28], TH intensity of graphene is enhanced by a FP cavity structure, the enhancement is about 117-fold. In our structure, the enhancement factor can be estimated by Eq. (13). In Figs. 6(c) and 6(d), the linear mode intensity at ${\lambda _{FF}}$ and ${\lambda _{TH}} $ are both enhanced by about 100 times, thus the TH intensity will be enhanced by about 10

^{8}. Similar enhancement factor have also been predicted theoretically in a 2D photonic crystals [26].

## 4. Conclusion

In conclusion, in this work we insert graphene to the interface of two PCs with opposite topological properties. By carefully designing the two PCs, topological edge modes appear at both the FF and TH frequency. Strong nonlinear responses are demonstrated in the double resonant structure both at THz and NIR range. The Kerr nonlinearity and TH generation are studied as examples of potential applications of our structure. At THz frequency, the linear and nonlinear conductivities of graphene are inherently large. Increasing *E _{F}* will blue shifts the resonant frequency, while increasing

*I*will cause red shifts. After the two effects are balanced, TH conversion efficiency about 2.5% can be achieved with

_{0}*I*= 10MW/cm

_{0}^{2}. At NIR wavelength, the nonlinear conductivities are about 7 orders smaller than those at THz frequency. Meanwhile, high electron temperature will result in low saturation intensities and the incident intensity is limited. The maximum

*η*is only about $5 \times {10^{ - 6}}$ with

*I*= 25MW/cm

_{0}^{2}and

*T*= 550 K. The structure is more suitable to work as a low-threshold saturable absorber for large incident intensities and large electron temperatures. In the work, all the incident intensities are chosen carefully to be below the saturation intensity, and finite temperature effect is always considered. We believe the results will have potential applications for nonlinear THz devices and nanophotonics at NIR wavelength.

_{e}## Funding

National Natural Science Foundation of China (11704259, 11874269, 61805112); Natural Science Foundation of Guangdong Province (2018A030313198).

## References

**1. **R. W. Boyd, * Nonlinear Optics* (Elsevier, 2003).

**2. **G. Li, S. Zhang, and T. Zentgraf, “Nonlinear photonic metasurfaces,” Nat. Rev. Mater. **2**(5), 17010 (2017). [CrossRef]

**3. **M. Kauranen and A. V. Zayats, “Nonlinear plasmonics,” Nat. Photonics **6**(11), 737–748 (2012). [CrossRef]

**4. **D. Smirnova and Y. S. Kivshar, “Multipolar nonlinear nanophotonics,” Optica **3**(11), 1241–1255 (2016). [CrossRef]

**5. **K. S. Novoselov, “Electric Field Effect in Atomically Thin Carbon Films,” Science **306**(5696), 666–669 (2004). [CrossRef]

**6. **A. H. Castro Neto, F. Guinea, N. M. R. Peres, K. S. Novoselov, and A. K. Geim, “The electronic properties of graphene,” Rev. Mod. Phys. **81**(1), 109–162 (2009). [CrossRef]

**7. **Q. Bao and K. P. Loh, “Graphene Photonics, Plasmonics, and Broadband Optoelectronic Devices,” ACS Nano **6**(5), 3677–3694 (2012). [CrossRef]

**8. **F. Bonaccorso, Z. Sun, T. Hasan, and A. C. Ferrari, “Graphene photonics and optoelectronics,” Nat. Photonics **4**(9), 611–622 (2010). [CrossRef]

**9. **R. R. Nair, P. Blake, A. N. Grigorenko, K. S. Novoselov, T. J. Booth, T. Stauber, N. M. R. Peres, and A. K. Geim, “Fine Structure Constant Defines Visual Transparency of Graphene,” Science **320**(5881), 1308 (2008). [CrossRef]

**10. **A. N. Grigorenko, M. Polini, and K. S. Novoselov, “Graphene plasmonics,” Nat. Photonics **6**(11), 749–758 (2012). [CrossRef]

**11. **L. A. Falkovsky and S. S. Pershoguba, “Optical far-infrared properties of a graphene monolayer and multilayer,” Phys. Rev. B **76**(15), 153410 (2007). [CrossRef]

**12. **S. A. Mikhailov, “Non-linear electromagnetic response of graphene,” Europhys. Lett. **79**(2), 27002 (2007). [CrossRef]

**13. **E. Hendry, P. J. Hale, J. Moger, A. K. Savchenko, and S. A. Mikhailov, “Coherent Nonlinear Optical Response of Graphene,” Phys. Rev. Lett. **105**(9), 097401 (2010). [CrossRef]

**14. **T. Gu, N. Petrone, J. F. McMillan, A. van der Zande, M. Yu, G. Q. Lo, D. L. Kwong, J. Hone, and C. W. Wong, “Regenerative oscillation and four-wave mixing in graphene optoelectronics,” Nat. Photonics **6**(8), 554–559 (2012). [CrossRef]

**15. **K. Alexander, N. A. Savostianova, S. A. Mikhailov, B. Kuyken, and D. Van Thourhout, “Electrically Tunable Optical Nonlinearities in Graphene-Covered SiN Waveguides Characterized by Four-Wave Mixing,” ACS Photonics **4**(12), 3039–3044 (2017). [CrossRef]

**16. **Y. Q. An, F. Nelson, J. U. Lee, and A. C. Diebold, “Enhanced Optical Second-Harmonic Generation from the Current-Biased Graphene/SiO2/Si(001) Structure,” Nano Lett. **13**(5), 2104–2109 (2013). [CrossRef]

**17. **N. Kumar, J. Kumar, C. Gerstenkorn, R. Wang, H.-Y. Chiu, A. L. Smirl, and H. Zhao, “Third harmonic generation in graphene and few-layer graphite films,” Phys. Rev. B **87**(12), 121406 (2013). [CrossRef]

**18. **G. Soavi, G. Wang, H. Rostami, D. G. Purdie, D. D. Fazio, T. Ma, B. Luo, J. Wang, A. K. Ott, D. Yoon, S. A. Bourelle, J. E. Muench, I. Goykhman, S. D. Conte, M. Celebrano, A. Tomadin, M. Polini, G. Cerullo, and A. C. Ferrari, “Broadband, electrically tunable third-harmonic generation in graphene,” Nat. Nanotechnol. **13**(7), 583–588 (2018). [CrossRef]

**19. **T. Jiang, D. Huang, J. Cheng, X. Fan, Z. Zhang, Y. Shan, Y. Yi, Y. Dai, L. Shi, K. Liu, C. Zeng, J. Zi, J. E. Sipe, Y.-R. Shen, W.-T. Liu, and S. Wu, “Gate-tunable third-order nonlinear optical response of massless Dirac fermions in graphene,” Nat. Photonics **12**(7), 430–436 (2018). [CrossRef]

**20. **H. Zhang, S. Virally, Q. Bao, L. K. Ping, S. Massar, N. Godbout, and P. Kockaert, “Z-scan measurement of the nonlinear refractive index of graphene,” Opt. Lett. **37**(11), 1856–1858 (2012). [CrossRef]

**21. **N. Vermeulen, D. Castelló-Lurbe, M. Khoder, I. Pasternak, A. Krajewska, T. Ciuk, W. Strupinski, J. Cheng, H. Thienpont, and J. V. Erps, “Graphene’s nonlinear-optical physics revealed through exponentially growing self-phase modulation,” Nat. Commun. **9**(1), 2675 (2018). [CrossRef]

**22. **K. Alexander, N. A. Savostianova, S. A. Mikhailov, D. Van Thourhout, and B. Kuyken, “Gate-Tunable Nonlinear Refraction and Absorption in Graphene-Covered Silicon Nitride Waveguides,” ACS Photonics **5**(12), 4944–4950 (2018). [CrossRef]

**23. **J. W. You, J. You, M. Weismann, and N. C. Panoiu, “Double-resonant enhancement of third-harmonic generation in graphene nanostructures,” Philos. Trans. R. Soc., A **375**(2090), 20160313 (2017). [CrossRef]

**24. **D. Kundys, B. Van Duppen, O. P. Marshall, F. Rodriguez, I. Torre, A. Tomadin, M. Polini, and A. N. Grigorenko, “Nonlinear Light Mixing by Graphene Plasmons,” Nano Lett. **18**(1), 282–287 (2018). [CrossRef]

**25. **M. Weismann and N. C. Panoiu, “Theoretical and computational analysis of second- and third-harmonic generation in periodically patterned graphene and transition-metal dichalcogenide monolayers,” Phys. Rev. B **94**(3), 035435 (2016). [CrossRef]

**26. **T. Wang and X. Zhang, “Improved third-order nonlinear effect in graphene based on bound states in the continuum,” Photonics Res. **5**(6), 629 (2017). [CrossRef]

**27. **M. A. Vincenti, D. de Ceglia, M. Grande, A. D’Orazio, and M. Scalora, “Third-harmonic generation in one-dimensional photonic crystal with graphene-based defect,” Phys. Rev. B **89**(16), 165139 (2014). [CrossRef]

**28. **C. Beckerleg, T. J. Constant, I. Zeimpekis, S. M. Hornett, C. Craig, D. W. Hewak, and E. Hendry, “Cavity enhanced third harmonic generation in graphene,” Appl. Phys. Lett. **112**(1), 011102 (2018). [CrossRef]

**29. **J. L. Cheng, N. Vermeulen, and J. E. Sipe, “Third-order nonlinearity of graphene: Effects of phenomenological relaxation and finite temperature,” Phys. Rev. B **91**(23), 235320 (2015). [CrossRef]

**30. **J. L. Cheng, N. Vermeulen, and J. E. Sipe, “Numerical study of the optical nonlinearity of doped and gapped graphene: From weak to strong field excitation,” Phys. Rev. B **92**(23), 235307 (2015). [CrossRef]

**31. **S. A. Mikhailov, “Quantum theory of the third-order nonlinear electrodynamic effects of graphene,” Phys. Rev. B **93**(8), 085403 (2016). [CrossRef]

**32. **M. Xiao, Z. Q. Zhang, and C. T. Chan, “Surface Impedance and Bulk Band Geometric Phases in One-Dimensional Systems,” Phys. Rev. X **4**(2), 021017 (2014). [CrossRef]

**33. **K. H. Choi, C. W. Ling, K. F. Lee, Y. H. Tsang, and K. H. Fung, “Simultaneous multi-frequency topological edge modes between one-dimensional photonic crystals,” Opt. Lett. **41**(7), 1644 (2016). [CrossRef]

**34. **W. Gao, X. Hu, C. Li, J. Yang, Z. Chai, J. Xie, and Q. Gong, “Fano-resonance in one-dimensional topological photonic crystal heterostructure,” Opt. Express **26**(7), 8634 (2018). [CrossRef]

**35. **L. Wang, T. Wang, S. Zhang, P. Xie, and X. Zhang, “Larger enhancement in four-wave mixing from graphene embedded in one-dimensional photonic crystals,” J. Opt. Soc. Am. B **34**(9), 2000 (2017). [CrossRef]

**36. **W. Born and M. Emil, * Principles of Optics* (Elsevier, 2013).

**37. **T. Zhan, X. Shi, Y. Dai, X. Liu, and J. Zi, “Transfer matrix method for optics in graphene layers,” J. Phys.: Condens. Matter **25**(21), 215301 (2013). [CrossRef]

**38. **R. Yu, V. Pruneri, and F. J. G. de Abajo, “Resonant Visible Light Modulation with Graphene,” ACS Photonics **2**(4), 550–558 (2015). [CrossRef]

**39. **W. Zouaghi, D. Voß, M. Gorath, N. Nicoloso, and H. G. Roskos, “How good would the conductivity of graphene have to be to make single-layer-graphene metamaterials for terahertz frequencies feasible?” Carbon **94**, 301–308 (2015). [CrossRef]

**40. **G. X. Ni, A. S. McLeod, Z. Sun, L. Wang, L. Xiong, K. W. Post, S. S. Sunku, B.-Y. Jiang, J. Hone, C. R. Dean, M. M. Fogler, and D. N. Basov, “Fundamental limits to graphene plasmonics,” Nature **557**(7706), 530–533 (2018). [CrossRef]

**41. **M. Jablan, H. Buljan, and M. Soljačić, “Plasmonics in graphene at infrared frequencies,” Phys. Rev. B **80**(24), 245435 (2009). [CrossRef]

**42. **K. O’Brien, H. Suchowski, J. Rho, A. Salandrino, B. Kante, X. Yin, and X. Zhang, “Predicting nonlinear properties of metamaterials from the linear response,” Nat. Mater. **14**(4), 379–383 (2015). [CrossRef]

**43. **M. S. Nezami, D. Yoo, G. Hajisalem, S.-H. Oh, and R. Gordon, “Gap Plasmon Enhanced Metasurface Third-Harmonic Generation in Transmission Geometry,” ACS Photonics **3**(8), 1461–1467 (2016). [CrossRef]

**44. **A. Marini, J. D. Cox, and F. J. G. de Abajo, “Theory of graphene saturable absorption,” Phys. Rev. B **95**(12), 125408 (2017). [CrossRef]

**45. **J. Guo, L. Jiang, Y. Jia, X. Dai, Y. Xiang, and D. Fan, “Low threshold optical bistability in one-dimensional gratings based on graphene plasmonics,” Opt. Express **25**(6), 5972–5981 (2017). [CrossRef]

**46. **H. A. Hafez, S. Kovalev, J.-C. Deinert, Z. Mics, B. Green, N. Awari, M. Chen, S. Germanskiy, U. Lehnert, J. Teichert, Z. Wang, K.-J. Tielrooij, Z. Liu, Z. Chen, A. Narita, K. Müllen, M. Bonn, M. Gensch, and D. Turchinovich, “Extremely efficient terahertz high-harmonic generation in graphene by hot Dirac fermions,” Nature **561**(7724), 507–511 (2018). [CrossRef]

**47. **Z. Mics, K.-J. Tielrooij, K. Parvez, S. A. Jensen, I. Ivanov, X. Feng, K. Müllen, M. Bonn, and D. Turchinovich, “Thermodynamic picture of ultrafast charge transport in graphene,” Nat. Commun. **6**(1), 7655 (2015). [CrossRef]

**48. **G. Soavi, G. Wang, H. Rostami, A. Tomadin, O. Balci, I. Paradeisanos, and E. a. A. Pogna, G. Cerullo, E. Lidorikis, M. Polini, and A. C. Ferrari, “Hot electrons modulation of third harmonic generation in graphene,” arXiv:1903.00989 [cond-mat.mes-hall] (2019).

**49. **S. A. Mikhailov, “Theory of the strongly nonlinear electrodynamic response of graphene: A hot electron model,” Phys. Rev. B **100**(11), 115416 (2019). [CrossRef]

**50. **K.-J. Tielrooij, N. C. H. Hesp, A. Principi, M. B. Lundeberg, E. A. A. Pogna, L. Banszerus, Z. Mics, M. Massicotte, P. Schmidt, D. Davydovskaya, D. G. Purdie, I. Goykhman, G. Soavi, A. Lombardo, K. Watanabe, T. Taniguchi, M. Bonn, D. Turchinovich, C. Stampfer, A. C. Ferrari, G. Cerullo, M. Polini, and F. H. L. Koppens, “Out-of-plane heat transfer in van der Waals stacks through electron–hyperbolic phonon coupling,” Nat. Nanotechnol. **13**(1), 41–46 (2018). [CrossRef]