Improved third-order nonlinear effect in graphene based on bound states in the continuum

Tiecheng Wang; Xiangdong Zhang

doi:10.1364/PRJ.5.000629

Abstract

The scattering matrix theory has been developed to calculate the third-order nonlinear effect in sphere-graphene-slab structures. By designing structural parameters, we have demonstrated that the incident electromagnetic wave can be well confined in the graphene in these structures due to the formation of a bound state in the continuum (BIC) of radiation modes. Based on such a bound state, third-harmonic (TH) generation and four-wave mixing (FWM) have been studied. It is found that the efficiency of TH generation in monolayer graphene can be enhanced about 7 orders of magnitude. It is interesting that we can design structure parameters to make all beams (the pump beam, probe beam, and generated FWM signal) be BICs at the same time. In such a case, the efficiency of FWM in monolayer graphene can be enhanced about 9 orders of magnitude. Both the TH and FWM signals are sensitive to the wavelength, and possess high Q factors, which exhibit very good monochromaticity. By taking suitable BICs, the selective generation of TH and FWM signals for S- and P-polarized waves can also be realized, which is beneficial for the design of optical devices.

Keywords

(190.2620) Harmonic generation and mixing.parametric processes (190.4410) Nonlinear optics

1. INTRODUCTION

In recent years, nonlinear optical properties of graphene have attracted much attention because recent investigations have shown that the third-order nonlinear optical response of the single-atom-layer graphene is particularly strong ^[1–5]. Many studies have focused on such an aspect, such as saturable absorption ^[6–11], optical limiting ^[12,13], two-photon absorption ^[14], four-wave mixing (FWM) ^[1,15,16], and third-harmonic (TH) generation ^[17–24]. This allows one to employ graphene in active photonic devices with improved functionality ^[3,18,19]. In these applications, the strength of interaction between graphene and electromagnetic (EM) waves plays a central role. In general, such an interaction is weak due to the single atom thickness of graphene. The weak interactions block efficient generation of nonlinear effects and potential applications of graphene in optoelectronic devices. Thus, various methods to improve the interaction between graphene and EM waves have been proposed ^{[5–13,25–27]}. For example, TH generation may be dramatically enhanced in monolayer graphene by inserting it in properly designed 1D photonic crystals (PCs) ^[28] or using ultrastrong localized surface plasmon resonances ^[29–31].

On the other hand, analogous to the localized electrons with energy larger than their potential barriers, light bound states in the continuum (BIC) have been observed experimentally in recent years ^[32–34]. The BICs are known as embedded trapped modes, which correspond to discrete eigenvalues coexisting with extended modes of a continuous spectrum ^[35–37]. For example, the BICs for the EM wave have been shown to exist in the dielectric gratings ^[35], waveguide structures ^{[32,36,38–40]}, the surface of the object ^[37,41], PC slabs ^{[33,34,42,43]}, and some open subwavelength nanostructures ^[44–46]. Photonic nanostructures with such BICs possess many advantages. First, they may be more easily implemented compared with the electron system. In the electronic quantum system, some artificially designed potentials cannot be readily realized. Second, they may have a variety of potential applications, such as biosensing, perfect filters, and so on.

Motivated by these investigations, in this work we explore the possibility to improve the efficiency of the third-order nonlinear effect in graphene using BICs. Two kinds of third-order nonlinear effect, TH generation and FWM, are considered. The method to calculate these nonlinear effects in sphere-graphene-slab structures and multilayer structures has been developed. We find that the efficiency of the TH generation and FWM can be dramatically enhanced in monolayer graphene based on BICs. The generated signals possess very good monochromaticity. The physical origin for such a phenomenon has also been analyzed. The rest of this paper is organized as follows: We will describe the theory and method in Section 2, present the results and discussions in Section 3, and give a summary in Section 4.

Sign up for Photonics Research TOC. Get the latest issue of Photonics Research delivered right to you！Sign up now

2. THEORY AND METHOD

We consider a three-layer structure consisting of monolayer dielectric spheres with a square lattice, a monolayer graphene, and a dielectric slab, as shown in Figs. 1(a) and 1(b); the monolayer graphene is placed at the interface between the monolayer dielectric spheres and the dielectric slab. The distance between two neighbor spheres is taken as $a$ , the radius, relative permittivity, and relative permeability of spheres are denoted by $r_{s}$ , $ϵ_{s}$ , and $μ_{s}$ , respectively; the spheres are immersed in the air. The thickness, relative permittivity, and relative permeability of the dielectric slab are represented by $d$ , $ϵ_{d}$ , and $μ_{d}$ , respectively. All these materials are taken to be nonmagnetic $μ_{s} = μ_{d} = 0$ . The thickness, conductivity, and third-order susceptibility of the monolayer graphene are marked by $d_{g}$ , $σ$ , and $χ^{(3)}$ ; the nonlinear surface current due to the third-order nonlinear effects of the monolayer graphene at the graphene location is ${\vec{J}}^{(3)} (3 ω) = - i d_{g} 3 ω {\vec{P}}^{(3)} (3 ω)$ , where $3 ω$ and ${\vec{P}}^{(3)} (3 ω)$ represent the angular frequency of the TH EM wave and the third-order polarization density.

Figure 1.(a) Diagram of the sphere-graphene-slab structure and the TH generation process. The spheres are arranged in a square lattice with lattice constant $a$ ; the thickness of the dielectric slab is $d$ , and the monolayer graphene is put at the interface between the monolayer dielectric spheres and the dielectric slab. When a plane wave with angular frequency $ω$ is incident on the structure at the incident angle $θ_{i}$ , the transmitted and reflected TH fields are generated due to the third-order nonlinear effect of graphene. Inset shows graphene layer with up medium denoted by 1 and down medium marked by 2. (b) Schematic of the degenerate FWM process in the three-layer structure. The pump beam at frequency $ω_{1}$ and the probe beam at frequency $ω_{2}$ are incident on the structure with arbitrary incident angles $θ_{1}$ and $θ_{2}$ . Because of the third nonlinear effect of the graphene, the transmitted and reflected beams at FWM frequency $ω_{3} = 2 ω_{1} - ω_{2}$ are generated with angle $θ_{3}$ . The angle $θ_{3}$ is determined by the in-plane phase-matching condition.

The surface conductivity of the graphene $σ$ is governed by the Kubo formula ^[47–54]; we take the conductivity as follows:

The surface conductivity $σ$ contains the intraband contribution $σ_{intra}$ and the interband contribution $σ_{inter}$ ; here, $e$ , $ℏ$ , $k_{B}$ , and $T$ are the electronic charge, Planck constant, Boltzmann constant, and temperature, respectively; $f_{d} (ϵ)$ , $E_{F}$ , and $γ$ represent the Fermi distribution, Fermi energy, and relaxation rate. We take the temperature and relaxation rate as $T = 300 K$ and $γ = 1.3 meV$ ^[55] in the following calculation. In the THz regime, when $E_{F} ≫ k_{B} T$ and $2 E_{F} > ℏ ω$ , the conductivity of graphene can be well approximated by the Drude-like expression $σ = \frac{σ_{0}}{π} \frac{4 E_{F}}{ℏ γ - j ℏ ω}$ ^[56], where $σ_{0} = e^{2} / (4 ℏ)$ is the well-known universal conductivity. Thus, the conductivity $σ$ can be tuned by changing the Fermi level via chemical doping or applying a gate voltage.

Based on the recent theoretical and experimental studies, the only nonzero components of the third-order susceptibility tensor are $χ_{x x x}^{(3)}$ , $χ_{y y y}^{(3)}$ , and $χ_{z z z}^{(3)}$ due to the isotropic nature. We assume $χ_{x x x}^{(3)} = χ_{y y y}^{(3)} = χ^{(3)}$ and neglect $χ_{z z z}^{(3)}$ for $χ_{z z z}^{(3)} ≪ χ_{x x x}^{(3)}$ ^[19] and take the susceptibility as $χ^{(3)} = 1.4 \times 10^{- 16} m^{2} / V^{2}$ ^[20,21] compatible with the recent study. When the irradiance on graphene is increased, the saturation effect ^[57] will be enhanced; the third-order susceptibility is modified and taken as $χ^{(3)} S$ , where $S = \frac{1}{1 + ({| E_{x} |}^{2} + {| E_{y} |}^{2}) / {| E_{sat} |}^{2}}$ is the saturation coefficient ^[58,59], $E_{x}$ and $E_{y}$ are the electric field on graphene. The saturation field $E_{sat}$ is taken from the experimental data ^[58]. In the following, we will provide the calculation method for TH generation and FWM in such a structure.

A. Theory for the TH Generation in the Sphere-Graphene-Slab Structure

When a plane wave ${\vec{E}}_{i} (\vec{r}, ω)$ with angular frequency $ω$ is incident on the structure at the incident angle $θ_{i}$ , as shown in Fig. 1(a), the transmitted and reflected TH fields are generated due to the third-order nonlinear effect of graphene. The fundamental frequency (FF) field [ ${\vec{E}}_{x} (\vec{r}, ω)$ , ${\vec{E}}_{y} (\vec{r}, ω)$ , and ${\vec{E}}_{z} (\vec{r}, ω)$ ] can be calculated by the layer-multiple-scattering method, of which details of the calculation process are described in Ref. ^[60]. Thus, the third-order polarization densities of graphene are expressed as $P_{x}^{(3)} (3 ω) = χ^{(3)} {[E_{x} {(\vec{r}, ω) |}_{z = 0}]}^{3}$ , $P_{y}^{(3)} (3 ω) = χ^{(3)} {[E_{y} {(\vec{r}, ω) |}_{z = 0}]}^{3}$ , and $P_{z}^{(3)} (3 ω) = 0$ ; here, $E_{x} {(\vec{r}, ω) |}_{z = 0}$ and $E_{y} {(\vec{r}, ω) |}_{z = 0}$ represent $x$ and $y$ components of the FF electric field on the graphene, which is located at the interface $z = 0$ ; the corresponding surface currents are expressed as $J_{x}^{(3)} (3 ω) = - 3 i d_{g} ω P_{x}^{(3)} (3 ω)$ and $J_{y}^{(3)} (3 ω) = - 3 i d_{g} ω P_{y}^{(3)} (3 ω)$ .

In the following, we present the calculated method for the TH generation in such a structure. Because the sphere layer and dielectric slab are considered as the linear materials, at the TH frequency we only consider the scattering matrix of graphene. We first consider that a monolayer graphene laid at the interface $z = 0$ , as shown in inset of Fig. 1(a), is radiated by a plane wave at the FF from the $x - z$ plane; the up/down medium of the graphene is denoted by 1/2. The TH electric and magnetic fields at the up/down medium of the graphene are marked by ${\vec{E}}_{1} (\vec{r}, 3 ω) / {\vec{E}}_{2} (\vec{r}, 3 ω)$ and ${\vec{H}}_{1} (\vec{r}, 3 ω) / {\vec{H}}_{2} (\vec{r}, 3 ω)$ with $x$ and $z$ components of the wave vector marked by $q_{x}$ and $q_{1 z}$ ( $q_{x}$ and $q_{2 z}$ ), the $x$ components of the wave vector $q_{x}$ in the up and down media are equal. At the boundary, we have $\vec{n} \times {[{\vec{E}}_{2} (\vec{r}, 3 ω) - {\vec{E}}_{1} (\vec{r}, 3 ω)] |}_{z = 0} = 0,$ (4a) $\vec{n} \times {[{\vec{H}}_{2} (\vec{r}, 3 ω) - {\vec{H}}_{1} (\vec{r}, 3 ω)] |}_{z = 0} = σ {\vec{E}}_{1} {(\vec{r}, 3 ω) |}_{z = 0} + {\vec{J}}^{(3)} (3 ω),$ (4b)where $\vec{n}$ represents the direction of the normal unit vector at the interface. The electric fields at the TH frequency in the up and down media can be expressed as ${\vec{E}}_{1} (\vec{r}, 3 ω) = {\vec{E}}_{1}^{i} (3 ω) e^{i (q_{x} x + q_{1 z} z - 3 ω t)} + {\vec{E}}_{1}^{r} (3 ω) e^{i (q_{x} x - q_{1 z} z - 3 ω t)},$ (5a) ${\vec{E}}_{2} (\vec{r}, 3 ω) = {\vec{E}}_{2}^{t} (3 ω) e^{i (q_{x} x + q_{2 z} z - 3 ω t)} .$ (5b)

Here ${\vec{E}}_{1}^{i} (3 ω)$ , ${\vec{E}}_{1}^{r} (3 ω)$ , and ${\vec{E}}_{2}^{t} (3 ω)$ represent the expanded coefficients of the electric TH fields. Inserting Eq. (5) into the boundary condition Eq. (4), we obtain the relations between the electric fields up and down the graphene: $E_{2 x}^{t} (3 ω) = T_{x}^{(1, 2)} E_{1 x}^{i} (3 ω) + C T_{x}^{(1, 2)},$ (6a) $E_{2 y}^{t} (3 ω) = T_{y}^{(1, 2)} E_{1 y}^{i} (3 ω) + C T_{y}^{(1, 2)},$ (6b) $E_{1 x}^{r} (3 ω) = R_{x}^{(1, 2)} E_{1 x}^{i} (3 ω) + C R_{x}^{(1, 2)},$ (6c) $E_{1 y}^{r} (3 ω) = R_{y}^{(1, 2)} E_{1 y}^{i} (3 ω) + C R_{y}^{(1, 2)},$ (6d)where $T_{x}^{(1, 2)}$ and $T_{y}^{(1, 2)}$ represent the transmission coefficients, and $T_{y}^{(1, 2)}$ and $R_{y}^{(1, 2)}$ are reflection coefficients, which are given by $T_{x}^{(1, 2)} = \frac{2 q_{2 z} / ϵ_{2}}{q_{1 z} / ϵ_{1} + q_{2 z} / ϵ_{2} + q_{1 z} q_{2 z} σ / (ω ϵ_{0} ϵ_{1} ϵ_{2})},$ (7a) $T_{y}^{(1, 2)} = \frac{2 q_{1 z}}{q_{1 z} + q_{2 z} + ω μ_{0} σ},$ (7b) $R_{x}^{(1, 2)} = \frac{q_{2 z} / ϵ_{2} - q_{1 z} / ϵ_{1} - q_{1 z} q_{2 z} σ / (ω ϵ_{0} ϵ_{1} ϵ_{2})}{q_{1 z} / ϵ_{1} + q_{2 z} / ϵ_{2} + q_{1 z} q_{2 z} σ / (ω ϵ_{0} ϵ_{1} ϵ_{2})},$ (7c) $R_{y}^{(1, 2)} = \frac{q_{1 z} - q_{2 z} - ω μ_{0} σ}{q_{1 z} + q_{2 z} + ω μ_{0} σ} .$ (7d)

Here $C T_{x}^{(1, 2)}$ , $C T_{y}^{(1, 2)}$ , $C R_{x}^{(1, 2)}$ , and $C R_{y}^{(1, 2)}$ are the nonlinear terms originating from the third-order nonlinear effect of graphene and are expressed as $C T_{x}^{(1, 2)} = \frac{- J_{x}^{(3)} (3 ω) q_{1 z} q_{2 z} σ / (ω ϵ_{0} ϵ_{1} ϵ_{2})}{q_{1 z} / ϵ_{1} + q_{2 z} / ϵ_{2} + q_{1 z} q_{2 z} σ / (ω ϵ_{0} ϵ_{1} ϵ_{2})},$ (8a) $C T_{y}^{(1, 2)} = \frac{- J_{y}^{(3)} (3 ω) ω μ_{0}}{q_{1 z} + q_{2 z} + ω μ_{0} σ},$ (8b) $C R_{x}^{(1, 2)} = \frac{- J_{x}^{(3)} (3 ω) q_{1 z} q_{2 z} σ / (ω ϵ_{0} ϵ_{1} ϵ_{2})}{q_{1 z} / ϵ_{1} + q_{2 z} / ϵ_{2} + q_{1 z} q_{2 z} σ / (ω ϵ_{0} ϵ_{1} ϵ_{2})} = C T_{x}^{(1, 2)},$ (8c) $C R_{y}^{(1, 2)} = \frac{- J_{y}^{(3)} (3 ω) ω μ_{0}}{q_{1 z} + q_{2 z} + ω μ_{0} σ} = C T_{y}^{(1, 2)},$ (8d)where $ϵ_{0}$ and $μ_{0}$ are the permittivity and permeability of the vacuum, respectively. In the following, we extend the layer-multiple-scattering method to include a nonlinear graphene layer. The electric fields outside the graphene layer are expressed as ${\vec{E}}_{in}^{+} (\vec{r}, 3 ω)$ , ${\vec{E}}_{tr}^{+} (\vec{r}, 3 ω)$ , and ${\vec{E}}_{rf}^{-} (\vec{r}, 3 ω)$ , which represent the incident wave along the $+ z$ direction (propagating downward) and the corresponding transmitted wave along $+ z$ direction and reflected waves along $- z$ direction (propagating upward). The corresponding waves with opposite directions are marked by ${\vec{E}}_{in}^{-} (\vec{r}, 3 ω)$ , ${\vec{E}}_{tr}^{-} (\vec{r}, 3 ω)$ , and ${\vec{E}}_{rf}^{+} (\vec{r}, 3 ω)$ , respectively. These waves can be expanded as a set of plane waves: ${\vec{E}}_{in (tr, rf)}^{\pm} (\vec{r}, 3 ω) = \sum_{l = 1}^{2} \sum_{m n} E_{in (tr, rf), l, m n}^{\pm} (3 ω) e^{i (k_{m n x} x + k_{m n y} y \pm k_{m n z} z)} {\hat{u}}_{l},$ (9)where $l = 1, 2$ represent the $x$ , $y$ coordinates; $(k_{m n x}, k_{m n y}) = (k_{x}, k_{y}) + (g_{m n x}, g_{m n y})$ is the Bragg wave vector and $k_{m n z} = \sqrt{q^{2} - k_{m n x}^{2} - k_{m n y}^{2}}$ . Here, $q$ is the wave vector, and $(k_{x}, k_{y})$ is the reduced wave vector, which lies in the surface Brillouin zone; $(g_{m n x}, g_{m n y}) = m {\vec{b}}_{1} + n {\vec{b}}_{2}$ is the general reciprocal vector, and ${\vec{b}}_{1}$ and ${\vec{b}}_{2}$ represent the reciprocal lattice vector. Outside the graphene, the relationships between the incident and the transmitted (reflected) waves are $E_{tr, l}^{s} (3 ω) = \sum_{l^{'} = 1}^{2} N_{l l^{'}}^{s s} E_{in, l^{'}}^{s} (3 ω) + D_{l}^{s s},$ (10a) $E_{r f, l}^{- s} (3 ω) = \sum_{l^{'} = 1}^{2} N_{l l^{'}}^{- s s} E_{in, l^{'}}^{s} (3 ω) + D_{l}^{- s s},$ (10b)where $s = + (-)$ denotes a wave propagating downward (upward); the matrixes $N^{s s^{'}}$ and $D^{s s^{'}}$ are given by $N_{l l^{'}}^{s s^{'}} = (\begin{matrix} N_{x}^{s s^{'}} \cos^{2} ϕ + N_{y}^{s s^{'}} \sin^{2} ϕ & (N_{x}^{s s^{'}} - N_{y}^{s s^{'}}) \sin ϕ \cos ϕ \\ (N_{x}^{s s^{'}} - N_{y}^{s s^{'}}) \sin ϕ \cos ϕ & N_{x}^{s s^{'}} \sin^{2} ϕ + N_{y}^{s s^{'}} \cos^{2} ϕ \end{matrix}),$ (11a) $D_{l}^{s s^{'}} = (\begin{matrix} \cos ϕ D_{x}^{s s^{'}} - \sin ϕ D_{y}^{s s^{'}} \\ \sin ϕ D_{x}^{s s^{'}} + \cos ϕ D_{y}^{s s^{'}} \end{matrix}),$ (11b)where $ϕ$ is the azimuthal angle of $q_{‖}$ with respect to the $x$ axis and $q_{‖}$ is the parallel component of $q$ . Here $N_{l}^{++} = T_{l}^{(1, 2)}$ , $N_{l}^{+ -} = R_{l}^{(2, 1)}$ , $N_{l}^{- +} = R_{l}^{(1, 2)}$ , and $N_{l}^{- -} = T_{l}^{(2, 1)}$ represent the elements in matrixes $N^{s s^{'}}$ ; $D_{l}^{++} = C T_{l}^{(1, 2)}$ , $D_{l}^{+ -} = C R_{l}^{(2, 1)}$ , $D_{l}^{- +} = C R_{l}^{(1, 2)}$ , and $D_{l}^{- -} = C T_{l}^{(2, 1)}$ are the elements in matrixes $D_{l}^{s s^{'}}$ . We can further calculate scattering matrixes $Q_{2}^{η} (η = I, II, III, IV)$ and $C_{2}^{κ} (κ = I, II)$ for the monolayer graphene, the subscript number 2 in these matrixes denotes the graphene layer (because it is the second layer in the structure). Similarly, we can also obtain $Q_{1}^{η}$ and $C_{1}^{κ}$ (the scattering matrixes of the spheres layer), $Q_{3}^{η}$ and $C_{3}^{κ}$ (the scattering matrixes of the dielectric layer). The incident fields from two sides of graphene can be expanded in plane waves ${\vec{E}}_{1}^{+} (\vec{r}, 3 ω) = \sum_{l = 1}^{2} \sum_{m n} E_{1, l, m n}^{+} (3 ω) e^{i (k_{m n x} x + k_{m n y} y + k_{m n z} z)} {\hat{u}}_{l}$ and ${\vec{E}}_{2}^{-} (\vec{r}, 3 ω) = \sum_{l = 1}^{2} \sum_{m n} E_{2, l, m n}^{-} (3 ω) e^{i (k_{m n x} x + k_{m n y} y - k_{m n z} z)} {\hat{u}}_{l}$ , the output fields ${\vec{E}}_{2}^{+} (\vec{r}, 3 ω)$ and ${\vec{E}}_{1}^{-} (\vec{r}, 3 ω)$ from two sides of graphene can be expanded in a similar way; here ${\vec{E}}_{1}^{+} (\vec{r}, 3 ω) = {\vec{E}}_{in}^{+} (\vec{r}, 3 ω)$ , ${\vec{E}}_{2}^{-} (\vec{r}, 3 ω) = {\vec{E}}_{in}^{-} (\vec{r}, 3 ω)$ , ${\vec{E}}_{2}^{+} (\vec{r}, 3 ω) = {\vec{E}}_{tr}^{+} (\vec{r}, 3 ω) + {\vec{E}}_{rf}^{-} (\vec{r}, 3 ω)$ , and ${\vec{E}}_{1}^{-} (\vec{r}, 3 ω) = {\vec{E}}_{rf}^{+} (\vec{r}, 3 ω) + {\vec{E}}_{tr}^{-} (\vec{r}, 3 ω)$ . Then the electrical fields outside the graphene can be expressed as $(\begin{matrix} E_{2}^{+} (3 ω) \\ E_{1}^{-} (3 ω) \end{matrix}) = (\begin{matrix} Q_{2}^{I} & Q_{2}^{II} \\ Q_{2}^{III} & Q_{2}^{IV} \end{matrix}) (\begin{matrix} E_{1}^{+} (3 ω) \\ E_{2}^{-} (3 ω) \end{matrix}) + (\begin{matrix} C_{2}^{I} \\ C_{2}^{II} \end{matrix}),$ (12)where the electrical fields up and down the graphene are written in the matrix forms $E_{1}^{\pm} (3 ω) = {(\begin{matrix} \dots & E_{1, x, m n}^{\pm} (3 ω) & E_{1, y, m n}^{\pm} (3 ω) & \dots \end{matrix})}^{T}$ and $E_{2}^{\pm} (3 ω) = {(\begin{matrix} \dots & E_{2, x, m n}^{\pm} (3 ω) & E_{2, y, m n}^{\pm} (3 ω) & \dots \end{matrix})}^{T}$ ; here, the superscript $T$ represents matrix transposition. In a similar way, the electric fields up the sphere layer and down the dielectric layer can also be expanded in plane waves ${\vec{E}}_{0}^{\pm} (\vec{r}, 3 ω) = \sum_{l = 1}^{2} \sum_{m n} E_{0, l, m n}^{\pm} (3 ω) e^{i (k_{m n x} x + k_{m n y} y + k_{m n z} z)} {\hat{u}}_{l}$ and ${\vec{E}}_{3}^{\pm} (\vec{r}, 3 ω) = \sum_{l = 1}^{2} \sum_{m n} E_{3, l, m n}^{\pm} (3 ω) e^{i (k_{m n x} x + k_{m n y} y + k_{m n z} z)} {\hat{u}}_{l}$ ; we also define the corresponding matrixes $E_{0}^{\pm} (3 ω) = {(\begin{matrix} \dots & E_{0, x, m n}^{\pm} (3 ω) & E_{0, y, m n}^{\pm} (3 ω) & \dots \end{matrix})}^{T}$ and $E_{3}^{\pm} (3 ω) = {(\begin{matrix} \dots & E_{3, x, m n}^{\pm} (3 ω) & E_{3, y, m n}^{\pm} (3 ω) & \dots \end{matrix})}^{T}$ .

In the following, we attempt to calculate the scattering matrixes of the three-layer structure, which are marked by $Q_{1, 2, 3}^{η}$ and $C_{1, 2, 3}^{κ}$ based on the scattering matrixes of every layer. First, we calculate the scattering matrixes of the first two successive layers (the spheres layer and graphene layer), which are denoted by $Q_{1, 2}^{η}$ and $C_{1, 2}^{κ}$ . The matrixes $Q_{1, 2}^{η}$ and $C_{1, 2}^{κ}$ can be obtained by combining the matrixes $Q_{1}^{η}$ and $C_{1}^{κ}$ of the sphere layer and the matrixes $Q_{2}^{η}$ and $C_{2}^{κ}$ of the graphene layer: $Q_{1,2}^{I} = Q_{2}^{I} {[1 - Q_{1}^{II} Q_{2}^{III}]}^{- 1} Q_{1}^{I},$ (13a) $Q_{1, 2}^{II} = Q_{2}^{II} + Q_{2}^{I} Q_{1}^{II} {[1 - Q_{2}^{III} Q_{1}^{II}]}^{- 1} Q_{2}^{IV},$ (13b) $Q_{1, 2}^{III} = Q_{1}^{III} + Q_{1}^{IV} Q_{2}^{III} {[1 - Q_{1}^{II} Q_{2}^{III}]}^{- 1} Q_{1}^{I},$ (13c) $Q_{1,2}^{IV} = Q_{1}^{IV} {[1 - Q_{2}^{III} Q_{1}^{II}]}^{- 1} Q_{2}^{IV},$ (13d) $C_{1, 2}^{I} = Q_{2}^{I} {[1 - Q_{1}^{II} Q_{2}^{III}]}^{- 1} [Q_{1}^{II} C_{2}^{II} + C_{1}^{I}] + C_{2}^{I},$ (14a) $C_{1,2}^{II} = Q_{n}^{IV} {[1 - Q_{2}^{III} Q_{1}^{II}]}^{- 1} [Q_{2}^{III} C_{1}^{I} + C_{2}^{II}] + C_{1}^{II},$ (14b)where the matrixes $Q_{1, 2}^{η}$ in Eq. (13) are consistent with those in Ref. ^[60]; $C_{1, 2}^{κ}$ are the corresponding nonlinear terms, which are provided for the first time. Then, in a similar way with that given above, we can obtain the scattering matrixes $Q_{1, 2, 3}^{η}$ and $C_{1, 2, 3}^{κ}$ of the three successive layers based on the scattering matrixes $Q_{1, 2}^{η}$ and $C_{1, 2}^{κ}$ of the first two successive layers and the scattering matrixes $Q_{3}^{η}$ and $C_{3}^{κ}$ of the dielectric layer. Thus, we can obtain the TH fields upside $[E_{3}^{+} (3 ω)]$ and downside $[E_{0}^{-} (3 ω)]$ the three-layer structure from the following relation: $(\begin{matrix} E_{3}^{+} (3 ω) \\ E_{0}^{-} (3 ω) \end{matrix}) = (\begin{matrix} Q_{1,2, 3}^{I} & Q_{1,2, 3}^{II} \\ Q_{1,2, 3}^{III} & Q_{1,2, 3}^{IV} \end{matrix}) (\begin{matrix} E_{0}^{+} (3 ω) \\ E_{3}^{-} (3 ω) \end{matrix}) + (\begin{matrix} C_{1,2, 3}^{I} \\ C_{1,2, 3}^{II} \end{matrix}) = (\begin{matrix} Q_{1,2, 3}^{I} & Q_{1,2, 3}^{II} \\ Q_{1,2, 3}^{III} & Q_{1,2, 3}^{IV} \end{matrix}) (\begin{matrix} 0 \\ 0 \end{matrix}) + (\begin{matrix} C_{1,2, 3}^{I} \\ C_{1,2, 3}^{II} \end{matrix}) .$ (15)

Based on the calculated TH fields in Eq. (15), the electric field intensities of the transmitted $[I_{t} (3 ω)]$ and reflected $[I_{r} (3 ω)]$ TH fields radiated from the three-layer structure can be calculated by $I_{t} (3 ω) = \frac{1}{2} ϵ_{0} \sum_{l, m n} [E_{3, l, m n}^{+} (3 ω)] {[E_{3, l, m n}^{+} (3 ω)]}^{*},$ (16a) $I_{r} (3 ω) = \frac{1}{2} ϵ_{0} \sum_{l, m n} [E_{0, l, m n}^{-} (3 ω)] {[E_{0, l, m n}^{-} (3 ω)]}^{*} .$ (16b)

The electric field intensity of the incident FF field $I_{i} (ω)$ can also be calculated based on the known incident FF electric field ${\vec{E}}_{i}^{+} (\vec{r}, ω) = \sum_{l = 1}^{2} \sum_{m n} E_{i, l, m n}^{+} (ω) e^{i (k_{m n x}^{i} x + k_{m n y}^{i} y + k_{m n z}^{i} z)} {\hat{u}}_{l}$ ; it is $I_{i} (ω) = \frac{1}{2} ϵ_{0} \sum_{l, m n} [E_{i, l, m n}^{+} (ω)] {[E_{i, l, m n}^{+} (ω)]}^{*} .$ (16c)

In the calculating the sums of Eq. (16) waves with imaginary $z$ component of wave vector (decaying waves) are dropped, as they do not transfer energy along the $z$ direction.

Then, the conversion efficiencies for the transmitted (F1), reflected (F2), and total (F3) TH fields can be obtained: $F 1 = \frac{I_{t} (3 ω)}{I_{i} (ω)}, F 2 = \frac{I_{r} (3 ω)}{I_{i} (ω)}, F 3 = \frac{I_{t} (3 ω) + I_{r} (3 ω)}{I_{i} (ω)} .$ (17)

In order to compare the TH fields generated in the three-layer structure with those in the freestanding graphene, the enhancements of the transmitted (G1), reflected (G2), and total (G3) TH field are introduced as $G 1 = \frac{I_{t} (3 ω)}{I_{0 t} (3 ω)}, G 2 = \frac{I_{r} (3 ω)}{I_{0 r} (3 ω)}, G 3 = \frac{I_{t} (3 ω) + I_{r} (3 ω)}{I_{0 t} (3 ω) + I_{0 r} (3 ω)},$ (18)where $I_{0 t} (3 ω)$ and $I_{0 r} (3 ω)$ represent the electric field intensities of the transmitted and reflected TH fields generated in the freestanding graphene. Based on the above relations, we can obtain the conversion efficiencies and the enhancements for TH fields using numerical calculations.

B. Theory for the FWM in the Sphere-Graphene-Slab Structure

Now we turn to the FWM of the three-layer structure shown in Fig. 1(b). This structure is illuminated by a pump beam ${\vec{E}}_{1}^{+} (\vec{r}, ω_{1})$ at frequency $ω_{1}$ and incident angle $θ_{1}$ ; it is also radiated by a probe beam ${\vec{E}}_{2}^{+} (\vec{r}, ω_{2})$ at frequency $ω_{2}$ and incident angle $θ_{2}$ . Due to the third-order effect of the graphene, the transmitted beam ${\vec{E}}_{3}^{+} (\vec{r}, ω_{3})$ and reflected beam ${\vec{E}}_{3}^{-} (\vec{r}, ω_{3})$ at the FWM frequency $ω_{3} = 2 ω_{1} - ω_{2}$ are generated; here, the subscripts 1, 2, and 3 represent the three waves. The fields for them are expressed as ${\vec{E}}_{1 (2,3)} (\vec{r}, ω_{1 (2,3)}) = {\vec{E}}_{1 (2,3)}^{+} (\vec{r}, ω_{1 (2,3)}) + {\vec{E}}_{1 (2,3)}^{-} (\vec{r}, ω_{1 (2,3)}) = \sum_{m n} {\vec{E}}_{1 (2,3) m n} (z) e^{i (k_{1 (2,3) m n x} x + k_{1 (2,3) m n y} y - ω_{1 (2,3)} t)} = \sum_{m n} [{\vec{E}}_{1 (2,3) m n}^{+} (z) + {\vec{E}}_{1 (2,3) m n}^{-} (z)] e^{i (k_{1 (2,3) m n x} x + k_{1 (2,3) m n y} y - ω_{1 (2,3)} t)},$ (19)where the symbols $+ (-)$ shown in Eq. (19) represent the positive (negative) real parts of the $z$ components of the wave vector. If three frequencies are not equal, the nonlinear polarization along the $x$ direction on the graphene can be expressed as $P_{x} (ω_{3}) = ϵ_{0} χ^{(3)} E_{1 x} (\vec{r}, ω_{1}) E_{2 x}^{*} (\vec{r}, ω_{2}) E_{1 x} {(\vec{r}, ω_{1}) |}_{z = 0} = ϵ_{0} χ^{(3)} \sum_{otmnrs} {E_{1 o t x} (z) E_{2 m n x}^{*} (z) E_{1 r s x} (z) |}_{z = 0} \times e^{i (k_{3 x} x + G_{o - m + r x} x + k_{3 y} y + G_{t - n + s y} y - ω_{3} t)} .$ (20)

In the same way, the nonlinear polarization along the $y$ direction, $P_{y} (ω_{3})$ , can also be obtained. The directions of the transmitted and reflected FWM waves are determined by the in-plane phase-matching condition: $k_{3 x (y)} = 2 k_{1 x (y)} - k_{2 x (y)}$ . For example, we consider a particular case that all the frequencies (the pump, probe, and FWM waves) are equal, that is, $ω_{1} = ω_{2} = ω_{3} = ω_{0}$ . In such a case, if the pump wave is incident normally ( $θ_{1} = 0 °$ ), the angles of the probe and FWM waves are equal ( $θ_{2} = θ_{3} = θ$ ) due to the phase-matching condition, the transmitted and reflected FWM fields are optical phase conjugations. By considering the terms that have the same transverse momentum, we divide the third-order nonlinear polarization densities on the graphene along the $x$ direction into three parts, i.e., $P_{1 x} (ω_{0})$ , $P_{2 x} (ω_{0})$ , and $P_{3 x} (ω_{0})$ : $P_{1 x} (ω_{0}) = ϵ_{0} χ^{(3)} (3 E_{1 x} E_{1 x}^{*} E_{1 x} + 6 E_{2 x} E_{2 x}^{*} E_{1 x} + 6 E_{3 x} E_{3 x}^{*} E_{1 x} + 6 E_{2 x} E_{1 x}^{*} E_{3 x}) = ϵ_{0} χ^{(3)} \sum_{otmnrs} P_{1 x, otmnrs} e^{i (G_{o - m + n x} x + G_{t - n + s y} y - ω_{0} t)},$ (21a) $P_{2 x} (ω_{0}) = ϵ_{0} χ^{(3)} (3 E_{2 x} E_{2 x}^{*} E_{2 x} + 6 E_{3 x} E_{3 x}^{*} E_{2 x} + 6 E_{1 x} E_{1 x}^{*} E_{2 x} + 3 E_{1 x} E_{3 x}^{*} E_{1 x}) = ϵ_{0} χ^{(3)} \sum_{otmnrs} P_{2 x, otmnrs} e^{i [(k_{2 x} + G_{o - m + r t - n + s, x}) x + (k_{2 y} + G_{o - m + r t - n + s, y}) y - ω_{0} t]},$ (21b) $P_{3 x} (ω_{0}) = ϵ_{0} χ^{(3)} (3 E_{3 x} E_{3 x}^{*} E_{3 x} + 6 E_{1 x} E_{1 x}^{*} E_{3 x} + 6 E_{2 x} E_{2 x}^{*} E_{3 x} + 3 E_{1 x} E_{2 x}^{*} E_{1 x}) = ϵ_{0} χ^{(3)} \sum_{otmnrs} P_{3 x, otmnrs} e^{i [(- k_{2 x} + G_{o - m + r t - n + s x}) x + (- k_{2 y} + G_{o - m + r t - n + s y}) y - ω_{0} t]},$ (21c)with $P_{1 x, otmnrs} = 3 E_{1 o t x} E_{1 m n x}^{*} E_{1 r s x} + 6 E_{2 o t x} E_{2 m n x}^{*} E_{1 r s x} + 6 E_{3 o t x} E_{3 m n x}^{*} E_{1 r s x} + 6 E_{2 o t x} E_{1 m n x}^{*} E_{3 r s x},$ (22a) $P_{2 x, otmnrs} = 6 E_{1 o t x} E_{1 m n x}^{*} E_{2 r s x} + 3 E_{2 o t x} E_{2 m n x}^{*} E_{2 r s x} + 6 E_{3 o t x} E_{3 m n x}^{*} E_{2 r s x} + 3 E_{1 o t x} E_{3 m n x}^{*} E_{1 r s x},$ (22b) $P_{3 x, otmnrs} = 3 E_{1 o t x} E_{2 m n x}^{*} E_{1 r s x} + 6 E_{1 o t x} E_{1 m n x}^{*} E_{3 r s x} + 6 E_{2 o t x} E_{2 m n x}^{*} E_{3 r s x} + 3 E_{3 o t x} E_{3 m n x}^{*} E_{3 r s x} .$ (22c)

The corresponding $P_{1 y} (ω_{0})$ , $P_{2 y} (ω_{0})$ , and $P_{3 y} (ω_{0})$ can be obtained in the same way; all the electrical fields here are taken on the graphene. The first terms in Eq. (21) represent a phenomenon called self-phase modulation ^[61,62], which results from dipole excitations induced by three photons, the second and third terms correspond to cross-phase modulation, the fourth terms lead to the generation of the optical phase conjugations. After the third-order nonlinear polarization densities are calculated, the surface currents of the graphene, $J_{x}^{(3)} (ω_{3})$ and $J_{y}^{(3)} (ω_{3})$ , can be obtained. Similar to the method for the TH generation shown in previous part, the scattering matrixes in the present case can also be obtained; the field intensities for the transmitted ( $I_{3 t}$ ) and reflected ( $I_{3 r}$ ) FWM waves can be obtained through numerical calculations. Furthermore, we can calculate the transmission and reflection coefficients, $T$ and $R$ , to characterize the efficiencies of FWM generations: $T = I_{3 t} / I_{2}, T = I_{3 r} / I_{2} .$ (23)

Here, $I_{2}$ represent the intensities of the incident probe beams.

3. RESULTS AND DISCUSSIONS

In this section, we present numerical results for the TH generation and FWM from the three-layer structures. The dielectric constants of spheres are taken as 6.25, which correspond to thallium iodide ^[63]. The lattice constant is taken as $a = 600 nm$ . We take ${Al}_{0.14} {Ga}_{0.86} As$ ^[64] as the dielectric slab, the corresponding dielectric constant is taken as $ϵ_{d} (ω) = (6.3 + 19.0 x) [f (χ) + \frac{1}{2} {(\frac{E_{0}}{E_{0} + Δ_{0}})}^{\frac{3}{2}} f (χ_{s 0})] + 9.4 - 10.2 x$ ^[64]; here $χ = ℏ ω / E_{0}$ , $χ_{s 0} = ℏ ω / (E_{0} + Δ_{0})$ , $E_{0} = 1.425 + 1.155 x + 0.37 x^{2} (eV)$ , $E_{0} + Δ_{0} = 1.765 + 1.115 x + 0.37 x^{2} (eV)$ , and $x = 0.14$ . The radius of the sphere and the thickness of the dielectric slab are designed as $r_{s} = 0.48 a$ and $d = 0.76 a$ , respectively. The third-order nonlinear effects of the spheres and dielectric slab are neglected.

An important feature is that the BICs around the graphene layer are easy to appear in these three-layer structures, which has been discussed in the previous investigation ^[43]. In Fig. 2(a), we plot the electric field intensities at the interface (graphene layer) between the sphere and the slab for $E_{f} = 0.23 eV$ and $E_{f} = 0.7 eV$ . Here the P-polarized incident wave with the intensity $1 MW / {cm}^{2}$ is radiated to the $X Y$ plane of the sample with $θ = 5 °$ , P- or S-polarized wave represents the TM or TE mode. The field localizations are observed at different wavelengths due to the formation of bound states. These BICs lead to the appearance of the resonant absorption. Figure 2(b) displays the corresponding results of the absorption as a function of wavelength. It can be seen clearly that the resonant absorption reaches 67% at $λ = 1824.91 nm$ when the Fermi energy is taken as 0.23 eV, while the absorption peak is only 2.7% at $λ = 1824.51 nm$ when the Fermi energy is 0.7 eV. This is because the imaginary part of the conductivity of graphene is small (near zero) at $E_{f} = 0.7 eV$ , which can be seen from the inset of Fig. 2(a). From the inset in Fig. 2(a), we can also see that $E_{f} = 0.23 eV$ and $E_{f} = 0.7 eV$ correspond to two variation characteristics of the graphene conductivity. In the following, we study the generation signals for the TH and FWM waves based on these BICs.

Figure 2.(a) The absorption as a function of the wavelength $λ$ when the Fermi energy is taken as $E_{f} = 0.23 eV$ (black line) and $E_{f} = 0.7 eV$ (red line). Inset shows the conductivity of the graphene. (b) The electric field distribution at the graphene location below the bottom of the spheres when the Fermi energy is taken as $E_{f} = 0.23 eV$ (black line) and $E_{f} = 0.7 eV$ (red line).

A. Numerical Results for the TH Generation Based on BICs

In Fig. 3(a), we plot the TH conversion efficiency from the structure shown in Fig. 1(a) as a function of the FF wavelength $λ$ at $E_{f} = 0.23 eV$ . The corresponding results at $E_{f} = 0.7 eV$ are shown in Fig. 3(b). Here the intensity and incident angle of the pump P-polarized beam are taken as $1 MW / {cm}^{2}$ and 5°, which corresponds to the case in Fig. 2. For comparison, the TH conversion efficiencies of the monolayer graphene in vacuum are also plotted in Figs. 3(a) and 3(b). Because of the mirror symmetry of the freestanding monolayer graphene, the corresponding conversion efficiencies of F1 and F2 are nearly equal, so the black and red dashed lines coincide with each other.

Figure 3.Transmitted (F1), reflected (F2), and total (F3) TH conversions in the three-layer structure and freestanding graphene as a function of the FF wavelength $λ$ when the Fermi energy is taken as (a) $E_{f} = 0.23 eV$ and (b) $E_{f} = 0.7 eV$ . The TH enhancement G3 of TH signal in the three-layer structure compared with that in the freestanding graphene as a function of the FF wavelength $λ$ for (c) different radius $r_{s}$ of the spheres or (d) various thickness $d$ of the slab, the Fermi energy in (c) and (d) is taken as $E_{f} = 0.23 eV$ .

The total TH conversion efficiency F3 reaches $6.1 \times 10^{- 7}$ at the wavelength 1824.91 nm for $E_{f} = 0.23 eV$ and $4.7 \times 10^{- 5}$ at the wavelength 1824.51 nm for $E_{f} = 0.7 eV$ due to the improved field localizations from BICs. The corresponding TH signals are enhanced about 5 and 7 orders of magnitude at two cases, respectively, as compared with the bare monolayer graphene. The efficiency of the forward optical phase conjugate wave and backward optical phase conjugate wave is nearly equal ^[16]. The calculations in Ref. ^[28] have shown that the TH generation can be enhanced 5 orders of magnitude using the defect in 1D PC. However, for such a way, the sample with dozens of layers has to be fabricated. Instead, only three layers are needed in our structure.

In addition to the high efficiency, the signals generated in this way have good monochromaticity due to high $Q$ factors of bound states. The monochromaticity can be described by the relative linewidth $Δ λ / λ_{Γ}$ in the signal spectrum; here $λ_{Γ}$ is the center wavelength of the peak and $Δ λ$ represents full width at half maximum (FWHM) of the peak. From Figs. 3(a) and 3(b), we obtain that $Δ λ / λ_{Γ}$ are 0.0356% and 0.0175% for the cases with $E_{f} = 0.23 eV$ and $E_{f} = 0.7 eV$ , respectively, which is much less than the relative linewidth (about 1% in Ref. ^[65]). This means that the TH conversion efficiency is low when the frequency is taken at the deviation of the BICs. For example, at 1822 nm the TH conversion in the three-layer structure is nearly as large as that in the freestanding graphene; at the 1815 nm, the TH conversion in the three-layer structure is 2 orders of magnitude smaller than that in the freestanding graphene.

Figures 3(c) and 3(d) display the calculated results of the TH enhancement G3 as a function of the FF wavelength $λ$ for different radius $r_{s}$ of the spheres and various thickness $d$ of the slab. The other parameters are identical with those in Fig. 3(a). It is clearly seen from Fig. 3(c) that the peak of G3 is slightly redshifted with the increase of the radius of the spheres. In Fig. 3(d) the peak of G3 is redshifted obviously with the increase of the thickness $d$ of the slab. We can utilize this property to tune the peak location of the TH signal by changing the thickness $d$ of the slab.

The above discussions are only for a fixed incident angle of the incident wave. In fact, the phenomena strongly depend on the incident angle of the wave. Figure 4 shows the calculated results for the TH enhancement as functions of the FF wavelength and the incident angle at $E_{f} = 0.23 eV$ . Here, the intensity of the pump beam is taken as $1 MW / {cm}^{2}$ . Figures 4(a) and 4(b) correspond to the transmitted enhancement for the S- and P-polarized pump beams, respectively; the corresponding results for the reflected enhancement are given in Figs. 4(c) and 4(d). It is seen clearly that both G1 and G2 are sensitive to the incident angle. There are common TH enhancement peaks around the wavelength 1825.57 nm for S- and P-polarized pump beams located at the normal incident angle.

Figure 4.TH enhancement of the (a) and (b) transmitted and (c) and (d) reflected TH beams as a function of the fundamental wavelength and the incident angle for S-(P) polarized pump beam. (e) and (f) show the absorption spectra of the three-layer structure. The bound states are denoted by S1, S2, and S3 (P1 and P2) for S-(P) polarized wave.

With the increase of the incident angle from 0° to 30°, the peaks split into three for the S-polarized pump beam (S1, S2, and S3) and do not split for the P-polarized pump beam (P1). At the same time, we find that S2 and P1 overlap each other, which are blueshifted with the increase of the incident angle. These phenomena originated from the formation of BICs. This can be seen clearly in Figs. 4(e) and 4(f), which display the corresponding absorptions as functions of the wavelength and the incident angle of the incident wave for S- and P-polarized waves, respectively. The angle-dependent resonant absorptions are clearly visible in Figs. 4(e) and 4(f) corresponding to the existence of trapped EM modes (BICs) in the structure. The change characters of angle-dependent BICs are analogous to those of TH generations, as shown in Figs. 4(a)–4(d). This means that the formation of BICs leads to angle dependences of TH generations.

B. Numerical Results for the FWM Based on BICs

For the FWM from the structure shown in Fig. 1(b), we consider two kinds of cases: Case I and Case II. For Case I, the waves keep the same frequency in the process, as shown in Fig. 5(a), while the waves are in different frequencies for Case II. For Case I, the FWM from bare graphene has been discussed in Ref. ^[15]. In Fig. 5, we present the comparison between our results with those from the bare graphene; here, the wavelengths are denoted by $λ$ .

Figure 5.(a) Schematic of the FWM progress in the three-layer structure when all the waves are at the same wavelength; the incident angles of the pump and probe beams are 0° and $θ$ , respectively; the angle of the generated FWM wave is equal to $θ$ because of phase-match condition. (b) and (c) correspond to $T$ and $R$ of the transmitted and reflected FWM beams as a function of the FF wavelength and incident angle of the probe beam. The $T$ and $R$ in the three-layer structure and freestanding graphene as a function of the fundamental wavelength at $θ = 16.09 °$ corresponding to the maxima when the Fermi energy is taken as (d) $E_{f} = 0.23 eV$ and (e) $E_{f} = 0.7 eV$ . The blue and pink lines show the absorption spectra under only one incident beam at the incident angles $θ = 0 °$ and $θ = 16.09 °$ , respectively.

Figures 5(b) and 5(c) show the transmission ( $T$ ) and reflection ( $R$ ) of the generated FWM signals as functions of the wavelength $λ$ and the angle $θ$ of the probe beam at $E_{f} = 0.23 eV$ . Here the P-polarized waves are used as the probe and pump beams; the intensities of the pump and probe beams are taken as $1$ and $0.1 MW / {cm}^{2}$ . These intensities lead to a negligible saturation effect of the graphene. We can see that three resonant peaks (denoted by 1, 2, and 3) appear in Figs. 5(b) and 5(c), respectively; these phenomena are originated from the formation of BICs of the pump and probe beams in this structure. The resonant modes 1 and 2 correspond to the resonant modes P1 and P2, as shown in Fig. 4(f); the vertical resonant mode 3 corresponds to the BIC of the pump beam with the normal incident angle. When all beams are in BICs, the maxima of the generated FWM signals can be obtained, as shown by bright dots in Figs. 5(b) and 5(c). This can be seen more clearly in Fig. 5(d).

In Fig. 5(d), we plot the transmission (black lines) and the reflection (black triangle lines) as a function of the wavelength $λ$ . Here, the incident angle of the probe beam is fixed at $θ = 16.09 °$ ; the other parameters are identical to those in Figs. 5(b) and 5(c), where a resonant peak [corresponding to the red point in Fig. 5(b)] around $λ = 1825.57 nm$ appears. In order to understand this phenomenon, we also plot the corresponding absorption spectra; the blue line corresponds to the absorption when the beam irradiates the structure with the normal incident angle; the pink line corresponds to the absorption spectrum, where the incident angle is taken as 16.09°. We can see that two absorption peaks for two cases correspond to the resonant peaks for the generated FWM signals, which also correspond to two BICs for the pump and probe beams. When the pump and probe beams are located at the BICs, the corresponding EM fields are localized at the graphene layer, and then the corresponding third-order polarization densities are largely enhanced. This is the reason why we can observe enhancement of FWM signals about 6 orders of magnitude around the resonant peak at $λ = 1825.57 nm$ . The relative linewidth of the resonant peak is 0.0346%, which means that good monochromaticity of the generated signal is observed again.

After the Fermi energy is tuned to $E_{f} = 0.7 eV$ , while other parameters are not changed, the corresponding results are shown in Fig. 5(e). In such a case, the FWM signal is enhanced about 9 orders of magnitude around the resonant peak, and the relative linewidth of the resonant peak reaches 0.00863%. This is because of small imaginary part and low loss of absorption, which corresponds to strong localization of the electric field, as described in Fig. 2(a).

For the bare graphene, the previous investigation has shown that the optical phase conjugation can be tuned in a large range by changing the phase difference ^[16]. In fact, the present structure still possesses such a property. The solid lines and dashed lines in Fig. 6 represent the transmission $T$ and reflection $R$ of the optical phase conjugations as a function of the phase difference between the two parts of the pump beam (the pump beam is divided into two parts with the same intensity $0.5 MW / {cm}^{2}$ ). The intensity and incident angle of the probe beam are taken as $0.1 MW / {cm}^{2}$ and 16.09°, respectively; the other parameters taken here correspond to the maximum in Fig. 5(c). A period modulation in the generation of optical phase conjugations is clearly observed; a maximum is reached when the phase difference is equal to 0°. A minimum appears when the phase difference is equal to 180°; the period is 360°. This phenomenon is similar to that in Ref. ^[16]; that is, we can also tune the optical phase conjugation in a large range by changing the phase difference in the present case.

Figure 6.Transmission ( $T$ ) and reflection ( $R$ ) as a function of the phase difference between the two parts of the pump beam with the same intensity.

For Case II, as shown in Fig. 7(a), we consider a pump beam normal to the surface of the structure at the wavelength $λ_{1}$ and a probe beam with an incident angle $θ_{2}$ and the wavelength $λ_{2}$ . Due to the third-order nonlinear effect of the graphene, the transmitted and reflected FWM signals at the angle $θ_{3}$ and wavelength $λ_{3}$ are generated. In this process, the following relations are satisfied: ${\begin{cases} ω_{3} = 2 ω_{1} - ω_{2} \Rightarrow \frac{2 π c}{λ_{3}} = 2 \frac{2 π c}{λ_{1}} - \frac{2 π c}{λ_{2}} \\ k_{3 x} = 2 k_{1 x} - k_{2 x} \Rightarrow \frac{2 π}{λ_{3}} \sin θ_{3} = 2 \frac{2 π}{λ_{1}} \cdot 0 - \frac{2 π}{λ_{2}} \sin θ_{2} \end{cases} .$ (24)

Figure 7.(a) Schematic of the FWM progress in the three-layer structure when the wavelengths of the pump and probe beams ( $λ_{1}$ and $λ_{2}$ ) are not equal, the incident angles of the pump and probe beams are $θ_{1} = 0 °$ and $θ_{2}$ , respectively. (b) and (c) correspond to the transmitted and reflected FWM beams as a function of the wavelength $λ_{1}$ and incident angle $θ_{2}$ of the probe beam. The wavelength of the probe beam is taken as $λ_{2} = 1760.29 nm$ .

Figures 7(b) and 7(c) display $T$ and $R$ of the generated FWM signals as a function of the wavelength of the pump beam $λ_{1}$ and the angle of the probe beam $θ_{2}$ ; the wavelength of the probe beam is taken as $λ_{2} = 1760.29 nm$ . The intensities of the pump and probe beams are taken as $1$ and $0.1 MW / {cm}^{2}$ ; the Fermi energy is $E_{f} = 0.23 eV$ . It is seen clearly that there are two vertical resonant modes (at $λ_{1} = 1705.11 nm$ and $λ_{1} = 1825.57 nm$ ) in every figure, which correspond to two BICs of the normal pump beam.

When the angles of the probe beam are taken as $θ_{2} = 60 °$ and $θ_{2} = 7.2 °$ , the probe beam is also in BICs. In such a case, strong FWM signals at $θ_{3} = 54.4 °$ (Point A) and $θ_{3} = 6.8 °$ (Point B) are generated [see bright dots in Figs. 7(b) and 7(c)]; the corresponding wavelength for these signals is $λ_{3} = 1653.29 nm$ . Points A and B are highlighted because the probe beam, pump beam, and generated FWM signal at these points are in BICs simultaneously; in such a case, very high generated efficiency can be obtained. The generated efficiency at point A is even more prominent because of the better electric field confinement of the probe beam. These can be seen more clearly in Fig. 8.

Figure 8. $T$ and $R$ in the three-layer structure and freestanding graphene as a function of the wavelength $λ_{1}$ with $θ_{2} = 60 °$ , which correspond to a peak when the Fermi energy is taken as (a) $E_{f} = 0.23 eV$ and (c) $E_{f} = 0.7 eV$ , respectively. (b) Blue and pink lines show the absorption as a function of the wavelengths of the pump or FWM beams ( $λ_{1}$ or $λ_{3}$ ) under one incident beam at the incident angles 0° and $θ_{3}$ , respectively; the two relations are put together via phase match condition.

In Fig. 8(a), we plot $T$ and $R$ as a function of $λ_{1}$ at $θ_{2} = 60 °$ and $E_{f} = 0.23 eV$ , which corresponds to the maxima in Figs. 7(a) and 7(b). The other parameters are identical with those in Figs. 7(a) and 7(b). It is clearly seen that two peaks appear, which are located at $λ_{1} = 1705.11 nm$ and $λ_{1} = 1825.57 nm$ . In order to clearly analyze the origin, in Fig. 8(b) we show the corresponding absorption spectrum of the pump beam (see blue line); two resonant peaks in the absorption spectrum correspond to those shown in Fig. 8(a). Furthermore, we also show the absorption spectrum of the generated FWM signal as a function of $λ_{3}$ (see red line). The resonant absorption peak appears around $λ_{3} = 1653.29 nm$ , which corresponds to $λ_{1} = 1705.11 nm$ . This means that the pump, probe, and FWM beams are located at BICs simultaneously in such a case, which leads to the enhancement of the generated FMW signal about 6 orders of magnitude. This is in contrast with the other cases; for example, the FWM signal is enhanced by about 2 orders of magnitude at $λ_{1} = 1825.57 nm$ because the pump beam is only in the BIC. When the Fermi energy is tuned to $E_{f} = 0.7 eV$ , the similar phenomena are observed; the corresponding results are shown in Fig. 8(c). Due to low loss and strong field localization, in such a case the FWM signal is enhanced by about 9 orders of magnitude at $λ_{1} = 1705.11 nm$ compared with that in the freestanding monolayer graphene; the relative line width of the corresponding resonant peak also reaches 0.0195%. The generation efficiency without the enhancement of BICs of the pump and FWM beams is shown away from the peaks in Figs. 8(a) and 8(c); it can be clearly seen that, at $λ_{1} = 1650 nm$ , the generation efficiency in the three-layer structure is nearly as large as that in the freestanding graphene; at $λ_{1} = 1750 nm$ , the generation efficiency in the three-layer structure is about 2 orders of magnitude smaller than that in the freestanding graphene.

We would like to point out that the purpose of the present paper is to provide theoretical proof of improved third-order nonlinear effect in graphene based on BICs. However, it is more convincing if the experimental observation can be provided. In fact, it is difficult to control the periodic pattern of the nanospheres; we have to optimize the structural design further so that the samples can be fabricated experimentally. These will be our next tasks.

4. SUMMARY

In summary, we have designed sphere-graphene-slab structures to improve the third-order nonlinear interaction between the EM wave and the graphene. The scattering matrix theory to calculate third-order nonlinear effects in these structures has been developed. Improved TH generation and FWM signals from the graphene in these structures has been discussed. Due to the field localization of the BICs, we have found that the efficiency of TH generation from the monolayer graphene in these structures can be enhanced about 7 orders of magnitude compared with that in bare graphene. By designing structure parameters, we can make all beams (the pump beam, probe beam, and generated FWM signal) be BICs at the same time. In such a case, the efficiency of FWM in the monolayer graphene can be enhanced about 9 orders of magnitude. It is interesting that the generated TH and FWM signals possess very good monochromaticity due to the high $Q$ factors of BICs, which are also sensitive to the incident angle of the wave. Furthermore, we can also obtain selective generations of TH and FWM signals for two polarized waves (S and P waves) by taking suitable BICs. The phenomena can appear at any wavelength from the mid-infrared to far-infrared band, which are very beneficial for the design of optical devices.

References

[1] E. Hendry, P. J. Hale, J. Moger, A. K. Savchenko, S. A. Mikhailov. Coherent nonlinear optical response of graphene. Phys. Rev. Lett., 105, 097401(2010).

[2] R. Wu, Y. Zhang, S. Yan, F. Bian, W. Wang, X. Bai, X. Lu, J. Zhao, E. Wang. Purely coherent nonlinear optical response in solution dispersions of graphene sheets. Nano Lett., 11, 5159-5164(2011).

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

[4] X. Yao, A. Belyanin. Giant optical nonlinearity of graphene in a strong magnetic field. Phys. Rev. Lett., 108, 255503(2012).

[5] M. Tokman, X. Yao, A. Belyanin. Generation of entangled photons in graphene in a strong magnetic field. Phys. Rev. Lett., 110, 077404(2013).

[6] Q. Bao, H. Zhang, Y. Wang, Z. Ni, Y. Yan, Z. X. Shen, K. P. Loh, D. Y. Tang. Atomic-layer graphene as a saturable absorber for ultrafast pulsed lasers. Adv. Funct. Mater., 19, 3077-3083(2009).

[7] W. B. Cho, J. W. Kim, H. W. Lee, S. Bae, B. H. Hong, S. Y. Choi, I. H. Baek, K. Kim, D.-I. Yeom, F. Rotermund. High-quality, large-area monolayer graphene for efficient bulk laser mode-locking near 1:25 μm. Opt. Lett., 36, 4089-4091(2011).

[8] D. Popa, Z. Sun, F. Torrisi, T. Hasan, F. Wang, A. Ferrari. Sub 200 fs pulse generation from a graphene mode-locked fiber laser. Appl. Phys. Lett., 97, 203106(2010).

[9] Y. W. Song, S. Y. Jang, W. S. Han, M. K. Bae. Graphene mode-lockers for fiber lasers functioned with evanescent field interaction. Appl. Phys. Lett., 96, 051122(2010).

[10] W. Tan, C. Su, R. Knize, G. Xie, L. Li, D. Tang. Mode locking of ceramic Nd: yttrium aluminum garnet with graphene as a saturable absorber. Appl. Phys. Lett., 96, 031106(2010).

[11] J.-L. Xu, X.-L. Li, Y.-Z. Wu, X.-P. Hao, J.-L. He, K.-J. Yang. Graphene saturable absorber mirror for ultra-fast-pulse solid-state laser. Opt. Lett., 36, 1948-1950(2011).

[12] J. Wang, Y. Hernandez, M. Lotya, J. N. Coleman, W. J. Blau. Broadband nonlinear optical response of graphene dispersions. Adv. Mater., 21, 2430-2435(2009).

[13] M. Feng, H. Zhan, Y. Chen. Nonlinear optical and optical limiting properties of graphene families. Appl. Phys. Lett., 96, 033107(2010).

[14] H. Yang, X. Feng, Q. Wang, H. Huang, W. Chen, A. T. Wee, W. Ji. Giant two-photon absorption in bilayer graphene. Nano Lett., 11, 2622-2627(2011).

[15] H. Harutyunyan, R. Beams, L. Novotny. Controllable optical negative refraction and phase conjugation in graphite thin films. Nat. Phys., 9, 423-425(2013).

[16] S. M. Rao, A. Lyons, T. Roger, M. Clerici, N. I. Zheludev, D. Faccio. Geometries for the coherent control of four-wave mixing in graphene multilayers. Sci. Rep., 5, 15399(2015).

[17] J. L. Cheng, N. Vermeulen, J. E. Sipe. Third order optical nonlinearity of graphene. New J. Phys., 16, 053014(2014).

[18] S.-Y. Hong, J. I. Dadap, N. Petrone, P.-C. Yeh, J. Hone, R. M. Osgood. Optical third-harmonic generation in graphene. Phys. Rev. X, 3, 021014(2013).

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

[20] S. A. Mikhailov. Quantum theory of third-harmonic generation in graphene. Phys. Rev. B, 90, 241301(2014).

[21] S. A. Mikhailov. Quantum theory of third-harmonic generation in graphene. Phys. Rev. B (erratum), 91, 039904(2014).

[22] J. L. Cheng, N. Vermeulen, J. E. Sipe. Third-order nonlinearity of graphene: effects of phenomenological relaxation and finite temperature. Phys. Rev. B, 91, 235320(2016).

[23] J. L. Cheng, N. Vermeulen, J. E. Sipe. Third-order nonlinearity of graphene: effects of phenomenological relaxation and finite temperature. Phys. Rev. B (erratum), 93, 039904(2016).

[24] A. V. Gorbach, E. Ivanov. Perturbation theory for graphene-integrated waveguide: cubic nonlinearity and third-harmonic generation. Phys. Rev. A, 94, 013811(2016).

[25] D. B. S. Soh, R. Hamerly, H. Mabuchi. Comprehensive analysis of the optical Kerr coefficient of graphene. Phys. Rev. A, 94, 023845(2016).

[26] Y. F. Song, L. Li, H. Zhang, D. Y. Shen, D. Y. Tang, K. P. Loh. Vector multi-soliton operation and interaction in a graphene mode-locked fiber laser. Opt. Express, 21, 10010-10018(2013).

[27] Z.-C. Luo, M. Liu, Z.-N. Guo, X.-F. Jiang, A.-P. Luo, C.-J. Zhao, X.-F. Yu, W.-C. Xu, H. Zhang. Microfiber-based few layer black phosphorus saturable absorber for ultra-fast fiber laser. Opt. Express, 23, 20030-20039(2015).

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

[29] J. Niu, M. Luo, Q. H. Liu. Enhancement of graphene’s third-harmonic generation with localized surface plasmon resonance under optical/electro-optic Kerr effects. J. Opt. Soc. Am. B, 33, 615-621(2016).

[30] A. V. Gorbach. Nonlinear graphene plasmonics: amplitude equation for surface plasmons. Phys. Rev. A, 87, 013830(2013).

[31] G. T. Adamashvili, D. J. Kaup. Optical surface breather in graphene. Phys. Rev. A, 95, 053801(2017).

[32] D. C. Marinica, A. G. Borisov, S. V. Shabanov. Bound states in the continuum in photonics. Phys. Rev. Lett., 100, 183902(2008).

[33] E. N. Bulgakov, A. F. Sadreev. Bound states in the continuum in photonic waveguides inspired by defects. Phys. Rev. B, 78, 075105(2008).

[34] Y. Plotnik, O. Peleg, F. Dreisow, M. Heinrich, S. Nolte, A. Szameit, M. Segev. Experimental observation of optical bound states in the continuum. Phys. Rev. Lett., 107, 183901(2011).

[35] M. I. Molina, A. E. Miroshnichenko, Y. S. Kivshar. Surface bound states in the continuum. Phys. Rev. Lett., 108, 070401(2012).