Theory of the strongly nonlinear electrodynamic response of graphene: A hot electron model

TL;DR

This paper develops a hot electron model to analyze the nonlinear electrodynamic response of graphene under strong electromagnetic radiation, providing detailed solutions and comparing with experimental data.

Contribution

It introduces a hot electron model for graphene's nonlinear response, deriving equations and solutions for various physical quantities under different conditions.

Findings

Absorption coefficient varies with radiation intensity at different frequencies.

Model predictions agree well with recent experimental data.

Identifies regimes of increasing and decreasing absorption with intensity.

Abstract

An electrodynamic response of graphene to a strong electromagnetic radiation is considered. A hot electron model (HEM) is introduced and a corresponding system of nonlinear equations is formulated. Solutions of this system are found and discussed in detail for intrinsic and doped graphene: the hot electron temperature, non-equilibrium electron and holes densities, absorption coefficient and other physical quantities are calculated as functions of the incident wave frequency $\omega$ and intensity $I$, of the equilibrium chemical potential $\mu_0$ and temperature $T_0$, scattering parameters, as well as of the ratio $\tau_\epsilon/\tau_{\rm rec}$ of the intra-band energy relaxation time $\tau_\epsilon$ to the recombination time $\tau_{\rm rec}$. The influence of the radiation intensity on the absorption coefficient $A$ at low ($\hbar\omega\lesssim 2|\mu_0|$, $dA/dI>0$) and high ($\hbar\omega\gtrsim 2|\mu_0|$, $dA/dI<0$) frequencies is studied. The results are shown to be in good agreement with recent experimental data.

Equations

$$E_{l\mathbf{k}}=(-1{)}^l\hslash v_F|\mathbf{k}|,$$

$$f_0(E,{\mu }_0,T_0)=\frac{1}{1+\exp \left(\frac{E-{\mu }_0}{T_0}\right)},$$

$$n_e^0=4.35\times 10^7{\text{ cm}}^{-2},n_h^0=3.12\times 10^{12}{\text{ cm}}^{-2},\text{(doped)}$$

$$n_e^0=n_h^0=8.14\times 10^{10}{\text{ cm}}^{-2},\text{(intrinsic)}.$$

$${\tau }_{ee}\ll {\tau }_p\ll {\tau }_{ϵ};$$

$$f(E,{\mu }_e,{\mu }_h,T)=\frac{\Theta (E)}{1+\exp \left(\frac{E-{\mu }_e}{T}\right)}+\frac{\Theta (-E)}{1+\exp \left(\frac{E-{\mu }_h}{T}\right)},$$

$$f_e(E,{\mu }_e,T)=\frac{1}{1+\exp \left(\frac{E-{\mu }_e}{T}\right)},E>0,$$

$$f_h(E,{\mu }_h,T)=\frac{1}{1+\exp \left(\frac{{\mu }_h-E}{T}\right)},E<0.$$

$$n_e({\mu }_e,T)=\frac{g_s{g}_v}{S}\sum _{\mathbf{k}}{f}_e(E_{2\mathbf{k}},{\mu }_e,T)=\frac{2T^2}{\pi (\hslash v_F{)}^2}{F}_1\left(\frac{{\mu }_e}{T}\right),$$

$$n_h({\mu }_h,T)=\frac{2T^2}{\pi (\hslash v_F{)}^2}{F}_1\left(-\frac{{\mu }_h}{T}\right),$$

$$F_n(z)={\int }_0^{\infty }\frac{x^{n}dx}{1+\exp \left(x-z\right)}.$$

$${\mathcal{E}}_e({\mu }_e,T)=\frac{g_s{g}_v}{S}\sum _{\mathbf{k}}{E}_{2\mathbf{k}}{f}_e(E_{2\mathbf{k}},{\mu }_e,T)=\frac{2T^3}{\pi (\hslash v_F{)}^2}{F}_2\left(\frac{{\mu }_e}{T}\right)$$

$${\mathcal{E}}_h({\mu }_h,T)=\frac{2T^3}{\pi (\hslash v_F{)}^2}{F}_2\left(-\frac{{\mu }_h}{T}\right).$$

$$\mathcal{E}({\mu }_e,{\mu }_h,T)={\mathcal{E}}_e({\mu }_e,T)+{\mathcal{E}}_h({\mu }_h,T)$$

$${\sigma }^{(1)}(\omega ,{\mu }_e,{\mu }_h,T)={\sigma }_{\mathrm{intra}}^{(1),e}(\omega ,{\mu }_e,T)+{\sigma }_{\mathrm{intra}}^{(1),h}(\omega ,{\mu }_h,T)+{\sigma }_{\mathrm{inter}}^{(1)}(\omega ,{\mu }_e,{\mu }_h,T),$$

$$\frac{{\sigma }_{\mathrm{intra}}^{(1),e}(\omega ,{\mu }_e,T)}{\frac{e^2{g}_s{g}_v}{16\hslash }}=\frac{i}{\pi \hslash T}{\int }_0^{\infty }\frac{EdE}{\omega +i{\gamma }_{\mathrm{intra}}(E)}\frac{1}{{\cosh }^2\left(\frac{E-{\mu }_e}{2T}\right)},$$

$$\frac{{\sigma }_{\mathrm{intra}}^{(1),h}(\omega ,{\mu }_h,T)}{\frac{e^2{g}_s{g}_v}{16\hslash }}=\frac{i}{\pi \hslash T}{\int }_0^{\infty }\frac{EdE}{\omega +i{\gamma }_{\mathrm{intra}}(E)}\frac{1}{{\cosh }^2\left(\frac{E+{\mu }_h}{2T}\right)},$$

$$\frac{{\sigma }_{\mathrm{inter}}^{(1)}(\omega ,{\mu }_e,{\mu }_h,T)}{\frac{e^2{g}_s{g}_v}{16\hslash }}=\frac{-i}{\pi }{\int }_0^{\infty }dE\left(\frac{1}{1+\exp \left(\frac{-E-{\mu }_h}{T}\right)}-\frac{1}{1+\exp \left(\frac{E-{\mu }_e}{T}\right)}\right)\frac{\hslash (\omega +i{\gamma }_{\mathrm{inter}})}{E^2-[\hslash (\omega +i{\gamma }_{\mathrm{inter}})/2{]}^2}.$$

$$\hslash {\gamma }_{\mathrm{intra}}(E)=\frac{|E|}{\frac{\zeta }{2}+\sqrt{1+\frac{E^4}{E_i^4}}-1}$$

$$\frac{{\sigma }_{intra}^{(1)}(\omega ,{\mu }_e,{\mu }_h,T)}{\frac{e^2{g}_s{g}_v}{16\hslash }}=\frac{4i}{\pi }\frac{T}{\hslash (\omega +i{\gamma }_p)}\left[F_0\left(\frac{{\mu }_e}{T}\right)+F_0\left(-\frac{{\mu }_h}{T}\right)\right]=\frac{4i}{\pi }\frac{T}{\hslash (\omega +i{\gamma }_p)}\ln \left[\left(1+e^{{\mu }_e/T}\right)\left(1+e^{-{\mu }_h/T}\right)\right].$$

$$\text{Re }\frac{{\sigma }_{\mathrm{inter}}^{(1)}(\omega ,{\mu }_e,{\mu }_h,T)}{\frac{e^2{g}_s{g}_v}{16\hslash }}=\frac{\sinh \left(\frac{\hslash |\omega |-({\mu }_e-{\mu }_h)}{2T}\right)}{\cosh \left(\frac{{\mu }_h+{\mu }_e}{2T}\right)+\cosh \left(\frac{\hslash |\omega |-({\mu }_e-{\mu }_h)}{2T}\right)}$$

$$\text{Im }\frac{{\sigma }_{\mathrm{inter}}^{(1)}(\omega ,{\mu }_e,{\mu }_h,T)}{\frac{e^2{g}_s{g}_v}{16\hslash }}=-\frac{\hslash \omega }{\pi }\mathcal{P}{\int }_0^{\infty }\frac{\sinh \left(\frac{E-({\mu }_e-{\mu }_h)/2}{T}\right)}{\cosh \left(\frac{{\mu }_h+{\mu }_e}{2T}\right)+\cosh \left(\frac{E-({\mu }_e-{\mu }_h)/2}{T}\right)}\frac{dE}{E^2-(\hslash \omega /2{)}^2}.$$

$$A(\omega ,{\mu }_e,{\mu }_h,T)=A^{\mathrm{intra}}(\omega ,{\mu }_e,{\mu }_h,T)+A^{\mathrm{inter}}(\omega ,{\mu }_e,{\mu }_h,T),$$

$$A^{\mathrm{intra}}(\omega ,{\mu }_e,{\mu }_h,T)=A_e^{\mathrm{intra}}(\omega ,{\mu }_e,T)+A_h^{\mathrm{intra}}(\omega ,{\mu }_h,T)=\frac{\frac{4\pi }{c}\text{Re }{\sigma }_{\mathrm{intra}}^{(1)}(\omega ,{\mu }_e,{\mu }_h,T)}{{\left|\frac{n+1}{2}+\frac{2\pi }{c}{\sigma }^{(1)}(\omega ,{\mu }_e,{\mu }_h,T)\right|}^2},$$

$$A^{\mathrm{inter}}(\omega ,{\mu }_e,{\mu }_h,T)=\frac{\frac{4\pi }{c}\text{Re }{\sigma }_{\mathrm{inter}}^{(1)}(\omega ,{\mu }_e,{\mu }_h,T)}{{\left|\frac{n+1}{2}+\frac{2\pi }{c}{\sigma }^{(1)}(\omega ,{\mu }_e,{\mu }_h,T)\right|}^2},$$

$$\frac{\partial {\mathcal{E}}_e({\mu }_e,T)}{\partial t}=A_e^{\mathrm{intra}}(\omega ,{\mu }_e,T)I-\frac{{\mathcal{E}}_e({\mu }_e,T)-{\mathcal{E}}_e({\mu }_e^0,T_0)}{{\tau }_{ϵ}},$$

$$\frac{\partial {\mathcal{E}}_h({\mu }_h,T)}{\partial t}=A_h^{\mathrm{intra}}(\omega ,{\mu }_h,T)I-\frac{{\mathcal{E}}_h({\mu }_h,T)-{\mathcal{E}}_h({\mu }_h^0,T_0)}{{\tau }_{ϵ}},$$

$$n_e({\mu }_e,T)=n_e({\mu }_e^0,T_0),n_h({\mu }_h,T)=n_h({\mu }_h^0,T_0).$$

$$\frac{\partial \mathcal{E}({\mu }_e,{\mu }_h,T)}{\partial t}=A^{\mathrm{intra}}(\omega ,{\mu }_e,{\mu }_h,T)I-\frac{\mathcal{E}({\mu }_e,{\mu }_h,T)-\mathcal{E}({\mu }_e^0,{\mu }_h^0,T_0)}{{\tau }_{ϵ}}.$$

$$n_e(t)=n_e^0+\delta n(t),n_h(t)=n_h^0+\delta n(t),$$

Full text

1. I Introduction

2. II Theory

3. II.1 The system in equilibrium

4. II.2 Hot electrons distribution functions

5. II.3 Electron and hole densities

6. II.4 Electron and hole energy densities

7. II.5 Conductivity

8. II.5.1 Intra-band conductivity

9. II.5.2 Inter-band dynamic conductivity

10. II.6 Absorption coefficient

11. II.7 Dynamics equations of the hot electron model

12. II.7.1 Energy relaxation

13. II.7.2 Recombination

14. II.7.3 Preliminary summary and discussion

15. II.8 Steady-state equations

16. III Results

17. III.1 Density-temperature diagrams

18. III.2 Power dependencies of different physical quantities

19. III.2.1 Doped graphene

20. III.2.2 Intrinsic graphene

21. III.3 Frequency dependencies of different physical quantities

22. IV Summary

23. A A model for the intra-band scattering rate

I Introduction

The nonlinear electrodynamic response of graphene attracted great attention in recent years. After the pioneering prediction Mikhailov (2007) of the strongly nonlinear electrodynamic properties of graphene, a large number of theoretical Mikhailov and Ziegler (2008); Dean and van Driel (2009, 2010); Smirnova et al. (2014); Savostianova and Mikhailov (2015); Wang et al. (2016); Rostami and Polini (2016); Cheng et al. (2017); Savostianova and Mikhailov (2017); Mikhailov (2012); Yao et al. (2014); Tokman et al. (2016); Cheng et al. (2014a); Marini et al. (2017); Savostianova and Mikhailov (2018); Mikhailov (2009a); Yao and Belyanin (2013); Tokman et al. (2014); Mikhailov (2011); Mikhailov and Beba (2012); Peres et al. (2014); Cox and de Abajo (2014, 2015); Cox et al. (2016); Mikhailov (2017a, b); Cheng et al. (2014b, 2015a); Mikhailov (2016, 2009b); Ishikawa (2010); Cheng et al. (2015b); Semnani et al. (2016); Mikhailov (2019) and experimental Dragoman et al. (2010); Bykov et al. (2012); Kumar et al. (2013); Hong et al. (2013); Soavi et al. (2018); Hendry et al. (2010); Gu et al. (2012); Alexander et al. (2017); König-Otto et al. (2017); Alexander et al. (2018); Bao et al. (2009); Zhang et al. (2009); Winnerl et al. (2011); Zheng et al. (2012); Bianchi et al. (2017); Zhang et al. (2012); Chen et al. (2013); Miao et al. (2015); Dremetsika et al. (2016); Vermeulen et al. (2016); Tomadin et al. (2018) papers have been published. Theoretically the higher harmonics generation Mikhailov (2007); Mikhailov and Ziegler (2008); Dean and van Driel (2009, 2010); Smirnova et al. (2014); Savostianova and Mikhailov (2015); Wang et al. (2016); Rostami and Polini (2016); Cheng et al. (2017); Savostianova and Mikhailov (2017), the frequency mixing Mikhailov (2012); Yao et al. (2014); Tokman et al. (2016), the direct current induced second harmonic generation Cheng et al. (2014a), the saturable absorption and Kerr effects Marini et al. (2017); Savostianova and Mikhailov (2018) have been studied. The nonlinear graphene response in magnetic fields Mikhailov (2009a); Yao and Belyanin (2013); Tokman et al. (2014) and the plasma wave related nonlinear effects Mikhailov (2011); Mikhailov and Beba (2012); Yao et al. (2014); Cox and de Abajo (2014, 2015); Cox et al. (2016); Mikhailov (2017a) have been also discussed in detail. A nonperturbative quasiclassical theory based on the relaxation time approximation and a quantum perturbation theory of all third order nonlinear effects have been developed in Ref. Mikhailov (2017b) and in Refs. Cheng et al. (2014b, 2015a); Mikhailov (2016) respectively. Experimentally the higher harmonics generation Dragoman et al. (2010); Bykov et al. (2012); Kumar et al. (2013); Hong et al. (2013); Soavi et al. (2018), the four-wave mixing Hendry et al. (2010); Gu et al. (2012); Alexander et al. (2017); König-Otto et al. (2017); Alexander et al. (2018), the radiation induced absorption changesBao et al. (2009); Zhang et al. (2009); Winnerl et al. (2011); Zheng et al. (2012); Bianchi et al. (2017), Kerr effect Zhang et al. (2012); Chen et al. (2013); Miao et al. (2015); Dremetsika et al. (2016); Vermeulen et al. (2016), the photoconductivity Tomadin et al. (2018) and other nonlinear phenomena have been observed. All of them demonstrated very large absolute values of the nonlinear optical parameters of graphene. The nonlinear electrodynamic properties of graphene can be used in many applications including broadband detection Ryzhii et al. (2012); Gan et al. (2013); Wang et al. (2013); Pospischil et al. (2013); Ryzhii et al. (2015), electrically tunable modulation of terahertz Liu et al. (2015); Kindness et al. (2018) and optical radiation Liu et al. (2011); Phare et al. (2015), mode-locked lasers Sun et al. (2010); Zhang et al. (2010); Popa et al. (2010, 2011); Bao et al. (2011) and other Ryzhii et al. (2007); Otsuji et al. (2012); Ryzhii et al. (2013).

The third-order fourth-rank conductivity tensor of graphene $\sigma_{\alpha\beta\gamma\delta}^{(3)}(\omega_{1},\omega_{2},\omega_{3};\mu_{0},T_{0})$, analytically calculated in the quantum theory Cheng et al. (2014b, 2015a); Mikhailov (2016), describes all possible third-order nonlinear effects for arbitrary polarizations and frequencies of the incident waves. Since $\sigma_{\alpha\beta\gamma\delta}^{(3)}$ was calculated within the perturbation theory, it depends, apart from the input frequencies $\omega_{1}$, $\omega_{2}$, $\omega_{3}$, on the equilibrium chemical potential $\mu_{0}$ and equilibrium temperature $T_{0}$. In many nonlinear response experiments, however, the incident radiation is so strong that the system gets excited far beyond the equilibrium state, and the use of parameters $\mu_{0}$, $T_{0}$ becomes not fully relevant. This problem can be partly circumvented by replacing $T_{0}$ in the expression for $\sigma_{\alpha\beta\gamma\delta}^{(3)}$ by an effective temperature $T$ which is considered as a fitting parameter and can be (much) larger than $T_{0}$; this way to interpret experimental data was used, e.g., in Refs. Winnerl et al. (2011); Soavi et al. (2018); Alexander et al. (2018); Tomadin et al. (2018). However, in general, not only the temperature, but also the chemical potentials of electron ($\mu_{e}$) and hole ($\mu_{h}$) gases should be considered to be different from $\mu_{0}$.

The description of a strongly nonequilibrium electron-hole plasma in terms of the quasi-equilibrium electron and hole Fermi gases with their own chemical potentials and temperatures Malic et al. (2011); Sun et al. (2012); Song et al. (2013); Tomadin et al. (2013) is justified if the electron-electron, electron-hole and hole-hole scattering processes (characterized by a typical scattering time $\tau_{ee}$) are more probable than the electron-phonon and electron-impurities ones. There exist theoretical arguments Song et al. (2013); Tomadin et al. (2013) and experimental evidences Lui et al. (2010); Breusing et al. (2011); Brida et al. (2013); Johannsen et al. (2013) that in typical graphene samples this situation is the case indeed.

Although the hot electron model (HEM) has been already used for interpretation of several nonlinear graphene experiments, a comprehensive theory which would analyze different physical situations and would give the opportunity to calculate $\mu_{e}$, $\mu_{h}$ and $T$ as a function of different input parameters of the problem is still absent. In this paper we develop such a theory. In Section II we introduce a HEM and formulate a system of nonlinear differential equations which allows to calculate $\mu_{e}$, $\mu_{h}$, $T$ and other physical quantities characterizing the electron-hole plasma in graphene in the strongly non-equilibrium state. In Section III we analyze solutions of this system of equations in doped and intrinsic graphene, as well as compare results of our theory with some experimental data. In Section IV the results are summarized and conclusions are drawn.

II Theory

II.1 The system in equilibrium

We consider a graphene monolayer lying at the plane $z=0$ on top of a dielectric with the dielectric constant $\kappa$ and the refractive index $n=\sqrt{\kappa}$. The energy spectrum of electrons ($l=2$) and holes ($l=1$) in graphene is

[TABLE]

where $v_{F}\approx 10^{8}$ cm/s is the Fermi velocity and the energy $E_{l{\bm{k}}}$ and the wave vector $\bm{k}$ are counted from one of the Dirac points. In equilibrium (without irradiation) the electron distribution function has the form (the Boltzmann constant $k_{B}=1$ is set to be unity everywhere)

[TABLE]

where $\mu_{0}$ and $T_{0}$ are the equilibrium chemical potential and temperature, the same for electrons and holes.

Below we will analyze two representative cases, with $\mu_{0}=-0.2$ eV (doped graphene) and $\mu_{0}=0$ eV (intrinsic graphene). If $\mu_{0}=-0.2$ eV then at room temperature $T_{0}=300$ K the equilibrium densities of electrons and holes,

[TABLE]

differ by almost five orders of magnitude. In intrinsic graphene at $T_{0}=300$ K the densities are

[TABLE]

The equilibrium chemical potential $\mu_{0}$ can be experimentally varied by the gate voltage.

II.2 Hot electrons distribution functions

Now we assume that graphene is irradiated by an external electromagnetic wave with the frequency $\omega$ and intensity $I$. The photon energy $\hbar\omega$ can be both larger and smaller than $2|\mu_{0}|$, and the intensity of radiation is assumed to be so large that the perturbation theory is inapplicable. The photo-excited electrons absorb the wave energy, due to the intra- and inter-band absorption processes, and relax their energy to the crystal lattice and to the substrate via different scattering processes. We denote the electron-electron (as well as hole-hole and electron-hole) scattering time as $\tau_{ee}$, the momentum and energy intra-band relaxation times, due to the electron scattering by lattice imperfections (phonons, impurities, etc.), as $\tau_{p}$ and $\tau_{\epsilon}$, and the inter-band energy relaxation (actually recombination) time as $\tau_{\rm rec}$ (the time $\tau_{\rm rec}$ will be discussed later in Section II.7). The momentum relaxation time $\tau_{p}$ is typically much smaller than $\tau_{\epsilon}$, $\tau_{p}\ll\tau_{\epsilon}$. Further, we will accept a hypothesis Song et al. (2013); Tomadin et al. (2013); Lui et al. (2010); Breusing et al. (2011); Brida et al. (2013); Johannsen et al. (2013) that the electron-electron scattering time is smaller than $\tau_{p}$,

[TABLE]

according to the literature, $\tau_{ee}$ is about a few tens of fs, while $\tau_{p}$ is at least 0.1 ps or larger. Under these conditions, shortly after the excitation quasi-equilibrium Fermi distributions

[TABLE]

with the electron ($\mu_{e}$) and hole ($\mu_{h}$) chemical potentials and the common temperature $T\neq T_{0}$, are formed in the conduction and valence bands. The distribution functions of electron ($f_{e}=f$) and holes ($f_{h}=1-f$) then read

[TABLE]

[TABLE]

It is also possible to consider the version of the theory in which the temperatures of the electron and hole gases, $T_{e}$ and $T_{h}$, are different. This corresponds to a situation in which electron-electron and hole-hole scattering is more likely than electron-hole scattering. As was shown in Ref. Sun et al. (2012) this is typically not the case, therefore we will restrict ourselves by the model with $T_{e}=T_{h}=T$.

We have introduced three unknown quantities $\mu_{e}$, $\mu_{h}$ and $T$, and now need equations which would determine their dependencies on the equilibrium parameters $\mu_{0}$ and $T_{0}$, as well as on the frequency and intensity of the incident radiation.

II.3 Electron and hole densities

The density of electrons and holes in the strongly non-equilibrium state (6) are determined by the distribution functions (7)–(8) in the usual way,

[TABLE]

[TABLE]

where $g_{s}=g_{v}=2$ are the spin and valley degeneracies, $S$ is the sample area, and the function $F_{n}(z)$ is defined as

[TABLE]

The equilibrium electron and hole densities (3)–(4) are determined by Eqs. (9)–(10) in which $\mu_{e}=\mu_{h}=\mu_{0}$ and $T=T_{0}$.

II.4 Electron and hole energy densities

The energy density of the electron and hole gases per unit area (per cm2) is determined by

[TABLE]

[TABLE]

The total energy of electrons and holes is

[TABLE]

II.5 Conductivity

The linear-response conductivity has three contributions (the derivation can be found, e.g., in Ref. Mikhailov (2016)): intra-band electron, intra-band hole and inter-band,

[TABLE]

where

[TABLE]

[TABLE]

and

[TABLE]

We discuss these contributions separately.

II.5.1 Intra-band conductivity

In order to calculate the intra-band conductivity (16)–(17) we need a model for the scattering rate $\gamma_{\rm intra}(E)\equiv 1/\tau_{p}(E)$, where $\tau_{p}(E)$, see Eq. (5), is the energy-dependent momentum relaxation time due to the scattering of electrons and holes with impurities, phonons and other lattice imperfections (but not with each other). As a first choice we use for $\gamma_{\rm intra}(E)$ the model

[TABLE]

which is discussed in detail in Appendix A. The quantities $\zeta$ and $E_{i}$ in (19) are fitting parameters: $\zeta$ is the minimal static conductivity of graphene in the Dirac point, in units $e^{2}/h$; $E_{i}\propto\sqrt{N_{i}}$ is a Coulomb energy associated with the density of impurities $N_{i}$, see (49). The parameters $\zeta$ and $E_{i}$ can be found by fitting the formula (53) to the experimental data for the gate voltage dependence of the static linear conductivity of graphene, see example in Figure 13. Thus found parameters $\zeta$ and $E_{i}$ and the model expression (19) are then used in formulas (16)–(17) for the high-frequency nonlinear conductivities of graphene. In the rest of the paper we use $\zeta=4$ and $E_{i}=30$ meV, which corresponds to the mobility of about 7260 cm2/Vs. For the relation between the energy $E_{i}$ and the sample mobility, as well as for further discussion of the model (19) see Appendix A.

Alternatively, we also use the energy-independent momentum relaxation rate model with $\gamma_{\rm intra}=\gamma_{p}=1/\tau_{p}$. Then the intra-band dynamic conductivity assumes the form

[TABLE]

The model (19) better reproduces typical experimental data on the static conductivity of graphene, therefore we use it in the main part of the paper. The model (20) was used in some experiments (e.g., Ref. Winnerl et al. (2011)) by interpreting the measured data; comparing our results with Ref. Winnerl et al. (2011) in Section III we also use the model (20) with the energy-independent momentum relaxation rate $\gamma_{\rm intra}$.

II.5.2 Inter-band dynamic conductivity

The inter-band conductivity (18) depends on the inter-band scattering rate $\gamma_{\rm inter}$. We assume that $\gamma_{\rm inter}\to 0$ since a finite $\gamma_{\rm inter}$ does not influence the final result under the condition $\hbar\gamma_{\rm inter}\ll T$ which is typically satisfied in experiments. Then Eq. (18) can be simplified so that the real part assumes the form

[TABLE]

and the imaginary part is expressed in terms of a principal value integral (denoted by ${\cal P}$),

[TABLE]

Notice that the real part of the inter-band conductivity can be negative if $\mu_{e}>\mu_{h}$ and $\hbar\omega<(\mu_{e}-\mu_{h})$. Physically this is due to the population inversion in the non-equilibrium state.

II.6 Absorption coefficient

We assume that graphene lies on the surface of a dielectric substrate with the dielectric constant $\kappa$ and the refractive index $n=\sqrt{\kappa}$, and the external radiation with the intensity $I$ is normally incident on the structure. The incident radiation is transmitted through (with the intensity $TI$), reflected from (the intensity $RI$) and absorbed in the graphene layer (the intensity $AI$). The absorbed part of the radiation energy is determined by the absorption coefficient $A$, sometimes also referred to as absorbance. The coefficient $A$ is determined by the Joule heating $\bm{j}\cdot\bm{E}$ and in the linear-response regime is proportional to the real part of the first-order conductivity $\sigma^{(1)}$. In the nonlinear regime we will assume that the nonlinearity mainly manifests itself in changing the chemical potentials ($\mu_{0}\to\mu_{e},\mu_{h}$) and electron temperature ($T_{0}\to T$) in formulas (15) – (18) and hence, in accordance with these equations, the absorption coefficient can be presented in the form

[TABLE]

with intra-,

[TABLE]

and inter-band,

[TABLE]

contributions. The denominators in Eqs. (24)–(25) contain the total conductivity. If the substrate is made out of silicon dioxide then its dielectric constant is $\kappa_{\rm{SiO}_{2}}=3.9$ and $n=\sqrt{\kappa_{\rm{SiO}_{2}}}=1.975$.

The absorption coefficient formulas (23) – (25) are approximate. In general the current $\bm{j}$ contains the higher contributions $j_{\alpha}^{(3)}=\sigma^{(3)}_{\alpha\beta\gamma\delta}E_{\beta}E_{\gamma}E_{\delta}$, $j_{\alpha}^{(5)}=\sigma^{(5)}_{\alpha\beta\gamma\delta\mu\nu}E_{\beta}E_{\gamma}E_{\delta}E_{\mu}E_{\nu}$, etc., where all higher-order conductivities should be considered as functions of non-equilibrium chemical potentials $\mu_{e}$, $\mu_{h}$ and temperature $T$. However at present the functions $\sigma^{(5)}$, $\sigma^{(7)}$, etc., are unknown and the function $\sigma^{(3)}_{\alpha\beta\gamma\delta}(\omega_{1},\omega_{2},\omega_{3};\mu_{0},T_{0})$ was calculated Cheng et al. (2014b, 2015a); Mikhailov (2016) only for the quasi-equilibrium case with $\mu_{0}$ and $T_{0}$. Therefore in this paper we restrict ourselves by the approach (23) – (25), postponing developing of more general theories for future publications.

II.7 Dynamics equations of the hot electron model

Now we are prepared to formulate the basic equations of our HEM. We will assume that the intra-band energy relaxation time $\tau_{\epsilon}$ is shorter than the inter-band recombination time $\tau_{\rm rec}$, $\tau_{\epsilon}\ll\tau_{\rm rec}$. This condition is typically satisfied in conventional semiconductors. In graphene the radiative recombination time, according to estimates in Refs. Vasko and Ryzhii (2008); Alymov et al. (2018), is around hundreds of nanoseconds at room temperature, while $\tau_{\epsilon}$ lies in the tens-of-ps range. This justifies the use of the condition $\tau_{\epsilon}\ll\tau_{\rm rec}$ below.

II.7.1 Energy relaxation

Assume that the system is excited by a powerful incident radiation with the intensity $I$. Since the recombination is a slow process, the electron and hole Fermi gases are independent from each other in that sense that they are characterized (at the time $t\gtrsim\tau_{ee}$ after the excitation is switched on) by their own chemical potentials $\mu_{e}$ and $\mu_{h}$, and the temperature $T$. At the longer time scale $t\sim\tau_{\epsilon}$ the charge carriers, having been scattered by phonons, impurities and other lattice imperfections, relax their energy to the lattice. We assume that the energy relaxation equations for electrons and holes can then be written, as it is usually done in semiconductor physics, in the form

[TABLE]

[TABLE]

meaning that the energy of hot electron and hole gases grows in time due to the intra-band absorption in each (conduction and valence) band and relaxes to their steady-state quasi-equilibrium energies ${\cal E}_{e}(\mu_{e}^{0},T_{0})$ and ${\cal E}_{h}(\mu_{h}^{0},T_{0})$ with the characteristic time scale $\tau_{\epsilon}$. We emphasize that the temperature of the relaxed quasi-equilibrium state in Eqs. (26)–(27) coincides with the lattice temperature $T_{0}$ since $\tau_{\epsilon}$ describes the relaxation processes between the charge carrier gases and the lattice. However, the chemical potentials $\mu_{e}^{0}$ and $\mu_{h}^{0}$ differ from the equilibrium chemical potential $\mu_{0}$ since the density of electrons and holes are still larger than those in equilibrium since $\tau_{\epsilon}\ll\tau_{\rm rec}$. The relation between $\mu_{e,h}$ and $T$, from one side, and $\mu_{e,h}^{0}$ and $T_{0}$, from the other side, is determined by the conservation of the electron and hole densities,

[TABLE]

The energy relaxation times $\tau_{\epsilon}$ in Eqs. (26)–(27) can, in principle, be different. We will assume, for simplicity, that they are the same. Then we can take a sum of Eqs. (26)–(27) and get the total energy relaxation equation

[TABLE]

II.7.2 Recombination

At a longer time scale $\sim\tau_{\rm rec}$ electrons and holes recombine. Taking into account that they are generated and recombine by pairs,

[TABLE]

we write the generation-recombination rate equation in the form

[TABLE]

The first (generation) term $G$ in the right hand side of (31) represents the number of electron-hole pairs generated per second on a unit area. It equals the radiation intensity $I$ (the radiation energy incident on a unit area per second), times the inter-band absorption coefficient (which gives the energy absorbed on a unit area per second due to the electron-hole generation processes), and divided by the photon energy (which results in the number of electron-hole pairs generated on a unit area per second). The second term $R$ in the right hand side of (31) is the recombination rate. The recombination is a nonlinear bi-particle process with the recombination rate being proportional to the product of electron and hole densities $n_{e}n_{h}$. The recombination term describes the relaxation to the equilibrium electron and hole densities $n_{e}^{0}n_{h}^{0}$. The recombination coefficient $\alpha_{\rm rec}$ is measured in units cm2/s and is independent of the particle densities. It is the second (in addition to $\tau_{\epsilon}$) parameter of the theory.

Apart from the recombination coefficient $\alpha_{\rm rec}$ one can also introduce a quantity $\tau_{\rm rec}$ which is measured in units of time and at low excitation levels has the meaning of the recombination time (in general the recombination process is not purely exponential and the meaning of $\tau_{\rm rec}$ is more complicated, see below). Assume that the radiation intensity $I(t)$ is switched off at the time moment $t=0$ and consider the time evolution of the electron and hole densities $\delta n(t)$ at $t>0$. Substituting (30) into equation (31) with $I=0$ and taking into account the initial condition $\delta n(0)=\delta n_{0}$ we get

[TABLE]

where $\Delta(t)$ is the relative change of the charge carrier density as compared to their total equilibrium density, and $\Delta_{0}=\Delta(0)$. The quantity

[TABLE]

has the dimension of time, depends on the total equilibrium density of electrons and holes and determines the time evolution of the electron-hole recombination (32). We emphasize that $\tau_{\rm rec}$ has the meaning of time over which the initial carrier density reduces by a factor of $\sim 2.72...$ only at very low excitation levels $\Delta_{0}\ll 1$,

[TABLE]

see Figure 1. At high excitation levels, $\Delta_{0}\gg 1$ or $\delta n_{0}\gg(n_{e}^{0}+n_{h}^{0})$, the density $\delta n(t)$ first very quickly decreases, with the time constant $\sim\tau_{\rm rec}/\Delta_{0}$, down to the values $\sim(n_{e}^{0}+n_{h}^{0})$, and then decays further exponentially, see inset to Figure 1. Quantitatively, $\Delta(t)$ decreases from its initial value $\Delta_{0}$ by a factor of two during the time $\tau_{\rm rec}/(\Delta_{0}+1)\ll\tau_{\rm rec}$ and by a factor of $\Delta_{0}$ during the time $\sim\tau_{\rm rec}$.

The generation-recombination rate equation (31) can be also rewritten in the form explicitly containing $\tau_{\rm rec}$,

[TABLE]

II.7.3 Preliminary summary and discussion

The recombination of charge carriers in graphene characterized by a more complicated than $e^{-t/\tau}$-decay was commonly observed in time-resolved pump-probe experiments, see, e.g., Refs. Dawlaty et al. (2008); George et al. (2008); Winnerl et al. (2011) and other. It was often interpreted by introducing two different time scales $\tau_{1}$ and $\tau_{2}$ where different $\tau$-s were attributed to physically different relaxation mechanisms. As seen from Section II.7.2 the seemingly double-$\tau$ time decay is actually described by a single formula (32) with only one decay-time parameter $\tau_{\rm rec}$. The reason of the more complicated behavior of $\delta n(t)$ is the intrinsically nonlinear nature of the electron-hole recombination process seen in Eqs. (31), (35). The corresponding two time constants are $\tau_{1}=\tau_{\rm rec}/(\Delta_{0}+1)$ and $\tau_{2}=\tau_{\rm rec}$. In the strong excitation limit $\Delta_{0}\gg 1$ the first time is much shorter than the second one, $\tau_{1}\ll\tau_{2}$; while in the weak excitation limit ($\Delta_{0}\ll 1$) the two times merge into one, $\tau_{1}\simeq\tau_{2}=\tau_{\rm rec}$, Figure 1.

The nonlinearity of the recombination process is known in the semiconductor physics, e.g., Ref. Bonch-Bruevich and Kalashnikov (1977). We have briefly reproduced here the nonlinear recombination equations (31), (35) and the corresponding derivation of Eq. (32) since in some recent papers (e.g., Ref. Soavi et al. (2018)) a strongly non-equilibrium ($\Delta_{0}\gg 1$) recombination dynamics has been improperly described by a linear recombination term $R\propto-\delta n/\tau_{\rm rec}$.

Equations (29) and (35) [or (29) and (31)] describe the dynamics of the electron temperature $T$, chemical potentials $\mu_{e}$, $\mu_{h}$, and all other physical quantities within our HEM. The energy densities ${\cal E}_{e}$, ${\cal E}_{h}$ in Eq. (29), as well as the intra- and inter-band absorption coefficients $A^{\rm intra}$, $A^{\rm inter}$ in Eqs. (29) and (35), depend on six unknown quantities $T$, $\mu_{e}$, $\mu_{h}$, $\mu_{e}^{0}$, $\mu_{h}^{0}$, and $\delta n$. The four missing equations can be found by inverting the relations (9) and (10), namely,

[TABLE]

[TABLE]

[TABLE]

[TABLE]

where $G(y)$ is the inverse function of $F_{1}(x)$,

[TABLE]

with $F_{1}$ defined in (11). Equations (36)–(39), together with (29) and (35), give six equations for six unknown variables. The quantities $\omega$, $I$, $T_{0}$ and $\mu_{0}$ (or $E_{F}$), as well as parameters $\tau_{\epsilon}$ and $\tau_{\rm rec}$ are assumed to be known input parameters.

II.8 Steady-state equations

In the rest of the paper we will analyze the steady-state solutions of Eqs. (29), (35). Setting $\partial/\partial t=0$ we get

[TABLE]

The second equation here does not contain the intensity $I$. It can be presented in the dimensionless form

[TABLE]

where the left- and right-hand sides read

[TABLE]

[TABLE]

Solving the nonlinear equation (42), together with (36)–(39), we can find the relation between the relative change of the density $\Delta=\delta n/(n_{e}^{0}+n_{h}^{0})$ and the relative change of the temperature $\delta T/T_{0}=(T-T_{0})/T_{0}$. Notice that equation (42) depends only on the ratio of the characteristic times $\tau_{\epsilon}/\tau_{\rm rec}$ and not on each of them separately. After the relation between $\Delta$ and $\delta T/T_{0}$ is found we substitute it into any of the equations (41) and relate these quantities to the intensity of radiation. This gives two equivalent formulas

[TABLE]

or

[TABLE]

where we have introduced a power density unit

[TABLE]

which does not depend on the charge carrier density. At room temperature the power density unit $I_{0}$ is about 400 W/cm2 if $\tau_{\epsilon}\simeq 1$ ps. This is a rather small value, i.e., the strongly nonlinear regime corresponds to $I\gg I_{0}$. Notice that Eqs. (42)–(46) are presented in the explicitly dimensionless form, which means, in particular, that scaling all energies by the same numerical factor will not change the final results. This makes them universal in a sense.

Now we can start analyzing results which our model gives. We will assume that $T_{0}=300$ K and $\mu_{0}=-0.2$ eV which corresponds to the equilibrium charge carriers density of about $3\times 10^{12}$ cm*-2*, see (3). All energy quantities will be given in eV.

III Results

III.1 Density-temperature diagrams

Figure 2 shows the density-temperature diagrams obtained by solving equation (42) at different frequencies, $T_{0}=300$ K, $\tau_{\epsilon}/\tau_{\rm rec}=0.1$, and two different values of the equilibrium chemical potenial $\mu_{0}=-0.2$ and [math] eV. Each point on each curve corresponds to a certain value of the input wave power density. The low-intensity regime $I/I_{0}\ll 1$ corresponds to the origin of the plots where $\Delta\ll 1$ and $\delta T/T_{0}\ll 1$. The higher the power, the larger is the relative changes of both the density and the temperature, so that both $\Delta$ and $\delta T/T_{0}$ grow with the increasing intensity. The rate of their growth depends however on the frequency.

First, we consider the case $\mu_{0}=-0.2$ eV, Figure 2(a). If $\hbar\omega\gtrsim 2|\mu_{0}|=0.4$ eV (black to blue curves), i.e., when the inter-band transitions dominate, the relative change of the density $\Delta$, at low intensities, is much stronger than the relative change of the temperature $\delta T/T_{0}$. When $I/I_{0}$ grows from zero up to the values $\sim 10^{5}$ (see arrows in Figure 2(a)), $\Delta$ increases up to $\Delta\simeq 1-3$, while the temperature remains practically unchanged, $\delta T/T_{0}\ll 1$. The values $\Delta\simeq 1-3$ correspond to practically equal densities of photo-excited electrons and holes; remind that at $\mu_{0}=-0.2$ eV and $T_{0}=300$ K the equilibrium electron and hole densities differed by almost five orders of magnitude, Eq. (3). At even larger intensities, $I/I_{0}\gtrsim 10^{5}$, the hot carrier temperature starts to grow too: in this situation the density of electrons and holes are large and close to each other and the both gases are heated by the intra-band absorption in the corresponding energy bands.

In the regime $\hbar\omega\ll 2|\mu_{0}|$ the tendency is opposite. The temperature $\delta T/T_{0}$ grows much faster than the relative change of the density $\Delta$ which is physically clear since at $\hbar\omega\ll 2|\mu_{0}|$ the intra-band transitions dominate. At very large power densities, however, both $\Delta$ and $\delta T/T_{0}$ becomes quite comparable to each other.

If the case of intrinsic graphene, Figure 2(b), the relative density changes are always much larger than the relative change of temperature. The absolute values of $\Delta$ are also much larger than in Figure 2(a). These two features are a simple consequence of the fact that at $\mu_{0}=0$ eV the inter-band transitions are always dominant and that in equilibrium the density of charge carriers (4) is much smaller than in the doped one. For example, $\Delta\simeq 100$ in Figure 2(b) corresponds to approximately the same values of the non-equlibrium electron and holes densities ($n_{e}=n_{h}\approx 1.63\times 10^{13}$ cm*-2*) as $\Delta\simeq 4$ in Figure 2(a) ($n_{e}\approx 1.25\times 10^{13}$ cm*-2*, $n_{h}\approx 1.56\times 10^{13}$ cm*-2*).

III.2 Power dependencies of different physical quantities

Now we consider how different physical quantities vary with the radiation intensity. We show results for doped ($\mu_{0}=-0.2$ eV) and intrinsic graphene ($\mu_{0}=0$ eV).

III.2.1 Doped graphene

Figure 3 exhibits the power dependencies of the chemical potentials $\mu_{e}$ and $\mu_{h}$, temperature $T$, photo-excited charge carrier density $\delta n$ and the absorption coefficient $A$, in the doped graphene sample with $\mu_{0}=-0.2$ eV, under the condition $\hbar\omega>2|\mu_{0}|$ when the inter-band transitions dominate. The ratio of relaxation times is assumed to be $\tau_{\epsilon}/\tau_{\rm rec}=0.1$ and we consider a relatively high-mobility sample with $E_{i}=30$ meV (this corresponds to $N_{i}\approx 6.4\times 10^{11}$ cm*-2* and $\mu\approx 7260$ cm2/Vs, see Fig. 12). At low intensities $I/I_{0}\lesssim 10^{4}$, left panels, the density $\delta n$ substantially changes, as expected, while the temperature remains practically unchanged. The chemical potential of electrons quickly grows, from the initial value $\mu_{0}=-0.2$ eV, and becomes positive at $I/I_{0}\simeq 200$; the chemical potential of holes becomes more negative and varies slowly. The absorption coefficient does not change with the intensity up to $I/I_{0}\simeq 10^{3}$, but starts to decrease when $I/I_{0}$ approaches the values of order of $10^{4}$. It is mainly due to the inter-band contribution, and is about 1%. This number differs from traditional 2.3% since we consider graphene lying on a SiO2 substrate, Eqs. (24)–(25).

When the intensity grows further, electron and holes gases get heated and the charge carrier temperature increases too. The relative changes of density and temperature becomes equal at $I/I_{0}\approx 5.4\times 10^{6}$, Figure 3(a), right panel. At the point, where the black and red curves intersect, $\Delta=\delta T/T_{0}\approx 2.14$, which corresponds to $n_{e}\approx 6.68\times 10^{12}$ cm*-2*, $n_{h}\approx 9.80\times 10^{12}$ cm*-2*, and $T\approx 3.14T_{0}\approx 942$ K. The chemical potential of electrons continues to grow reaching the values $\simeq+0.26$ eV at $I/I_{0}\simeq 10^{6}$. The chemical potential of holes gets more negative and becomes equal to $\simeq-0.33$ eV at $I/I_{0}\simeq 10^{6}$. Since the occupation of electrons and holes states around $\hbar\omega/2\simeq 0.3$ becomes much less asymmetric as compared to equilibrium the inter-band absorption starts to fall down at $I/I_{0}\simeq 10^{4}$ and becomes about $0.04$% at $I/I_{0}\simeq 10^{6}$. This strong reduction of the absorption coefficient is usually referred to as the saturable absorption effect. The intra-band absorption remains small as compared to the inter-band one.

Now we consider the case $\hbar\omega<2|\mu_{0}|$ where the intra-band absorption plays the crucial role at low intensities. Figure 4 shows the power dependencies of different physical quantities under the same conditions as in Figure 3 but at $\hbar\omega=0.1$ eV. Now the temperature substantially grows at low intensities while the charge carrier density remains almost unchanged up to $I/I_{0}\simeq 10^{5}$, Figure 4(a), right panel. The chemical potential of electrons sharply grows at the radiation power $I/I_{0}\lesssim 700$ but then saturates at the much lower level $\mu_{e}\simeq-0.1$ eV than in Figure 3. The chemical potential of holes remains almost constant slightly decreasing in the absolute value. The absorption coefficient is about 0.12% at low intensities and remains practically constant up to $I/I_{0}\simeq 10^{5}$, Figure 4(c). It is mainly due to the intra-band contribution which is much smaller than in the previous example since the frequency lies in the gap between the intra- and inter-band absorption areas, $\gamma(\mu_{0})\ll\omega\ll 2|\mu_{0}|/\hbar$, where $\hbar\gamma(\mu_{0})\approx 4.5$ meV under our conditions. At $I/I_{0}\gtrsim 10^{5}$ the absorption coefficient starts to grow (the induced absorption Winnerl et al. (2011)), mainly due to the inter-band contribution which becomes essential since the occupation of energy levels at $E\simeq-\hbar\omega/2$ is no longer negligible due to the heating of the hole gas. At $I/I_{0}\gtrsim 10^{5}$ also the chemical potentials of both electrons and holes start to substantially grow making the distribution of charge carriers over the bands more uniform; at $I/I_{0}$ between $\simeq 10^{6}$ and $10^{8}$ the chemical potential of electrons even becomes positive, Figure 4(b), right panel. At even larger intensities $I/I_{0}\gtrsim 10^{8}$ the chemical potential of holes moves to the conduction band, $\mu_{h}>0$, while that of electrons becomes negative again, Figure 4(b), right panel.

In Figure 5 we further illustrate our results by showing the electron distribution function in the valence and conduction bands at different radiation intensities. Here one clearly sees a qualitative difference between the charge carrier distributions at $\hbar\omega>2|\mu_{0}|$, Figure 5(a), and at $\hbar\omega<2|\mu_{0}|$, Figure 5(b). In the inter-band absorption case $\hbar\omega>2|\mu_{0}|$ the slope of the $f(E)$ curves, and hence the charge carrier temperature $T$, remains practically unchanged up to intensities $I/I_{0}\simeq 10^{5}$. In contrast, the chemical potentials vary quite strongly: already at $I/I_{0}\simeq 10^{3}$ (red curve) the chemical potential of electrons is positive and the occupation of the conduction band is quite large. At higher intensities, $I/I_{0}\gtrsim 10^{6}$, when the densities of photo-excited electrons and holes become comparable, the temperature starts to grow too due to the intra-band absorption in each band, and the slope of the $f(E)$ curves decreases.

In the intra-band absorption case $\hbar\omega<2|\mu_{0}|$, Figure 5(b), the slope of the curves noticeably decreases already at $I/I_{0}\gtrsim 10^{4}$, and the occupation of the conduction band is much weaker than in the previous case (the chemical potential of electrons remains negative). Only at the very high intensity $I/I_{0}\simeq 10^{7}$ (magenta curve) $\mu_{e}$ becomes slightly positive.

III.2.2 Intrinsic graphene

Figure 6 shows the power dependencies of different physical quantities in intrinsic graphene with $\mu_{0}=0$ and $\hbar\omega=0.4$ eV. Other parameters ($T_{0}$, $E_{i}$ and $\tau_{\epsilon}/\tau_{\rm rec}$) are the same as in Figures 3, 4 and 5. Since at $\mu_{0}=0$ the condition $\hbar\omega>2|\mu_{0}|$ is always satisfied the curves shown in Figure 6 are qualitatively similar to those from Figure 3. The inter-band transitions dominate, therefore the relative change of the density $\Delta$ is always much larger than the relative change of temperature, Figure 6(a). Only at $I/I_{0}\gtrsim 10^{6}$ the temperature $T$ starts to noticeably grow, Figure 6(a), right panel. The chemical potentials of electrons and holes are symmetric, $\mu_{e}=-\mu_{h}$, and achieve the values of order $|\mu_{e}|=|\mu_{h}|\simeq 0.1$ eV at $I/I_{0}\gtrsim 10^{4}$ and $\simeq 0.2$ eV at $I/I_{0}\gtrsim 10^{5}-10^{7}$, Figure 6(b). The absorption curves show the saturable absorption effect: the absorption coefficient $A$ falls down from the value $\simeq 1$ % at $I/I_{0}\lesssim 10^{4}$ down to the values $\lesssim 0.12$ % at $I/I_{0}\gtrsim 10^{6}$. The contribution of the intra-band absorption to $A$ is negligibly small.

Figure 7 shows the non-equilibrium distribution function of electrons and holes at $\hbar\omega=0.4$ eV and different power levels. In accordance with Figure 6 the temperature remains low ($T\simeq T_{0}$) at the intensities up to $I/I_{0}\simeq 10^{4}$, while the chemical potentials $|\mu_{e,h}|$ grow. At higher intensities the temperature increases too (the cyan and, especially, magenta curves).

III.3 Frequency dependencies of different physical quantities

Another way to clarify the physics of the discussed phenomena is to analyze how the frequency dependencies of the absorption coefficient $A(\omega)$ and other physical quantites are modified under the influence of the strong radiation power. This is especially important question since the spectra $A(\omega)$ can be directly experimentally measured Li et al. (2008).

Figure 8 shows the power-dependent absorption spectra in (a) doped ($|\mu_{0}|=0.2$ eV) and (b) intrinsic ($\mu_{0}=0$ eV) graphene, at $\tau_{\epsilon}/\tau_{\rm rec}=0.1$, $E_{i}=30$ meV and $T_{0}=300$ K. The black curves in both panels exhibit the known linear response absorption curves (e.g., Refs. Mikhailov (2016); Li et al. (2008)) corresponding to the equilibrium chemical potential, $\mu_{e}=\mu_{h}=\mu_{0}$, and equilibrium temperature, $T=T_{0}$. One sees that the growing power substantially reduces the absorption in graphene. This saturable absorption effect is the case already at $\hbar\omega\gtrsim|\mu_{0}|$ (not $2|\mu_{0}|$!) in Figure 8(a) and at $\hbar\omega\gtrsim 0.04$ eV in Figure 8(b). The influence of the radiation power is noticeable at $I/I_{0}\simeq 10^{3}$ in the doped graphene, Figure 8(a), and at even lower powers ($I/I_{0}\simeq 10^{1}-10^{2}$) in the intrinsic graphene, Figure 8(b). Quantitatively, the suppression of $A(\omega)$ is very strong; for example, in doped graphene at $\hbar\omega=0.6$ eV and $I/I_{0}=10^{6}$ the absorption is only 0.042%, i.e. it is reduced by a factor of $\sim 24$. In the intrinsic graphene at $\hbar\omega=0.3$ eV and $I/I_{0}=10^{6}$ the absorption is reduced down to 0.085%, i.e. by a factor of $\sim 12$. The saturable absorption effect at high frequencies ($\hbar\omega\gtrsim|\mu_{0}|$) was experimentally observed in many experiments, see, e.g., Refs. Bao et al. (2009); Zhang et al. (2009); Winnerl et al. (2011); Zheng et al. (2012); Bianchi et al. (2017).

At lower frequencies ($\hbar\omega\lesssim|\mu_{0}|$) our model predicts an essentially different behavior. The absorption spectrum weakly depends on the radiation power at $\hbar\omega\lesssim|\mu_{0}|$ and the radiation may lead to a slight increase of the absorption. In doped graphene, Figure 8(a), this effect is rather small; for example, at $\hbar\omega=0.1$ eV the absorption coefficient is about 11.79% at $I/I_{0}\ll 1$ and increases by $\simeq 0.3$% at $I/I_{0}=10^{5}$ and by $\simeq 0.45$% at $I/I_{0}=10^{6}$. In intrinsic graphene, Figure 8(b), the growth of absorption is stronger and can achieve $2-3$% at $I/I_{0}$ up to $\simeq 10^{4}$: for example, at $\hbar\omega=0.01$ eV the absorption is about 2.3% at $I/I_{0}\ll 1$, 4.8% at $I/I_{0}=10^{3}$ and 5.8% at $I/I_{0}=10^{4}$. The specific numbers of the absorption change depend of course on the chosen parameters of the structure.

Physically the growth of absorption at $\hbar\omega\lesssim 2|\mu_{0}|$ and its reduction at $\hbar\omega\gtrsim 2|\mu_{0}|$ are explained by the radiation induced redistribution of electrons over quantum states in the conduction and valence bands, Figure 9. Notice that this qualitative picture, which appeared in many publications (see, e.g., Ref. Winnerl et al. (2011)), implies that the hot electron temperature $T$, as well as the chemical potentials of electrons and holes $\mu_{e}$ and $\mu_{h}$, essentially depend on the photon energy. Within our HEM we can quantitatively evaluate these dependencies. In Figure 10 we plot the hot electron temperature $T$ as a function of the photon energy for a few sets of experimental parameters. The four curves in the upper right corner show the dependencies $T(\omega)$ for parameters corresponding to Figure 8(a) ($|\mu_{0}|=0.2$ eV, $T_{0}=300$ K, $E_{i}=30$ meV and $\tau_{\epsilon}/\tau_{\rm rec}=0.1$) and four different power levels. The arrow labeled as $2|\mu_{0}|$ indicates the position of the double chemical potential. One sees that at $\hbar\omega\gtrsim 2|\mu_{0}|$ all curves tend to the equilibrium temperature value $T\to T_{0}=300$ K. Even if at low frequencies the hot electron temperature exceeds $T_{0}$ by more than one order of magnitude, at $\hbar\omega\simeq 2|\mu_{0}|$ it is already almost equal to $T_{0}$; for example, for $I/I_{0}=10^{5}$ (green curve) the temperature $T$ drops from 4486 K at $\hbar\omega=10$ meV down to 380 K at $\hbar\omega=2|\mu_{0}|=0.4$ eV.

In a recent experiment on multilayer epitaxial graphene Winnerl et al. (2011) the influence of a strong (pump) radiation on the transmission coefficient of the weak (probe) wave was studied, and a few-percent increase (decrease) of absorption was observed at $\hbar\omega\lesssim 2|\mu_{0}|$ ($\hbar\omega\gtrsim 2|\mu_{0}|$). Graphene was weakly doped ($\mu_{0}$ was evaluated to be $\simeq 13$ meV) and the experiment was performed at 10 K. The pump radiation with the fluence up to $\simeq 1$ $\mu$J/cm2 reduced the absorption coefficient by a few percent at $\hbar\omega=30$ meV and increased it by a few percent at $\hbar\omega=20$ meV. The authors interpreted their results applying a simplified HEM which assumed that the chemical potentials remain unchanged under the action of radiation, $\mu_{e}=\mu_{h}=\mu_{0}$, and the hot electron temperature does not depend on the radiation frequency, $T(\omega)=$const. A similar effect was observed and the same interpretation was applied to its explanation in Ref. Alexander et al. (2018), where the nonlinear absorption in graphene was measured under different equilibrium conditions with $\mu_{0}=-0.4$ eV and $T_{0}=300$ K.

As we have seen above, in reality the quantities $T$, $\mu_{e}$, and $\mu_{h}$ cannot be considered as frequency independent and a more general theory should be applied. In Figure 11(a) we plot the absorption spectra calculated for $|\mu_{0}|=16$ meV, $T_{0}=10$ K, $\tau_{p}=300$ fs, and $\tau_{\epsilon}/\tau_{\rm rec}=0.01$; the results were found to be weakly dependent on $\tau_{\epsilon}/\tau_{\rm rec}$ as long as this parameter is small as compared to unity. To calculate the graphene conductivity and the absorption coefficient we used the model of the energy independent momentum scattering time (20) with $\tau_{p}$ taken from Ref. Winnerl et al. (2011). One sees that the curves corresponding to different power levels intersect at one point lying approximately at $\hbar\omega\approx 28$ meV. This value is smaller than $2|\mu_{0}|=32$ meV (in contrast to the results of the simplified HEM, see Ref. Winnerl et al. (2011)); this difference results from strong frequency dependence of $T(\omega)$, shown in the lower left corner of Figure 10. Right and left from 28 meV the absorption decreases and increases respectively, and we have chosen the chemical potential $|\mu_{0}|=16$ meV to get approximately equal (in absolute values) changes of the absorption coefficient at $\hbar\omega=20$ and $30$ meV. This quantity is slightly larger than $|\mu_{0}|=13$ meV extracted in Ref. Winnerl et al. (2011) from the comparison of experimental data with the simplified HEM; the reason is again due to the frequency dependent (independent) hot electron temperature in the full (simplified) HEM. In general, one sees that our HEM gives reasonable results which can be used for analysis of different nonlinear optics experiments.

As seen from Figure 11(a) the reduction of absorption right from the intersection point (28 meV) is in general stronger than its growth left from this point. Qualitatively this is also in agreement with the experimental results of Ref. Winnerl et al. (2011) (see Figures 3(a),(b) there). Physically the increase of absorption at low frequencies is due to the larger intra-band contribution to $A$. This contribution can be increased in samples with a lower mobility: in graphene layers with a smaller values of $\tau_{p}$, i.e., with a lower mobility, the hot electron temperature and the absorption increase at low frequencies should be larger. This is confirmed indeed in Figure 11(b) which show the absorption spectra for the same parameters as in Figure 11(a) but for three times smaller $\tau_{p}=100$ fs. One sees that the absorption change increase by several times at $\hbar\omega=20$ meV while at $\hbar\omega\gtrsim 30$ meV changes are less dramatic. Figure 10 also confirms that at frequencies around $2|\mu_{0}|$ the temperature $T$ is larger in samples with $\tau_{p}=100$ fs than in those with $\tau_{p}=300$ fs.

IV Summary

Experiments on the nonlinear graphene optics are very often performed at so strong excitation powers that the perturbation theory fails to adequately describe their results. The hot electron model presented in this paper allows to calculate the most important parameters of highly non-equilibrium charge carriers in graphene – the chemical potentials of electrons and holes, as well as their effective temperature – thus enabling to correctly describe its response to the powerful electromagnetic radiation. The model is physically transparent and employs essentially one fitting parameter – the ratio $\tau_{\epsilon}/\tau_{\rm rec}$ of the intra-band energy relaxation time $\tau_{\epsilon}$ to the inter-band recombination time $\tau_{\rm rec}$. The derived system of strongly nonlinear differential equations (29), (35)–(39) allows to calculate all physical quantities characterizing the nonlinear graphene response as functions of the incident wave frequency and power, equilibrium temperature, doping level, sample mobility, dielectric environment and so on. The developed theory reasonably describes available experimental data. Together with the already published perturbative theories of the nonlinear graphene response Cheng et al. (2014b, 2015a); Mikhailov (2016) the work done here paves the way to a more accurate interpretation of nonlinear optics experiments and to the development of new optoelectronic devices for visible, infrared and terahertz spectral ranges.

Acknowledgements.

This work has received funding from the European Union s Horizon 2020 research and innovation programme Graphene Core 2 under Grant Agreement No. 785219.

Appendix A A model for the intra-band scattering rate

In order to calculate the intra-band conductivity we need a model for the energy dependent intra-band momentum relaxation rate $\gamma_{\rm intra}(E)$. It is known Das Sarma et al. (2011) that the most important scattering mechanism of electrons in graphene is the charged impurity scattering, and that at high energies $\gamma_{\rm intra}(E)$ is proportional to $N_{i}/|E|$, where $N_{i}$ is the impurity density. This can be written in the formTrushin and Schliemann (2007); Hwang and Das Sarma (2009)

[TABLE]

where

[TABLE]

is the Coulomb energy associated with the impurity density $N_{i}$, $\kappa$ is the effective dielectric constant of the medium surrounding the graphene layer, and $\alpha$ is a number of order unity. In Ref. Hwang and Das Sarma (2009) the formula $\alpha=2\sqrt{I_{0}(r_{s})}$ was derived, where

[TABLE]

is the effective fine structure constant of graphene and

[TABLE]

If graphene lies on the surface of silicon dioxide, $\kappa_{{\rm SiO}_{2}}=3.9$, the effective dielectric constant is $\kappa=(\kappa_{{\rm SiO}_{2}}+1)/2=2.45$ and $r_{s}\approx 0.896$. Then $\alpha\approx 0.36$.

The formula (48) is not valid at small energies (below $E_{i}$). We assume that it can be generalized as follows

[TABLE]

where $\zeta$ is a number of order unity. Substituting the model expression (52) into the equilibrium ($\mu_{e}=\mu_{h}=\mu_{0}$) static ($\hbar\omega=0$) zero-temperature ($T=0$, $\mu_{0}=E_{F}$) intra-band conductivity (16) or (17) we obtain

[TABLE]

where $n_{s}=k_{F}^{2}/\pi=E_{F}^{2}/\pi(\hbar v_{F})^{2}$ is the charge carrier density at $T=0$ and

[TABLE]

is another numerical factor. As seen from (53) the factor $\zeta$ determines the minimal conductivity of graphene at the Dirac point; in typical experiments $\zeta\simeq 4$. At large densities, $\beta n_{s}\gg N_{i}$, the conductivity in (53) is proportional to $n_{s}$, $\sigma^{(1)}_{\rm intra}(0,E_{F},0)=e\mu n_{s}$, which gives the relation between the low-temperature mobility $\mu$ and the density of impurities,

[TABLE]

Figure 12 illustrates the relations (49) and (55) between the impurity density $N_{i}$, the energy $E_{i}$ and the mobility $\mu$.

Varying two adjustable parameters $\zeta$ and $E_{i}$ one can now fit the expression (53) to experimental data on the density (or gate-voltage) dependence of the intra-band static conductivity and then use thus found parameters $\zeta$ and $E_{i}$ for calculations of the high-frequency linear and nonlinear response. In Figure 13 we illustrate this procedure by fitting (53) to some of the experimental data from Ref. Chen et al. (2008). In that paper the authors measured the graphene conductivity as a function of gate voltage in samples intentionally doped by potassium atoms. The black and magenta symbols in Figure 13 show the data from Figure 2 of Ref. Chen et al. (2008), for pristine graphene (denoted as 0 s) and for the same sample after the doping during 12 s. The data from Ref. Chen et al. (2008) are replotted as a function of charge carrier density

[TABLE]

where the gate capacitance per unit area $c_{g}=1.15\times 10^{-4}$ F/m2 is taken from Ref. Chen et al. (2008). One sees that the curves (53) excellently reproduce experimental data at reasonable values of the fitting parameters $\zeta$ and $N_{i}$. The found values of $N_{i}$ are about $N_{i}\approx 0.54\times 10^{12}$ cm*-2* for pristine graphene (black symbols) and $N_{i}\approx 3.8\times 10^{12}$ cm*-2* for the doping time of 12 s (magenta symbols). The maximum impurity density (for 18 s doping time) was estimated in Ref. Chen et al. (2008) as $(1.4-1.8)\times 10^{-3}$ potassium per carbon, or $N_{i}^{\rm max}\simeq(5.3-6.9)\times 10^{12}$ cm*-2*, which agrees very well with the numbers obtained from our fit.

In the main text we use the energy $E_{i}$ as a fitting parameter, instead of $N_{i}$. If $N_{i}$ varies in the range $\sim(0.5-5)\times 10^{12}$ cm*-2*, the energy $E_{i}$ lies in the interval from $\sim 26.5$ to $\sim 83.9$ meV and the mobility in the interval from $9300$ to $930$ cm2/Vs.