Comparison of Linear and Nonlinear Properties of Graphene and Silicene in Terahertz Range

Volodymyr Grimalsky¹, Jesus Escobedo-Alatorre¹, Yuriy Rapoport², Anatoliy Kotsarenko³

¹Center for Investigations on Engineering and Applied Science (CIICAp), Institute for Investigations on Basic and Applied Science (IICBA), Autonomous University of State Morelos (UAEM), Cuernavaca, Mexico.
²Space Radio-Diagnostic Research Centre, University of Warmia and Mazury (UWM), Olsztyn, Poland.
³Faculty of Engineering, Autonomous University of Carmen (UNACAR), Ciudad del Carmen, Mexico.

DOI: 10.4236/jemaa.2025.171001                

Abstract

Resonant linear and nonlinear properties in terahertz range of 2D materials graphene and silicene placed into a bias magnetic field are investigated theoretically on the base of the quasi-classical kinetic theory. When the electromagnetic frequency is close to the cyclotron one, the linear conductivity increases two orders. Under the resonant frequencies nonlinearity becomes essential at low magnitudes of terahertz electric fields. In absence of a bias magnetic field the nonlinear dependences of the surface electric currents on terahertz electric field are practically the same simulated from kinetics and electron hydrodynamics with nonzero “kinetic” electron effective mass. Graphene possesses higher values of nonlinearity of the resonant conductivity, whereas in absence of a bias magnetic field, the electron nonlinearity is higher in silicene.

Keywords

Terahertz Range, Graphene, Silicene, Resonant Conductivity, Nonlinearity

1. Introduction

Low dimensional materials [1]-[3] can be used as the basic elements to create the controlling and nonlinear metamaterial media in terahertz (THz) range [4]-[6] 0.1 - 30 THz. Nowadays there exists a variety of two-dimensional (2D) materials like graphene, silicene, germanene etc. [1]-[3]. The layered metamaterial structures “dielectric – 2D material -dielectric…” and the periodic sets of graphene strips on dielectric substrates are perspective to realize nonlinear wave propagation and interaction in THz range, like frequency conversion and generation of higher harmonics [7] [8], self-action, switching [9]-[12], and solitons [13] [14]. In THz range under the temperatures T ≥ 20 K in graphene and silicene the nonlinearity is due to pseudo relativistic dispersion law [7] [10] [11] [14] [15] of electrons. It is important for applications that in these materials the dynamic nonlinearity in THz range possesses the reactive, or non-dissipative, character [7], and the collision frequencies are small there ν = 10¹² − 5·10¹² s⁻¹. The following inequality is valid: ν = ω, where ω ≥ 10¹³ s⁻¹ is the circular frequency of the electromagnetic (EM) fields in THz range.

It is interesting to compare the properties of 2D electron gas in graphene and silicene in THz range. The band structure in these materials is similar near the Dirac cone, but the Fermi velocity is essentially smaller in the silicene [1] [2], whereas usually the collision frequency is smaller in the graphene [2].

The application of external magnetic fields to the structures with the solid-state plasmas changes qualitatively the wave propagation and yields additional possibilities to realize linear and nonlinear resonant phenomena. Due to the electron band structure, in graphene and silicone the cyclotron frequencies are in THz range ω_B = 10¹³ − 10¹⁴ s⁻¹ when the bias magnetic field is about 0.5 - 2 T [12] [15]. In 2D materials the cyclotron frequencies depend not only on the bias magnetic field but also on 2D electron concentrations [12].

There are different methods to analyze the electron properties of 2D materials in THz range. They are the direct quantum kinetic approach, including the Kubo one [14], the quasi-classical kinetics based on the Boltzmann equation [7] [8] [11], and the electron hydrodynamics [8]. At the temperatures T ≥ 20 K the simulations on the base of the Boltzmann equation are adequate to investigate the conductivity in 2D materials [7] [8].

In this paper, the resonant linear and nonlinear properties of 2D electron gas in graphene and silicene are investigated. Based on the Boltzmann kinetic approach, the linear and nonlinear expressions for the 2D tensor conductivity are derived. The quasi-linear method is used to analyze nonlinear resonant conductivity when the bias magnetic field is applied. For realistic values of 2D electron concentrations and collision frequencies, the values of the resonant conductivity coincide with ones obtained from the simpler hydrodynamic approach. In addition, the nonlinear dependences of the surface density of the electric current on THz electric field practically coincide when obtained both from the kinetics and from the hydrodynamics in the absence of a bias magnetic field. This makes possible to use the simpler hydrodynamic method to analyze the nonlinear propagation of THz EM waves through the layered structures “dielectric-graphene or silicene—dielectric”.

2. Basic Equations

The consideration is based on the Boltzmann kinetic equation [7] [8] [11] [12] [15] for the electron distribution function f, SI units:

$\frac{\partial f}{\partial t}+v\cdot \frac{\partial f}{\partial r}+e\left(E+v\times B\right)\cdot \frac{\partial f}{\partial p}=St\left\{f\right\}.$ (1)

Here p is the electron quasi-momentum, St{f} is the collision integral that below is considered within the collision frequency ν approximation $St\left\{f\right\}=-\nu \left(f-f_{00}\right)$ , f₀₀ is the equilibrium Fermi distribution function. The bias magnetic field B₀ is directed along OZ axis. This equation is valid in THz range, where the inequality takes place: ν < ω; ω ≥ 5⋅10¹² s⁻¹ is the frequency of EM field, ν = 10¹² − 5·10¹² s⁻¹ is the electron collision frequency. In THz range it is possible to neglect by interband electron transitions when the 2D electron concentration is n₂₀ > 10¹⁰ cm⁻². The effects due to the quantizing magnetic field are not considered here, so the bias magnetic field is not high: B₀ < 2 T.

Below the case of the two-component THz electric field is considered in the XOY plane: E = (E_x, E_y, 0). The kinetic Equation (1) is rewritten as:

$\frac{\partial f}{\partial t}+e\left(E_x\frac{\partial f}{\partial p_x}+E_y\frac{\partial f}{\partial p_y}\right)+e{B}_0\left(v_y\frac{\partial f}{\partial p_x}-v_x\frac{\partial f}{\partial p_x}\right)=-\nu \cdot \left(f-f_{00}\right).$ (2)

The stationary solution of Equation (1) is the Fermi distribution function:

$f_{00}\left(p\right)={\left(1+\exp \left(\frac{E\left(p\right)-E_F}{E_T}\right)\right)}^{-1}.$ (3)

Here $E_F$ is the Fermi energy in equilibrium, $E\left(p\right)$ is the electron energy, $E_T\equiv k_{B}T$ , T is a temperature in Kelvins. Our goal is to derive the dependence of the surface density of the electric current on THz electric field.

The following notations are used:

$p_F=\frac{E_F}{v_F},p_T=\frac{E_T}{v_F},m^{*}\equiv \frac{p_F}{v_F}>m_g^{*},{\omega }_B\equiv \frac{e{B}_0}{m^{*}}.$ (4)

Here p_F is the Fermi momentum, $m^{*}$ is so-called “kinetic” effective mass, which is utilized below, ω_B is the electron cyclotron frequency.

The equation for p_F is:

$n_{20}=\frac{4}{4{\pi }^2{\hslash }^2}{\displaystyle \underset{-\infty }{\overset{+\infty }{\int }}{\displaystyle \underset{-\infty }{\overset{+\infty }{\int }}{f}_{00}\left(p\right){\text{d}}^{2}p}}$ . (5)

Here, n₂₀ = 10¹⁰ – 10¹³ cm⁻² is the steady 2D concentration of the electron gas.

The dispersion law for 2D electrons in the graphene and silicene is [1] [2] [7]:

$\begin{array}{l}E=v_F{\left(p^2+p_g^2\right)}^{1/2},p^2\equiv p_x^2+p_y^2;\\ p_g\equiv m_g^{*}{v}_F.\end{array}$ (6)

Here, v_F is so-called the Fermi velocity, it is v_F ≈ 10⁶ m/s for graphene and v_F ≈ 4.5⋅10⁵ m/s for silicene [2]. For the sake of generality, the small effective mass $m_g^{*}$ ≤ 0.001m_e is introduced in Equation (6). The nonzero effective mass $m_g^{*}$ can occur in 2D materials placed onto some substrates, like SiO₂.

The electron velocity for the quasi-classical electron motion is obtained as:

$\begin{array}{l}v\equiv \frac{\partial E}{\partial p};v_x\equiv \frac{\partial E}{\partial p_x}=\frac{v_F{p}_x}{{\left(p^2+p_g^2\right)}^{1/2}},\\ v_y=\frac{v_F{p}_y}{{\left(p^2+p_g^2\right)}^{1/2}},p^2\equiv p_x^2+p_y^2;\\ p_x=p\cdot \cos \phi ,p_y=p\cdot \sin \phi .\end{array}$ (7)

Below the quasi-linear approach is applied [16] to derive the nonlinear conductivity. The solution of Equation (2) is searched in the form:

$\begin{array}{l}f=f_0\left(t,p\right)+\frac{1}{2}\left(\tilde{f}\left(p,\phi \right)\cdot \exp \left(-i\omega t\right)+c.c.\right);\\ E_{x,y}=\frac{1}{2}\left({\tilde{E}}_{x,y}\exp \left(-i\omega t\right)+c.c.\right).\end{array}$ (8)

Here f₀ is the almost constant, or basic, part of the distribution function. Another part of it oscillates with the frequency ω. In the quasi-linear approach, the inverse action of the high-frequency oscillations to the basic part of the wave function f₀ is taken into account. The generation of higher harmonics is not considered because it is essential under satisfying matching conditions, which are realized in a specified geometry [8].

From Equation (2) there is the following equation for $\tilde{f}$ :

$-\frac{e{B}_0{v}_F}{{\left(p^2+p_g^2\right)}^{1/2}}\frac{\partial \tilde{f}}{\partial \phi }+\nu \tilde{f}=-\frac{e}{p}\left(E_x{p}_x+E_y{p}_y\right)\frac{\partial f_0\left(p\right)}{\partial p}.$ (9)

$\begin{array}{l}E\cdot p=p\cdot \left({\tilde{E}}_{+}\exp \left(-i\phi \right)+{\tilde{E}}_{-}\exp \left(i\phi \right)\right),\\ {\tilde{E}}_{+,-}=\frac{1}{2}\left({\tilde{E}}_x\pm i{\tilde{E}}_y\right).\end{array}$

The equation for the perturbation of the distribution function can be rewritten as:

$\begin{array}{l}\left(-i\omega +\nu -\frac{{\omega }_B\cdot {\left(p_F^2+p_g^2\right)}^{1/2}}{{\left(p^2+p_g^2\right)}^{1/2}}\frac{\partial }{\partial \phi }\right)\tilde{f}\\ =-\frac{e{p}_{\perp }}{p}\left({\tilde{E}}_{+}\cdot \exp \left(-i\phi \right)+{\tilde{E}}_{-}\cdot \exp \left(i\phi \right)\right)\frac{\partial f_0\left(p\right)}{\partial p}.\end{array}$ (10)

The expression for the surface density of the electric current is [7]:

$i_{sx}=\frac{4e}{4{\pi }^2{\hslash }^2}{\displaystyle \underset{-\infty }{\overset{+\infty }{\int }}{\displaystyle \underset{-\infty }{\overset{+\infty }{\int }}{v}_{x}f\left(p\right){\text{d}}^{2}p}}.$ (11)

Then in Equation (11) the solution of Equation (10) in 2D case is used, and the formula for the surface density of current is:

$\begin{array}{c}{i}_{sx}=-\frac{i{e}^2{v}_F}{\pi {\hslash }^2}\cdot \{{\tilde{E}}_{+}{\displaystyle \underset{0}{\overset{+\infty }{\int }}\frac{p^2}{{\left(p^2+p_g^2\right)}^{1/2}\left(\omega +i\nu \right)-{\omega }_B{\left(p_F^2+p_g^2\right)}^{1/2}}\cdot \frac{\partial f_0}{\partial p}\text{d}p}\\ +{\tilde{E}}_{-}{\displaystyle \underset{0}{\overset{+\infty }{\int }}\frac{p^2}{{\left(p^2+p_g^2\right)}^{1/2}\left(\omega +i\nu \right)+{\omega }_B{\left(p_F^2+p_g^2\right)}^{1/2}}\cdot }\frac{\partial f_0}{\partial p}\text{d}p\}.\end{array}$ (12)

The analogous expression is for i_sy.

It is possible to release in the surface current density the resonant part i_s+ and the non-resonant one i_s- for EM fields of different circular polarizations:

$\begin{array}{l}{i}_{sx}\equiv {\sigma }_{1s}{E}_x+{\sigma }_{Bs}{E}_y\equiv i_{s+}+i_{s-};\\ i_{sy}\equiv -{\sigma }_{Bs}{E}_x+{\sigma }_{1s}{E}_y\equiv \frac{1}{i}\left(i_{s+}-i_{s-}\right).\end{array}$ (13)

The quasi-linear equation for f₀ is:

$\frac{\partial f_0}{\partial t}+\frac{e}{4}\langle \left({\tilde{E}}_x\frac{\partial {\tilde{f}}^{*}}{\partial p_x}+{\tilde{E}}_y\frac{\partial {\tilde{f}}^{*}}{\partial p_y}+c.c.\right)\rangle +\nu \cdot \left(f_0-f_{00}\right)=0.$ (14)

Then the following relation is used that results from relation for ${\tilde{E}}_{+}$ , Equation (9), the last line:

$\langle \left({\tilde{E}}_x\frac{\partial {\tilde{f}}^{*}}{\partial p_x}+{\tilde{E}}_y\frac{\partial {\tilde{f}}^{*}}{\partial p_y}+c.c.\right)\rangle =\frac{1}{p}\left[{\tilde{E}}_{+}\cdot \frac{\partial }{\partial p}\left(p\cdot {\tilde{f}}^{*}\right)+c.c.\right]$ . (15)

Here $\langle \dots \rangle$ is the symbol of averaging over the oscillation period 2π/ω.

As a result, the equation for f₀ is:

$\frac{\partial f_0}{\partial t}+\frac{e}{4p}\left[ie{\tilde{E}}_{+}\frac{\partial }{\partial p}\left(p\frac{\partial f_0}{\partial p}\frac{{\tilde{E}}_{+}{}^{*}}{\omega -i\nu -{\omega }_B\cdot Q^{1/2}}\right)+c.c.\right]+\nu \cdot \left(f_0-f_{00}\right)=0.$ (16)

The following notation is used:

$Q\equiv \frac{p_F^2+p_g^2}{p^2+p_g^2}.$ (17)

Thus, the equation for f₀ is:

$\frac{\partial f_0}{\partial t}-\frac{e{\left|{\tilde{E}}_{+}\right|}^2}{2p}\frac{\partial }{\partial p}\left(p\frac{\partial f_0}{\partial p}\frac{\nu }{{\left(\omega -{\omega }_B{Q}^{1/2}\right)}^2+{\nu }^2}\right)+\nu \cdot \left(f_0-f_{00}\right)=0.$ (18)

This equation possesses the diffusive character in the space of quasi-momenta [16]. The boundary conditions for f₀ are [16]:

$f_0\left(p\to +\infty \right)\to 0,\frac{\partial f_0}{\partial p}\left(p=0\right)=0.$ (19)

The expressions for the resonant and non-resonant surface conductivities are:

${\sigma }_{s+}\equiv {\sigma }_{1s}-i{\sigma }_{Bs};{{\displaystyle }}^{}{\sigma }_{s-}\equiv {\sigma }_{1s}+i{\sigma }_{Bs}.$ (20)

From Equation (12), the following expression for the resonant surface conductivity is written down:

${\sigma }_{s+}=-\frac{i{e}^2{v}_F}{\pi {\hslash }^2}{\displaystyle \underset{0}{\overset{+\infty }{\int }}\frac{p^2\frac{\partial f_0}{\partial p}}{\left(\omega +i\nu -{\omega }_B\cdot Q^{1/2}\right){\left(p^2+p_g^2\right)}^{1/2}}\text{d}p}.$ (21)

In the limiting case T → 0 there are simplified expressions from the kinetic theory for the linear 2D conductivity components [12] that coincide with ones obtained from the hydrodynamic approximation with the “kinetic” effective mass m^*:

$\begin{array}{l}{\sigma }_{1s}=\frac{i{e}^2{n}_{20}\cdot \left(\omega +i\nu \right)}{{\left(\omega +i\nu \right)}^2-{\omega }_B^2};{\sigma }_{Bs}=-\frac{e^2{n}_{20}{\omega }_B}{{\left(\omega +i\nu \right)}^2-{\omega }_B^2};\\ {\sigma }_{s+}=\frac{i{e}^2{n}_{20}}{\omega -{\omega }_B+i\nu }.\end{array}$ (22)

The nonzero “kinetic” effective mass m^* used in Equations (22) is taken from Equation (6).

In the hydrodynamic approach the surface density of current is expressed through the hydrodynamic velocity v [12]:

$\begin{array}{l}{i}_{sx,y}=e{n}_{20}{v}_{x,y};v_{+,-}\equiv \frac{1}{2}\left(v_x\pm i{v}_y\right);\\ i_{s+}=e{n}_{20}{v}_{+}\equiv {\sigma }_{s+}{\tilde{E}}_{+}.\end{array}$ (23)

In the kinetic approach the expression for σ_s+ can be written down by means of complex integrals, Equation (21). It possesses the resonant dependence on frequency ω. In the hydrodynamic approach the linear formula for σ_s+ is simple and is obtained from the equation for the electron velocity:

$\begin{array}{l}\left(-i\omega +\nu \right)v_{+}+i{\omega }_B{v}_{+}\approx \frac{e}{m^{*}}{\tilde{E}}_{+};\\ \left(-i\omega +\nu \right)v_{-}-i{\omega }_B{v}_{-}\approx \frac{e}{m^{*}}{\tilde{E}}_{-}.\end{array}$ (24)

It is seen that EM wave with the circular polarization ${\tilde{E}}_{+}$ is subject to the resonant interaction with 2D electron gas; another EM wave with the non-resonant circular polarization ${\tilde{E}}_{-}$ is not considered here.

Therefore, the hydrodynamic expression for the linear resonant conductivity is

${\sigma }_{s+}\approx \frac{i{e}^2{n}_{20}}{m^{*}}\frac{1}{\omega -{\omega }_B+i\nu }.$ (25)

Here ω_B ≡ eB₀/m^* is the cyclotron frequency for a particle with the “kinetic” effective mass m^*, also see Equation 4.

3. Simulations

3.1. Linear Resonant Conductivity

In this Subsection, the linear conductivity of graphene and silicene has been simulated. The dependences of the resonant component of the linear surface conductivity σ_s+ on frequency ω for graphene and silicene have been calculated with using Equations (21), (25) for the kinetics and the hydrodynamics correspondingly. In the linear case, the distribution function f₀ is equilibrium one, f₀ = f₀₀, Equation (3).

The typical dependences are presented in Figure 1 and Figure 2. The bias magnetic field is B₀ = 1 T, the electron collision frequency is ν = 2⋅10¹² s⁻¹. The results do not change when $m_g^{\ast }$ ≤ 0.001 m_e, see Equation (6). Parts a) are the dependences for the graphene, b) are ones for the silicene.

In Figure 1, 2D electron concentration is n₂₀ = 5⋅10¹⁰ cm⁻². The curves 1, 2, 3 are at T = 20, 30, and 50 K simulated from the kinetic approach. The curves 4, 5 are from the hydrodynamic approach T = 0, the collision frequency is the same for curve 4 ν = 2⋅10¹² s⁻¹, but it is increased ν = 4.8⋅10¹² s⁻¹ for curve 5. It is seen that the resonant frequencies are ω ≈ 3.93⋅10¹³ s⁻¹ for the graphene and ω ≈ 1.79⋅10¹³ s⁻¹ for the silicene there. They coincide with the cyclotron frequency ω_B, Equations (4), (25).

Figure 1. Dependences of the linear components of the resonant surface tensor conductivity σ_s+, namely the real (σ_s+', solid lines) and imaginary (σ_s+'', dot lines) parts, on frequency ω for graphene and silicene in the bias magnetic field B₀ = 1 T. The surface electron concentration is n₁₀ = 5⋅10¹⁰ cm⁻², the collision frequency is ν = 2⋅10¹² s⁻¹, except for the curve 5. Curves 1, 2, 3 are from kinetics, ones 4, 5 are from the hydrodynamics. Curves 1, 2, 3 are at the temperatures T = 20 K, 30 K, and 50 K, correspondingly. Curve 5 is for the increased collision frequency ν = 4.8⋅10¹² s⁻¹. Part (a) is for graphene, (b) is for silicene.

In Figure 2, 2D electron concentration is higher, n₂₀ = 3⋅10¹¹ cm⁻². The curves 1, 2, 3 are at T = 30, 50, 80 K, the kinetic approach. The resonant frequencies are ω ≡ ω_B ≈ 1.62⋅10¹³ s⁻¹ for the graphene and ω ≈ 0.73⋅10¹³ s⁻¹for the silicene. There is a good coincidence with curves 1 and 5, namely the last one is from the hydrodynamic approach with the slightly increased collision frequency ν = 2.5⋅10¹² s⁻¹.

Figure 2. The same as in Figure 1, but the surface electron concentration is higher, n₁₀ = 3⋅10¹¹ cm⁻², and the temperatures ate T = 30 K, 50 K, and 80 K. Curve 5 is from the hydrodynamic approach, but for the increased collision frequency ν = 2.5⋅10¹² s⁻¹.

The results of our simulations have shown that both in graphene and in silicene the dependences of the components of the linear resonant conductivity simulated from the hydrodynamic approach are of good agreement with ones obtained from the kinetics, when 2D electron concentrations are n₂₀ ≥ 3·10¹¹ cm⁻² and the collision frequencies are ν ≥ 5·10¹¹ s⁻¹. The absolute values of the resonant components of conductivity increase ≥100 times near the cyclotron frequency ω_B, see Equation (25). Because the Fermi velocity v_F in silicene is smaller than one in graphene, the resonant frequencies are smaller in the silicene, when the 2D electron concentration n₂₀ and the bias magnetic field B₀ are the same.

Our simulations have demonstrated that at smaller values of 2D electron concentrations n₂₀ < 1.5⋅10¹¹ cm⁻² there is an essential discrepancy between the values of the resonant surface conductivity calculated from kinetic and hydrodynamic approaches. But a comparison of curves 1 and 5 in Figure 1 shows that they practically coincide when the hydrodynamic approach is used with the increased value of the collision frequency.

3.2. Nonlinear Resonant Conductivity

In this Subsection, the quasi-linear approach is used to simulate the nonlinear dependences of the resonant 2D conductivity. The resulting equation for the basic distribution function f₀ is one with the nonlinear diffusion, Equations (18), (19). The following qualitative result is mentioned. Due to the nonlinearity, i.e. an influence of THz electric field $\left|{\tilde{E}}_{+}\right|$ on the basic distribution function f₀, the smoothing of f₀ occurs, which is equivalent to increasing the effective “kinetic” mass m^*.

In Figure 3 and Figure 4, there are the dependences of the basic distribution function f₀ on the quasi-momentum p at different values of THz electric field $\left|{\tilde{E}}_{+}\right|$ . The parameters correspond to ones used in Figure 1 and Figure 2; the chosen THz frequencies are equal to the resonant ones in Figure 1 and Figure 2, namely to the cyclotron frequency ω_B. The quasi-momentum p is normalized to p_T, see Equation (4). In the nonlinear cases the stationary distribution function f₀ is presented at the time moment t = 10⁻¹¹ s ? ν^-1.

In Figure 3 and Figure 4, curve 1 is the equilibrium distribution function f₀₀. Curves 2, 3, 4 are for the values of THz electric field $\left|{\tilde{E}}_{+}\right|$ = 10⁴ V/m, 3⋅10⁴ V/m, and 10⁵ V/m, correspondingly.

An influence of THz electric field $\left|{\tilde{E}}_{+}\right|$ results in the difference of the basic distribution function f₀ from the equilibrium one f₀₀ under relatively small values of $\left|{\tilde{E}}_{+}\right|$ ≥ 10⁴ V/m and, thus, in essential electron nonlinearity under the cyclotron resonance condition ω ≈ ω_B.

Dependences of the resonant surface conductivity on THz electric field are presented in Figure 5 and Figure 6. Curves 1 - 4 correspond to ones in Figure 3 and Figure 4. Note that curves 1 in Figure 1 and Figure 2 coincide with curves 1 in Figure 5 and Figure 6; it is the same linear case.

From Figure 5 and Figure 6, it is seen that under the cyclotron resonance conditions the electron nonlinearity is higher in graphene under equal 2D electron concentrations n₂₀. The nonlinearity is more essential at smaller values of n₂₀, compare Figure 5 and Figure 6.

Figure 3. Dependences of the basic component of the distribution function f₀ on the quasi-momentum for graphene and silicene in the magnetic field B₀ = 1 T. The surface electron concentration is n₁₀ = 5⋅10¹⁰ cm⁻², the THz frequency is ω= 3.93⋅10¹³ s⁻¹ for graphene and 1.79⋅10¹³ s⁻¹ for silicene. The collision frequency is ν = 2⋅10¹² s⁻¹. The initial electron temperature is T = 20 K. Curve 1 is the equilibrium distribution function f₀₀, 2 is f₀ at the THz fields $\left|{\tilde{E}}_{+}\right|$ = 10⁴ V/m, 3 is at $\left|{\tilde{E}}_{+}\right|$ = 3⋅10⁴ V/m, 4 is at $\left|{\tilde{E}}_{+}\right|$ = 10⁵ V/m. Part (a) is for graphene, (b) is for silicene.

Figure 4. The same as in Figure 3, but n₂₀ = 3⋅10¹¹ cm⁻², T = 30 K. The frequencies are 1.62⋅10¹³ s⁻¹ for graphene and 0.73⋅10¹³ s⁻¹ for silicene.

Figure 5. Nonlinear dependences of the resonant surface conductivity σ_s+. The real parts are given in solid lines, the imaginary ones in dot lines. The surface electron concentration is n₁₀ = 5⋅10¹⁰ cm⁻², the THz frequency is ω= 3.93⋅10¹³ s⁻¹ for graphene and 1.79⋅10¹³ s⁻¹ for silicene. Curve 1 is for the linear case; 2, ,3, 4 are for the values of THz electric field 10⁴ V/m, 3⋅10⁴ V/m, and 10⁵ V/m. Part (a) is for graphene, (b) is for silicene.

Figure 6. The same as in Figure 5, but n₂₀ = 3⋅10¹¹ cm⁻², T = 30 K. The frequencies are 1.62⋅10¹³ s⁻¹ for graphene and 0.73⋅10¹³ s⁻¹ for silicene.

3.3. Fully Nonlinear Approach without Bias Magnetic Field

Here, the nonlinear properties of graphene and silicene are compared under non-resonant conditions in the absence of a bias magnetic field B₀ = 0. THz EM field possesses the linear polarization E_x. The values of the frequency ω of THz electric field correspond to the resonant cases presented in Figure 1 and Figure 2.

The nonlinear properties of graphene and silicene have been investigated by means of both the kinetic approach [7] and the quantum one [14]. In the collisionless approximation the solution of the kinetic equation can be obtained by the method of characteristics [7]:

$\begin{array}{l}f\left(p\right)=f_0\left(p-p_A\right);p_A\equiv eA,E=\frac{\partial A}{\partial t};\\ f_0\left(p\right)={\left(1+\exp \left(\frac{E\left(p\right)-E_F}{E_T}\right)\right)}^{-1}.\end{array}$ (26)

Then the nonlinear dependence of the surface current density has been simulated from Equation (11) by means of direct integration [11].

To compare the results of the kinetics, also the direct quantum approach has been applied. The following expressions for the surface density of current in graphene were derived [14] from the fully quantum consideration at T = 0, i.e. when 2D electron gas is fully degenerated:

$i_{sx}\approx \frac{e{v}_F}{\pi }{\left(\frac{p_F}{\hslash }\right)}^2\tanh \Psi ;\Psi \equiv \frac{eA}{p_F}.$ (27)

It is possible to approximate this dependence as [14]:

$i_{sx}\approx \frac{e{v}_F}{\pi }{\left(\frac{p_F}{\hslash }\right)}^2\frac{\Psi }{{\left(1+{\Psi }^2\right)}^{1/2}}.$ (28)

Equation (28) also can be derived from the nonlinear hydrodynamic equation for the quasi-particles with the “kinetic” effective mass m^*, Equation (4), and the pseudo-relativistic dispersion law where the role of the velocity of light takes the Fermi velocity v_F [11] [14].

The nonlinear dependences of the surface current density are presented in Figure 7. These dependences are given for the resonant frequencies ω ≈ ω_B where the maximum values of the real parts of the linear surface conductivity have been obtained in the presence of the bias magnetic field given in Figure 2. In Figure 7, 2D electron concentration is n₂₀ = 3⋅10¹¹ cm⁻². The curves 1, 2, 3 are simulated from the kinetic approach at the electron temperatures T = 30, 50, and 80 K. Curve 4 is obtained from the quantum approach, Equation (27), curve 5 is from the simplified Equation (28).

Figure 7. Nonlinear dependences of the density of the surface electric current on THz electric field E_x for graphene, part a), and silicene, part b), in the absence of a bias magnetic field. The THz frequencies are ω = 1.62⋅10¹³ s⁻¹, part a), and ω = 0.73⋅10¹³ s⁻¹, part b). 2D electron concentration is n₀ = 3⋅10¹¹ cm⁻². Curves 1, 2, 3 are obtained from the kinetic approach under the temperatures T = 30 K, 50 K, and 80 K; curve 4 is from the quantum approach, Equation (27); curve 5 is from the quantum approach, simplified formula, Equation (28).

The nonlinear dependences of the surface current densities on THz electric field are sharper in the silicene than in the graphene. Therefore, in the absence of a bias magnetic field the electron nonlinearity can be higher in silicene. But the non-resonant nonlinearity manifests at the magnitudes of THz electric fields of about 10⁶ V/m, 2 orders higher than for the resonant nonlinearity in the bias magnetic fields.

Thus, the simpler hydrodynamic approach can be used to investigate the resonant nonlinear propagation of THz EM waves in the layered structures “dielectric-graphene (silicene)-dielectric…”, when the values of 2D electron concentrations are not small n₂₀ ≥ 3·10¹¹ cm⁻². Our simulations have shown that it is valid at the electron temperatures T ≤ 120 K.

4. Discussion

Our simulations have demonstrated that in the bias magnetic field both graphene and silicene possess a strong resonant dependence of the surface conductivity on frequency and essential nonlinearity at low magnitudes of THz electric fields of about 100 V/cm. As a nonlinear material, graphene is preferential there. The non-resonant nonlinearity in the absence of a bias magnetic field is more expressed in silicene, due to lower values of the Fermi velocity. Nonlinearity cannot be considered moderate in all cases.

Now photonic crystals, metamaterials, and other layered structures, which are designed for nonlinear applications in THz range, use 2D materials like graphene and silicene. It is of interest to combine the resonant nonlinearity due to a bias magnetic field with the geometrical resonances due to the thicknesses of dielectric layers. Various approaches can be used to analyze these complex structures. Due to the complexity, simple but adequate methods are needed to study essentially nonlinear phenomena there. The most suitable approach is the quasi-classical kinetic one, but also the application of simpler methods like electron hydrodynamics can be useful.

5. Conclusion

Expressions of resonant linear dependences of the surface conductivity in graphene and silicene at terahertz frequencies have been simulated from the kinetic approach, where the finite temperatures of 2D electron gas are considered. These dependences practically coincide with those obtained from a simpler hydrodynamic approach with the nonzero electron “kinetic” effective mass. This “kinetic” mass can be obtained from the direct quantum approach. At the resonant frequencies the essential nonlinearity occurs both in graphene and silicene at low magnitudes of terahertz electric fields of circular polarization, but the resonant nonlinearity in graphene is higher. In absence of a bias magnetic field nonlinear dependences of surface current densities on the terahertz electric field can also be obtained from the hydrodynamic approach, where the pseudo-relativistic dependence of the “kinetic” effective mass on the electron velocity is used. The Fermi velocity plays the role of the speed of light. The non-resonant nonlinearity is higher in silicene.

Acknowledgements

Y.R. acknowledges the National Science Centre, Poland, for funding initiative under project number UMO-2023/49/B/ST10/03465. V.G., J.E.-A., and A.K. are grateful to SEP-CONAHCyT (Mexico) for partial funding.

Conflicts of Interest

The authors declare no conflicts of interest regarding the publication of this paper.

References

[1]	Cahangirov, S., Sahin, H., Le Lay, G. and Rubio, A. (2017) Introduction to the Physics of Silicene and Other 2D Materials. Springer.
[2]	Spencer, M.J.S. and Morishita, T. (2016) Silicene, Structure, Properties and Applications. Springer.
[3]	Trivedi, S., Srivastava, A. and Kurchania, R. (2014) Silicene and Germanene: A First Principle Study of Electronic Structure and Effect of Hydrogenation-Passivation. Journal of Computational and Theoretical Nanoscience, 11, 781-788. https://doi.org/10.1166/jctn.2014.3428
[4]	Perenzoni, M. and Paul, D.J. (2014) Physics and Applications of Terahertz Radiation. Springer.
[5]	Rieh, J.-S. (2021) Introduction to Terahertz Electronics. Springer. https://doi.org/10. 1007/978-3-030-51842-4
[6]	Biswas, A., Banerjee, A., Acharyya, A., Inokawa, H. and Roy, J.N. (2020) Emerging Trends in Terahertz Solid-State Physics and Devices Sources, Detectors, Advanced Materials, and Light-Matter Interactions. Springer. https://doi.org/10.1007/978-981-15-3235-1
[7]	Mikhailov, S.A. (2007) Non-Linear Electromagnetic Response of Graphene. Europhysics Letters (EPL), 79, Article 27002. https://doi.org/10.1209/0295-5075/79/27002
[8]	Rapoport, Y., Grimalsky, V., Lavrinenko, A.V. and Boardman, A. (2017) Double Resonant Excitation of the Second Harmonic of Terahertz Raditation in Dielectric-Graphene Layered Metamaterials. Journal of Optics, 19, Article 095104. https://doi.org/10.1088/2040-8986/aa7f54
[9]	Kaipa, C.S.R., Yakovlev, A.B., Hanson, G.W., Padooru, Y.R., Medina, F. and Mesa, F. (2012) Enhanced Transmission with a Graphene-Dielectric Microstructure at Low-Terahertz Frequencies. Physical Review B, 85, 245407-245413. https://doi.org/10.1103/physrevb.85.245407
[10]	Rapoport, Y., Grimalsky, V., Iorsh, I., Kalinich, N., Koshevaya, S., Castrejon-Martinez, C., et al. (2013) Nonlinear Reshaping of Terahertz Pulses with Graphene Metamaterials. JETP Letters, 98, 503-506. https://doi.org/10.1134/s002136401321011x
[11]	Grimalsky, V., Koshevaya, S., Escobedo-Alatorre, J. and Rapoport, Y. (2017) Nonlinear Properties of Electron Gas in n-InSb and Graphene in THz Range under Finite Temperatures. 2017 IEEE 37th International Conference on Electronics and Nanotechnology (ELNANO), Kyiv, 18-20 April 2017, 37-41. https://doi.org/10.1109/elnano.2017.7939815
[12]	Grimalsky, V., Koshevaya, S., Rapoport, Y., Tretiak, N., Yanovsky, F. and Escobedo-Alatorre, J. (2019) Resonant Properties of Electron Gas in n-InSb and Graphene Layers in Magnetic Fields for THz Multilayered Dielectric-Plasma-Like Metamaterials. 2019 IEEE 39th International Conference on Electronics and Nanotechnology (ELNANO), Kyiv, 16-18 April 2019, 164-168. https://doi.org/10.1109/elnano.2019.8783772
[13]	Nesterov, M.L., Bravo-Abad, J., Nikitin, A.Y., García-Vidal, F.J. and Martin-Moreno, L. (2013) Graphene Supports the Propagation of Subwavelength Optical Solitons. Laser & Photonics Reviews, 7, L7-L11. https://doi.org/10.1002/lpor.201200079
[14]	Dong, H., Conti, C. and Biancalana, F. (2011) Terahertz Relativistic Spatial Solitons in Doped Graphene Metamaterials. arXiv:1107.5803 [physics. optics].
[15]	Grimalsky, V., Escobedo-Alatorre, J., Rapoport, Y. and Kotsarenko, A. (2024) Resonant Linear and Nonlinear Properties of Graphene and Silicene in Terahertz Range. 2024 25th International Microwave and Radar Conference (MIKON), Wroclaw, 1-4 July 2024, 345-348. https://doi.org/10.23919/mikon60251.2024.10633984
[16]	Sitenko, A. and Malnev, V. (1994) Plasma Physics Theory. Chapman & Hall.