Numerical modeling of nonlinear graphene-based devices at terahertz frequencies

Abstract

Rigorous mathematical models of nonlinear terahertz (THz) graphene-based devices are developed by solving nonlinear Maxwell’s equations with boundary conditions simultaneously with a model of graphene surface conductivity as a nonlinear function on the time-varying electric field. The numerical approach is based on the projection method. Nonlinear phenomena, related to the parametric instability due excitations by the incident pumping TEM-wave in the multilayer graphene-dielectric nanostructures, are investigated numerically at THz frequencies.

Keywords

Nonlinear parametric graphene-based devices nanostructure bias electric field numerical approach projection method

1. Introduction

Nonlinear effects have attracted attention in graphene, where a number of phenomena, including second and third harmonic generation, frequency mixing spatial, self-phase modulation and optical Kerr effect, have been theoretically discussed and recently were verified experimentally [1]. Four-wave mixing is an important parametric process that has numerous applications, such as parametric amplification, phase conjugation, wavelength conversion, optical sampling, signal regeneration [2].

In order to model graphene-comprising waveguides [3,4] and graphene-cavity [5] in the linear regime the commercial Finite-Element Method (FEM) and Finite-Difference Time-Domain Method (FDTDM) are utilized as a numerical tool. Nonlinear propagation in 3D nanophotonic graphene-comprising waveguides of several types is tackled by introducing the nonlinear polarization and the current as a perturbation to the linear solution [3]. The four-wave mixing effect in the hybrid graphene sheets waveguides [4] and graphene-silicon slow-light photonic crystal waveguides [6] was investigated by solving the modified coupled mode equations.

In contrast to these simplified nonlinear coupled-mode theory simulations [3–6] the alternative approach is the solution of 3D nonlinear boundary problems without any simplification of the Maxwell’s equations (Eqs.) and boundary conditions. In this work a numerical approach for the analysis of nonlinear, parametric interactions of electromagnetic waves with 3D graphene nanostructures in arbitrary configurations was developed to solve the nonlinear Maxwell’s Eqs. rigorously in order to recognize and predict nonlinear physical phenomena and effects aimed to design reconfigurable nonlinear graphene-based THz devices taking into account the constrained geometries. We have created a computing algorithm, using the numerical methods of autonomous blocs with the Floquet channels [9] for modeling and CAD of nonlinear graphene-based THz devices.

2. Rigorous mathematical model

The rigorous mathematical models of nonlinear graphene-contained THz devices are based on the solution of the nonlinear diffraction boundary problems for Maxwell’s Eqs.: $\begin{eqnarray}\displaystyle curl\vec{H}(t)={\varepsilon}_{0}{\varepsilon}_{b}\frac{\partial \vec{E}(t)}{\partial t}+{\sigma}(|\overrightarrow{E}(t)|)\vec{E}(t),\quad curl\vec{E}(t)=-{\mu}_{0}{\mu}_{b}\frac{\partial \overrightarrow{H}(t)}{\partial t}, & & \displaystyle\end{eqnarray}$ (1) where E (t) and H (t) are the electric and magnetic field intensity vectors; 𝜇_b, ϵ_b are a permittivity and a permeability of graphene; σ(|E (t)|) is a nonlinear graphene conductivity, σ = σ_s∕d, here σ_s is a graphene surface conductivity, d is a thickness of graphene sheet (we assume d = 1 nm), with electrodynamic boundary conditions, simultaneously with a model of graphene surface conductivity σ_s as a nonlinear function of the time variant bias electric field determined from the Kubo formula [7].

The nonlinear dependence σ_s(| E (t)|) is described by using the uniform approximation by the high (six) order even polynomial: $\begin{eqnarray}\displaystyle {\sigma}_{s}(|\overrightarrow{E}(t)|)={\sigma}_{0}+{\sigma}_{2}|\overrightarrow{E}(t)|^{2}+{\sigma}_{4}|\overrightarrow{E}(t)|^{4}+{\sigma}_{6}|\overrightarrow{E}(t)|^{6}, & & \displaystyle\end{eqnarray}$ (2) where the coefficients of the approximation are defined from Kubo formula [7].

Let’s represent the fields E (t), H (t) in the form of series of terms of the combination frequencies ω_m: $\begin{eqnarray}\displaystyle \overrightarrow{E}(t)=\mathop{\sum }_{m=-\infty }^{\infty }\overrightarrow{E}({\omega}_{m})\exp (i{\omega}_{m}t),\quad \overrightarrow{H}(t)=\mathop{\sum }_{m=-\infty }^{\infty }\overrightarrow{H}({\omega}_{m})\exp (i{\omega}_{m}t). & & \displaystyle\end{eqnarray}$ (3) Introducing (2) in (1) and taking into account (3), we obtain the system of the stationary nonlinear Maxwell’s Eqs. at the combination frequencies ω_m.: $\begin{eqnarray}curl\overrightarrow{H}({\omega}_{m})=i{\omega}_{m}{\varepsilon}_{0}\dot{{\varepsilon}}_{b}({\omega}_{m})\overrightarrow{E}({\omega}_{m})+\overrightarrow{J}({\omega}_{m}),\quad curl\overrightarrow{E}({\omega}_{m})=-i{\omega}_{m}{\mu}_{0}{\mu}_{b}({\omega}_{m})\overrightarrow{H}({\omega}_{m}),\end{eqnarray}$ (4) where $\begin{eqnarray}\displaystyle \dot{{\varepsilon}}_{b}={\varepsilon}_{b}-i\frac{{\sigma}_{0}}{{\varepsilon}_{0}{\omega}_{m}d}, & & \displaystyle \nonumber\\ \displaystyle \vec{J}_{b}({\omega}_{m})=\mathop{\sum }_{q=1}^{\infty }A_{q}\vec{E}({\omega}_{q}), & & \displaystyle\end{eqnarray}$ (5) $\begin{eqnarray}\displaystyle A_{q} & = & \displaystyle \dot{{\sigma}}_{2}({\omega}_{m})\mathop{\sum }_{r_{1}=0}^{3}E_{z}({\omega}_{r_{1}})\mathop{\sum }_{r_{2}=0}^{3}E_{z}({\omega}_{r_{2}}){\gamma}_{2}\nonumber\\ \displaystyle \text{} & \text{} & \displaystyle +∼\dot{{\sigma}}_{4}({\omega}_{m})\mathop{\sum }_{r_{1}=0}^{3}E_{z}({\omega}_{r_{1}})\mathop{\sum }_{r_{2}=0}^{3}E_{z}({\omega}_{r_{2}})\mathop{\sum }_{r_{3}=0}^{3}E_{z}({\omega}_{r_{3}})\mathop{\sum }_{r_{4}=0}^{3}E_{z}({\omega}_{r_{4}}){\gamma}_{4}\nonumber\\ \displaystyle \text{} & \text{} & \displaystyle +∼\dot{{\sigma}}_{6}({\omega}_{m})\mathop{\sum }_{r_{1}=0}^{3}E_{z}({\omega}_{r_{1}})\mathop{\sum }_{r_{2}=0}^{3}E_{z}({\omega}_{r_{2}})\mathop{\sum }_{r_{3}=0}^{3}E_{z}({\omega}_{r_{3}})\mathop{\sum }_{r_{4}=0}^{3}E_{z}({\omega}_{r_{4}})\mathop{\sum }_{r_{5}=0}^{3}E_{z}({\omega}_{r_{5}})\mathop{\sum }_{r_{6}=0}^{3}E_{z}({\omega}_{r_{6}}){\gamma}_{6}.\nonumber\end{eqnarray}$ Here m is a summarize index determined at the ensemble $\{k,n\}$ , $\begin{eqnarray}\displaystyle {\gamma}_{i} & = & \displaystyle \left\{\begin{array}{@{}l@{}}0,∼\text{if}∼{\omega}_{r,q}+{\omega}_{0,l_{1}}+{\omega}_{0,l_{2}}+\ldots +{\omega}_{0,l_{i}}\neq {\omega}_{k,n}\\ 1,∼\text{if}∼{\omega}_{r,q}+{\omega}_{0,l_{1}}+{\omega}_{0,l_{2}}+\ldots +{\omega}_{0,l_{i}}={\omega}_{k,n},∼i=2,4,6;\\ \end{array}\right.\nonumber\\ \displaystyle \text{} & \text{} & \displaystyle r_{1}=\left\{0,l_{1}\right\},r_{2}=\left\{0,l_{2}\right\},\ldots r_{6}=\left\{0,l_{6}\right\}.\nonumber\end{eqnarray}$ The expression for a current (5) was obtained assuming that the external bias electric field E ₀ = E ₀ z ₀ and the time variant bias electric field intensity E _z(t) (applied normal to the surface of graphene sheets) are not dependent on the coordinates because of small thickness of a monoatomic graphene layer.

3. Numerical approach

We consider the nonlinear THz graphene-based devices as “waveguide transformers” (WT) [8] with the Floquet channels [9] on the inputs cross-sections S _𝛽. The descriptors of WT [8] are determined by solving the 3D diffraction boundary problems, taken into account the boundary conditions.

We find the solution of the 3D diffraction problem for the stationary Maxwell’s Eq. (4) in the form of Fourier’s series using eigenwaves $\{\boldsymbol{E}_{\boldsymbol{n}}\}$ , $\{\boldsymbol{H}_{\boldsymbol{n}}\}$ of the rectangular cavity (inside the basic region of WT): $\begin{eqnarray}\vec{E}({\omega}_{m})=\mathop{\sum }_{n=1}^{N_{bas}}\tilde{a}_{n}({\omega}_{m})\vec{E}_{n},\quad \vec{H}({\omega}_{m})=\mathop{\sum }_{n=1}^{N_{bas}}\tilde{b}_{n}({\omega}_{m})\vec{H}_{n},\quad m=1,2,\ldots M,\end{eqnarray}$ (6) where N _bas is a number of taken basis functions, M is a number of taken combinational frequencies, or eigenwaves $\{\boldsymbol{e}_{l({\beta})}\},\{\boldsymbol{h}_{l({\beta})}\}$ of Floquet channels [9] on the inputs cross-sections S _𝛽 of WT: $\begin{eqnarray}\vec{E}_{{\beta}}({\omega}_{m})=\mathop{\sum }_{l=1}^{L}\left(c_{l({\beta})}^{+}({\omega}_{m})+c_{l({\beta})}^{-}({\omega}_{m})\right)\vec{e}_{l({\beta})},\quad \vec{H}_{{\beta}}({\omega}_{m})=\mathop{\sum }_{l=1}^{L}\left(c_{l({\beta})}^{+}({\omega}_{m})-c_{l({\beta})}^{-}({\omega}_{m})\right)\vec{h}_{l({\beta})},\end{eqnarray}$ (7) where $C_{l({\beta})}^{+}({\omega}_{m})$ , $C_{l({\beta})}^{-}({\omega}_{m})$ are the magnitudes of incident and reflected modes, l are the indices of eigenwaves of Floquet channels, L is a number of taken modes in Floquet channel.

Taking into account (5), we represent the Maxwell’s Eqs. (4) in the projection integral form [10]: $\begin{eqnarray}\displaystyle \text{} & \text{} & \displaystyle \oint _{S_{{\Sigma}}}(\vec{E}({\omega}_{m})\times \vec{H}_{k}^{\ast })\cdot d\vec{S}=-i∼{\omega}∼{\mu}_{0}∼\int _{V_{0}}{\mu}∼\vec{H}({\omega}_{m})\cdot \vec{H}_{k}^{\ast }∼dV+i∼{\omega}_{k∼}{\varepsilon}_{0}∼\int _{V_{0}}\vec{E}({\omega}_{m})\cdot \vec{E}_{k}^{\ast }∼dV,\nonumber\\ \displaystyle \text{} & \text{} & \displaystyle \oint _{S_{{\Sigma}}}(\vec{H}({\omega}_{m})\times \vec{E}_{k}^{\ast })\cdot d\vec{S}=i∼{\omega}∼{\varepsilon}_{0}∼\int _{V_{0}}{\varepsilon}∼\vec{E}({\omega}_{m})\cdot \vec{E}_{k}^{\ast }∼dV-i∼{\omega}_{k}∼{\mu}_{0}\int _{V_{0}}\vec{H}({\omega}_{m})\cdot \vec{H}_{k}^{\ast }∼dV\nonumber\\ \displaystyle \text{} & \text{} & \displaystyle \quad +∼\mathop{\sum }_{q=1}^{\infty }A_{q}\int _{V_{0}}∼E({\omega}_{q})\cdot \vec{E}_{k}^{\ast }∼dV,\quad k=1,2,\ldots ,N_{\text{bas}},\quad m=1,2,\ldots ,M,\end{eqnarray}$ (8) were S _𝛴 = S ₁ ∪ S ₃, with the non-asymptotic radiation boundary conditions [11] on the inputs cross-sections S _𝛼 of WT: $\begin{eqnarray}\displaystyle \begin{array}{@{}l@{}}\displaystyle \frac{1}{2}\int _{S_{{\alpha}}}(\vec{E}_{{\alpha}}({\omega}_{m})\times ∼\vec{h}_{k({\alpha})}^{\ast })\cdot d\vec{S}_{{\alpha}}-\frac{1}{2}\int _{S_{{\alpha}}}(\vec{e}_{k({\alpha})}\times \vec{H}_{{\alpha}}^{\ast }({\omega}_{m}))\cdot d\vec{S}_{{\alpha}}=c_{k({\alpha})}^{-}({\omega}_{m}),\\ \quad k=1,2,\ldots ,N,\quad m=1,2,\ldots ,M.\end{array} & & \displaystyle\end{eqnarray}$ (9) Substituting the Fourier’s series (6), (7) into (8) we arrive to systems of matrix Eqs. at combinational frequencies ω_m: $\begin{eqnarray}\displaystyle \begin{array}{@{}l@{}}\displaystyle \mathbf{N}\cdot \mathbf{c}^{-}({\omega}_{m})-\mathbf{A}\cdot \tilde{\mathbf{a}}({\omega}_{m})+\mathbf{B}\cdot \tilde{\mathbf{b}}({\omega}_{m})=-\mathbf{N}\cdot \mathbf{c}^{+}({\omega}_{m}),\\ \displaystyle \mathbf{M}\cdot \mathbf{c}^{-}({\omega}_{m})+\mathbf{C}\cdot \tilde{\mathbf{a}}({\omega}_{m})-\mathbf{D}\cdot \tilde{\mathbf{b}}({\omega}_{m})+\mathop{\sum }_{q=1}^{L}A_{q}\mathbf{F}\cdot \mathbf{c}^{-}({\omega}_{q})=\mathbf{M}\cdot \mathbf{c}^{+}({\omega}_{m})-\mathop{\sum }_{q=1}^{L}A_{q}\mathbf{F}\cdot \mathbf{c}^{+}({\omega}_{q}),\\ \displaystyle \mathbf{I}\cdot \mathbf{c}^{-}({\omega}_{m})-\mathbf{R}\cdot \tilde{\mathbf{a}}({\omega}_{m})+\mathbf{U}\cdot \tilde{\mathbf{b}}({\omega}_{m})=\mathbf{0},\quad m=1,2,\ldots M,\end{array} & & \displaystyle\end{eqnarray}$ (10) where the components of a (ω_m), b (ω_m) and c ⁺(ω_m), c ⁻(ω_m) are the Fourier coefficients from (6), (7), correspondingly; I is a unity matrix, A , B , M , C , D , F , R , U are the matrices of vector functions [12].

When the magnitudes $C_{l({\beta})}^{+}({\omega}_{m})$ of the modes, incident on the input cross-sections S _𝛽 of WT, are known, the unknown magnitudes $C_{l({\beta})}^{-}({\omega}_{m})$ of reflected (in the local co-ordinate systems) modes of Floquet channels on the input cross-sections S _𝛼 of WT at combinational frequency ω_m can be found by solving the systems of linear algebraic Eqs. (SLAE) (10).

4. Results of numerical simulation of the parametric THz device based on graphene nanostructure

Let us consider a specific example for the investigation into parametric interactions of electromagnetic waves with graphene-based nanostructures at THz frequencies. Numerical calculations were performed for the case of the parametric THz device based on the multilayer graphene-dielectric nanostructure (Fig. 1), containing N graphene sheets, layers of dielectric SiO ₂ (ϵ_a = 2.2) of thickness d ₁ = 10 nm, and a layer of dielectric 𝛼 − Ge (ϵ_a = 18.5) of thickness d ₂ = 2800 nm. We consider the parametric THz device (Fig. 1) as a WT [8] with the Floquet channels [9] on the inputs cross-sections S _𝛽, 𝛽 = 1,2, …4. (The local coordinate systems are used on the cross-sections S _𝛽 of WT [8]).

Let’s assume that two electromagnetic waves: the signal wave of frequency ω₁ and the pumping wave of frequency ω₂ are incident on the input cross-sections S ₁ and S ₂ of WT, correspondingly (Fig. 1). There are TEM-modes, having magnitudes $C_{1(1)}^{+}({\omega}_{1})$ and $C_{1(2)}^{+}({\omega}_{2})$ .

Fig. 1.

Parametric graphene-based THz device as a “waveguide transformer” with the Floquet channels on the inputs cross-sections S ₁, S ₂, S ₃, S ₄: the signal TEM-wave (frequency ω₁, magnitude $C_{1(1)}^{+}({\omega}_{1})$ ), the pumping TEM-wave (frequency ω₂, magnitude $C_{1(2)}^{+}({\omega}_{2})$ ); transmitted (or excited) TEM-wave (frequency ω₁, magnitude $C_{1(3)}^{-}({\omega}_{1})$ .

The external bias electric field E ₀ = E ₀ z ₀ and the electric field E _z(ω₂) of the pumping TEM-wave should be applied normal to the surface of graphene sheets (Fig. 1). The electric field of pumping wave E _z(t) = E _zm cos ω₂ t (a large signal) acts as a time variant bias [13], modifying a chemical potential 𝜇_c. The graphene surface conductivity is a nonlinear function σ_s(𝜇_c(E _z0 + E _z(t)) versus the chemical potential 𝜇_c and, consequently, versus the time variant bias electric field E _z0 + E _z(t) = E _z0 + E _zm cos ω₂ t, expressed by Kubo formula [7].

Using the computational algorithm based on the decomposition approach by autonomous blocks with Floquet channels [9] the results of electrodynamic calculation of the modulus of scattering parameter |S ₃₁| of multilayer graphene-dielectric nanostructure (Fig. 1), depending on the frequency, for different values of chemical potential 𝜇_c (external bias electric field E _z0) were obtained and shown in Fig. 2 (dashed-line curve: N = 39; solid-line curve: N = 79).

Fig. 2.

Scattering parameter |S ₃₁| of multilayer graphene-dielectric nanostructure depending on the frequency for different values of chemical potential 𝜇_c (external bias electric field E _0z): curves1 −𝜇_c = 0 (E _z0 = 0); 2 −𝜇_c = 0.2 eV (E _z0 = 0.5 V/nm); 3 −𝜇_c = 0.4 eV (E _z0 = 1.0 V/nm); dashed-line curve N = 39; solid-line curve N = 79).

As shown in Fig. 2, upon increasing the value of chemical potential 𝜇_c the location of the maximum of the transmission coefficient |S ₃₁| through the multilayer graphene-dielectric nanostructure (Fig. 1) is moving and, consequently, the resonance frequency ω₀ of this band pass filter can be tuned significantly by changing the external bias electric field E _z0.

If the stable regime of the regenerative parametric amplification is analyzed (the magnitude of signal TEM-wave $C_{1(1)}^{+}({\omega}_{1})\neq 0$ ), then the system of Eqs. (11) is not uniform, and it has a solution at the SLAE’s determinant is nonzero. The results of computer analysis of the normalized magnitude $|C_{1(3)}^{-}({\omega}_{1})|/|C_{1(1)}^{+}({\omega}_{1})|$ , where $|C_{1(3)}^{-}({\omega}_{1})|$ is the magnitude of the transmitted TEM-mode on output cross-section S ₃ of parametric THz amplifier (Fig. 1), $C_{1(1)}^{+}({\omega}_{1})$ is the magnitude of signal TEM-wave of frequency f ₁ = ω₁∕2π = 14.15 THz, depending on the magnitude $C_{1(2)}^{+}({\omega}_{2})$ of pumping TEM-wave (f ₂ = ω₂∕2π = 28.3 THz) for 𝜇 = 0.2 eV (E _z0 = 0.5 V/nm) are shown in Fig. 3.

Fig. 3.

Normalized magnitude $|C_{1(3)}^{-}({\omega}_{1})|/|C_{1(1)}^{+}({\omega}_{1})|$ of transmitted TEM-wave on the output cross-section of parametric graphene-based THz amplifier with respect to the magnitude $|C_{1(1)}^{+}({\omega}_{1})|$ of signal TEM-wave, depending on the magnitude $C_{1(2)}^{+}({\omega}_{2})$ of pumping TEM-wave: f ₁ = 14.15 THz, f ₂ = 28.3 THz, N = 39; 𝜇_c = 0.2 eV (E _z0 =0.5 V/nm).

If the regime of the parametric generation is analyzed $(C_{1(1)}^{+}({\omega}_{1})=0)$ , then system of Eqs (10) is homogeneous, and it has nonzero solution at the zero value of the determinant of this SLAE. That’s why it is easy to see the physical meaning of numerical solutions. By computing the determinant of SLAE (10), the nonlinear regime of the onset of labile oscillating mode causing the wave instability due the parametric excitation by the incident pumping TEM-wave in the parametric THz graphene-based device (Fig. 1) was simulated. It is shown that the nonlinear phenomena, related to the parametric instability, appear for the case N ≥ 79 graphene sheets in the multilayer graphene-dielectric nanostructure (Fig. 1).

The results of analysis of instability regions, depending on the magnitude $|C_{1(2)}^{+}({\omega}_{2})|$ of pumping TEM-wave and the normalized frequency (ω₀∕ω₂)² with ω₀ is the resonance frequency of the multilayer (N = 79) graphene-dielectric nanostructure (f ₀ = ω₀∕2π = 14.15 THz (Fig. 2)), are shown in Fig. 4 for 𝜇_c = 0.2 eV (E _z0 = 0.5 V/nm).

Fig. 4.

Instability regions for parametric excitation in the graphene-based THz generator, depending on the magnitude $|C_{1(2)}^{+}({\omega}_{2})|$ of the pumping TEM-wave and the normalized frequency ω₀∕ω₂, were ω₀ is the resonance frequency of the multilayer graphene-dielectric nanostructure, f ₀ = 14.15 THz; N = 79; 𝜇_c = 0.2 eV (E _z0 = 0.5 V/nm).

In the regions of instability the determinant of SLAE (10) is equal to zero (or close to zero), consequently, the multilayer (N = 79) graphene-dielectric nanostructure (Fig. 1) becomes a parametric generator at frequencies ω₀ = nω₂∕2, where n =1,2,3, … (f ₀ = f ₂∕2 = 14.15 THz for the first order n = 1 and f ₀ = f ₂ = 28.3 THz for the second order n = 2 processes of the parametric excitation (Fig. 4)), tunable by changing the resonance frequency f ₀ of graphene nanostructure.

Using rigorous mathematical models taking into account the constrained geometries it is shown that the threshold magnitudes of parametric instability in the parametric THz graphene-based devices depend on the value of chemical potential, as a result, it can be tuned by the external bias electric field and also controlled by engineering the configuration and sizes of graphene nanostructures.

In contrast to previous works [3–6,14–16], the numerical approach, developed in this paper, permits to investigate the instabilities of fields of electromagnetic waves in strongly nonlinear bounded 3D graphene nanostructures; to determine the new solutions of nonlinear Maxwell’s equations; to follow the changes of the wave instabilities depending on the values of control parameters and their sensitivity to the transition regimes of nonlinear graphene-based THz devices. Using these rigorous mathematical models new results are expected to predict the nonlinear diffraction, multistability, generation of higher order time harmonics, solitons etc. in 3D photonics graphene-based devices.

Using this technique, reliable engineering methods for the CAD for the numerical computation of nonlinear electromagnetic properties of graphene-based nanomaterials, applied in the modern high technology, and future nonlinear photonic and optoelectronic graphene-based devices (frequency mixers, frequency multipliers, generators of higher order time harmonics) may be developed.

Footnotes

Acknowledgements

This work is carried out within the theme “Spin” No. 01.2.006 01201463330, with partial support of grant from Russian Science Foundation and No. 17-12-01002 (Section 4).

References

Glazov

M.M.

and Ganichev

S.D.

, High frequency electric field induced nonlinear effects in graphene, Physics Reports 535 (2014), 101.

Cai

Deng

Huang

Xia

and Wang

, Enhanced parametric frequency conversion in a compact silicon-graphene microring resonator, Optics Express 23 (2015), 18679–18685.

Chatzidimitriou

Pitilakis

and Kriezis

E.lE.

, Rigorous calculation of nonlinear parameters in graphene-comprising waveguides, J. Appl. Phys. 118(2) (2015), 023105.

Sun

Qiao

and Sun

, Direct generation of graphene plasmonic polaritons at THz frequencies via four wave mixing in the hybrid graphene sheets waveguides, Optics Express 22, 23 (2014), 22880–27891.

Petrone

McMillan

J.F.

van der Zande

G.Q.

Kwong

D.L.

Hone

and Wong

C.W.

, Regenerative oscillation and four-wave mixing in graphene optoelectronics, Nat. Photonics 6 (2012), 554–559.

Zhou

McMillan

J.F.

Petrone

van der Zande

Hone

J.C.

Kwong

D.L.

Feng

Zhou

and Wong

C.W.

, Enhanced four-wave mixing in graphene-silicon slow-light photonic crystal waveguides, Applied Physics Letters 105 (2014), 091111.

Hanson

G.W.

, Dyadic Green’s functions and guided surface waves for a surface conductivity model of graphene, J. Appl. Phys. 103 (2008), 064302.

Nikol’skii

V.V.

and Nikol’skaya

T. I.

, Decomposition Approach in Electrodynamics Problems, Nauka, Moscow, 1983, (in Russian).

Golovanov

O.A.

, Solution of applied problems of electromagnetics based on the use of autonomous blocks with virtual floquet channels, J. Commun. Techn. Electron. 51 (2006), 1423–1430.

10.

Nikolskii

V.V.

, Projection methods in electrodynamics, in: Collection of papers on applied electrodynamics, VisshayaShkola, Moscow, 1977, (in Russian).

11.

Golovanov

O.A.

and Makeeva

G.S.

, The method of autonomous blocks with magnetic nanoinclusions and floquet channels applied for simulation of magnetic nanostructures with allowance for the exchange and boundary conditions, J. Commun. Techn. Electron. 54 (2009), 1345.

12.

Golovanov

O.A.

, Numerical algorithm for solving the problems of diffraction for microwave waveguides with nonlinear media, J. Comm. Tech. and Electron. 35 (1990), 1853–1863.

13.

Mizuji

E.G.

Abdolali

Aghamohamadi

Danaeifar

Hashemi

and Tehrani

K.N.

, Realization of THz band mixer using graphene, Advanced Electromagnetics 3 (2014), 15–19.

14.

Petrone

McMillan

Zande

A.V.

Kwong

Hone

and Wong

C.W.

, Regenerative oscillation and four-wave mixing in graphene optoelectronics, Nature Photonics 6 (2012), 554–559.

15.

Yao

B.C.

Feng

Q.Y.

Cao

X.L.

Zhou

X.Y.

Rao

Y.J.

Gong

Zhang

W.L.

Wang

Z.G.

Chen

Y.F.

and Chiang

K.S.

, Generation of cascaded four-wave-mixing with graphene-coated microfiber, Photonics Research 3 (2015), A64–A68.

16.

Hoseini

and Malekmohammad

, All-optical high performance graphene-photonic crystal switch, Optics Communications 383 (2017), 159–164.