Nonlinear electric conductivity and THz-induced charge transport in graphene

Recent developments of laser technologies have enabled the study of light-induced nonequilibrium electron dynamics in matter [1–7]. Among various intriguing materials, graphene has been attracting a large amount of interest owing to its unique electronic structure, the so-called Dirac cone. Various light-induced phenomena in graphene have been studied, such as the light-induced topological phase [8], high-order harmonic generation [9–13], and light-induced anomalous Hall effects [8, 14, 15]. Furthermore, nonlinear light-induced current injection in graphene has been investigated toward the achievement of ultrafast optoelectronic devices [16–18].

To control electric current and charge transport by light, a deep understanding of the field-induced nonequilibrium electron dynamics in matter is indispensable. The nonlinear conductivity of graphene in the THz regime has been investigated, and its field-induced transparency has been experimentally reported [19–24], reflecting the reduction in the conductivity caused by a strong field. Based on the semi-classical kinematic theory, the conductivity reduction has been explained through the change in the carrier scattering rate [25, 26]. By contrast, an enhancement of the conductivity has been suggested in a strong field regime via the Landau–Zener tunneling, which is based on the quantum nature of electrons in solids [27]. Therefore, the microscopic physics behind the nonlinear conductivity of graphene still needs to be clarified via a full quantum mechanical description in order to provide a comprehensive understanding of its nonlinear opto-electronic properties.

In this work, the nonlinear electronic current in graphene is theoretically studied under strong static fields using the quantum master equation approach. As a result of the fully-quantum mechanical simulations, it is shown that the conductivity reduction in graphene can be understood in terms of the depletion of the effective carriers without considering the phenomenological change of the scattering rates. Furthermore, the depletion of carriers is found to be overcome by the Landau–Zener tunneling in the strong field regime, resulting in an enhancement of the effective conductivity as well as the induced current. On account of the highly nonlinear behavior of the conductivity, a method to induce charge transport via few-cycle THz laser pulses is proposed, providing a possible foundation for realizing ultrafast optoelectronics with graphene.

Our calculation method has been described in detail elsewhere [13, 15, 28], so we briefly describe the theoretical modeling here. To describe the light-indued electron dynamics in graphene, the following quantum master equation was employed [15, 28],

$$\begin{equation}\frac{\mathrm{d}}{\mathrm{d}t}\enspace {\rho }_{\boldsymbol{k}}\left(t\right)=\frac{1}{i\hslash }\left[{\rho }_{\boldsymbol{k}}\left(t\right),{H}_{\boldsymbol{k}}\left(t\right)\right]+\hat{D}\left[{\rho }_{\boldsymbol{k}}\left(t\right)\right],\end{equation} \tag{ 1 }$$

where ρ _k (t) is the one-body reduced density matrix at each k -point, $\hat{D}\left[{\rho }_{\boldsymbol{k}}\left(t\right)\right]$ is the relaxation operator, and H _k (t) is the Hamiltonian. The simple two-band approximation was employed for the electronic structure of graphene as

$$\begin{equation}{H}_{\boldsymbol{k}}\left(t\right)={v}_{\mathrm{F}}{\sigma }_{x}\left[{k}_{x}+{A}_{x}\left(t\right)\right]+{v}_{\mathrm{F}}{\sigma }_{y}\left[{k}_{y}+{A}_{y}\left(t\right)\right],\end{equation} \tag{ 2 }$$

where σ_j are the Pauli matrices, k_j are the jth components of the Bloch wave vector k , and A_j (t) is the j component of the vector potential A (t), which is related to the external electric field, according to $\boldsymbol{E}\left(t\right)=-\dot {\boldsymbol{A}}\left(t\right)$ in the dipole approximation. The Fermi velocity v_F was set to 1.12 × 10⁶ m s⁻¹ in accordance with a previous ab initio calculation [29]. The relaxation operator $\hat{D}\left[{\rho }_{\boldsymbol{k}}\left(t\right)\right]$ was constructed using the relaxation time approximation [30]. $\hat{D}\left[{\rho }_{\boldsymbol{k}}\left(t\right)\right]$ depends on the longitudinal relaxation time T₁, transverse relaxation time T₂, electron temperature T_e, and chemical potential μ. T₁, T₂, and T_e were set to 100 fs, 20 fs, and 300 K, respectively, whereas μ was treated as a tunable parameter. Our treatment of the relaxation operator has been described in detail elsewhere [13, 15, 28]. As will be shown later (see figures 2(c)–(f)), we consider the electron dynamics only in a narrow region around the Dirac point. Therefore, the linear band approximation of equation (2) is expected to be valid.

The electron dynamics in graphene was computed under a static electric field described as A (t) = −E₀ e _x t, where E₀ is the strength of the applied field, and e _x is the unit vector along the x-direction. After sufficient time has elapsed, the system reaches a nonequilibrium steady state due to the balance between the field-induced excitation and intrinsic relaxation. The electric current in the nonequilibrium steady state can be evaluated as

$$\begin{equation}\boldsymbol{J}\left({E}_{0}\right)=\underset{t\to \infty }{\mathrm{lim}}\enspace \frac{\left(-1\right)}{{\left(2\pi \right)}^{2}}\int \mathrm{d}\boldsymbol{k}\enspace \mathrm{Tr}\left[{\rho }_{\boldsymbol{k}}\left(t\right)\frac{\partial {H}_{\boldsymbol{k}}\left(t\right)}{\partial \boldsymbol{A}\left(t\right)}\right].\end{equation} \tag{ 3 }$$

Due to the circular symmetry of the Hamiltonian in equation (2), the induced current has only the x-component, i.e., J (E₀) = J(E₀) e _x . With this notation, the effective conductivity is here denoted as σ(E₀) = J(E₀)/E₀. It should be noted that, in the weak field limit, the effective conductivity σ(E₀) approaches the linear conductivity i.e., ${\sigma }_{0}={\mathrm{lim}}_{{E}_{0}\enspace \to 0}\sigma \left({E}_{0}\right)$.

Assuming that the THz fields vary slowly in time and the induced electron dynamics is well approximated by the quasi-static description, the THz and dc conductivity of graphene were investigated based on the effective conductivity σ(E₀). This approximation becomes accurate when the THz field frequency ω_THz is sufficiently smaller than the intrinsic relaxation rates such as 1/T₁ and 1/T₂. Figure 1 shows the computed conductivity σ(E₀) as a function of the field strength E₀ for different values of the chemical potential μ. In the relatively weak field regime, the conductivity decreases upon increasing the applied field strength for all the investigated chemical potentials. Since the energy loss of the external field is provided by Joule heating, namely, $J\left({E}_{0}\right){E}_{0}=\sigma \left({E}_{0}\right){E}_{0}^{2}$, the conductivity reduction is directly related to the field-induced transparency of graphene observed in the experiments [19–24]. Using a semi-classical model, this conductivity reduction has been previously explained in terms of an enhancement of the carrier scattering rate [25, 26]. However, the microscopic mechanism of this conductivity reduction has not been yet clarified on the basis of a quantum description. Remarkably, the proposed fully-quantum model can describe the reduction of the electric conductivity without considering the change of the relaxation times, T₁ and T₂. Hence, these results indicate the existence of yet another microscopic mechanism behind the field-induced transparency of graphene, which is different from the change in the scattering rate observed in the classical description.

Zoom In Zoom Out Reset image size

The inset of figure 1 illustrates the nonlinear conductivity σ(E₀) at the charge neutrality point (μ = 0). While the conductivity decreases as the field strength increases up to around E₀ = 1 MV m⁻¹, the conductivity starts to increase almost linearly. In a previous work, an increase in the conductivity via the Landau–Zener tunnel mechanism was suggested [27]. However, this enhancement mechanism has not been investigated in combination with relaxation effects. Therefore, the interplay between the tunneling mechanism and relaxation needs to be considered in order to develop a comprehensive understanding of the nonlinear carrier transport in graphene.

To obtain further insight into these phenomena, the contribution from the intraband current was evaluated. For this purpose, the instantaneous eigenstates were defined as H _k (t)|u_{b, k + A (t)}⟩ = _{b, k + A (t)}|u_{b, k + A (t)}⟩, where b denotes the band index, i.e., valence (b = v) or conduction (b = c) bands. The intraband current is then defined as

$$\begin{equation}{\boldsymbol{J}}^{\text{intra}}\left({E}_{0}\right)=\sum\limits _{b=v,c}\underset{t\to \infty }{\mathrm{lim}}\enspace \frac{\left(-1\right)}{{\left(2\pi \right)}^{2}}\int \mathrm{d}\boldsymbol{k}\enspace \frac{\partial {{\epsilon}}_{b,\boldsymbol{k}+\boldsymbol{A}\left(t\right)}}{\partial \boldsymbol{k}}{n}_{b,\boldsymbol{k}+\boldsymbol{A}\left(t\right)},\end{equation} \tag{ 4 }$$

where the band population n_{b, k + A (t)} is defined as n_{b, k + A (t)} = ⟨u_{b, k + A (t)}|ρ _k (t)|u_{b, k + A (t)}⟩.

Figure 2(a) shows the conductivities computed using the full current J (E₀) and intraband current J ^intra(E₀). The contribution from the intraband current dominates the total conductivity in the whole range of field strengths investigated. Therefore, the band velocity ∂_{b, k} /∂ k and band population n_{b, k} play a central role in the nonlinear conductance of graphene.

Zoom In Zoom Out Reset image size

Due to the circular symmetry of the model, the band velocity exhibits a circular symmetry. Therefore, the intraband current in equation (4) originates from the non-symmetric population distribution in the nonequilibrium steady state. To elucidate the population distribution under field, the total conduction population, N_tot = ∫d k n_{c, k} , and corresponding asymmetric population, ${N}_{\text{asym}}=\int \mathrm{d}\boldsymbol{k}\left\vert {n}_{c,\boldsymbol{k}}-{n}_{c,-\boldsymbol{k}}\right\vert /2$ were evaluated. Figure 2(b) shows the total and asymmetric populations as a function of the field strength E₀ for μ = 0. In the zero field limit, the total population N_tot reaches a finite value due to thermally excited carriers, whereas the asymmetric population N_asym vanishes. In the weak field regime (below E₀ = 1 MV m⁻¹), the asymmetric population shows a faster increase than the total population N_tot upon increasing the field strength. However, N_asym is always smaller than N_tot by definition, resulting in a slowdown of the increase in N_asym. As a result, the field-normalized asymmetric population N_asym/E₀ (see the inset of figure 2(b)) is characterized by a similar reduction to that of the conductivity shown in figure 2(a). Hence, the conductivity reduction can be understood in terms of the suppression of the field-induced asymmetric population. The suppression of the asymmetric distribution can be further understood based on the depletion of the thermal carriers: in the weak field regime, the field-induced asymmetric distribution is formed as a consequence of the displacement of the thermal carriers caused by the applied field. However, once the field strength becomes relatively large, a large proportion of the thermal carriers are already displaced, and a larger asymmetric population cannot be formed. As a result, the effective conductivity is suppressed due to the suppression of the asymmetric distribution. The same mechanism for the doped carriers should explain the conductivity reduction for the finite chemical potential shifts.

In contrast to the weak and moderate field regimes, the total conduction population N_tot increases significantly upon increasing the field strength in the strong field regime (E₀ > 1 MV m⁻¹), as shown in figure 2(b). This can be understood based on the carrier generation through the Landau–Zener tunneling. Consequently, the depletion of the conduction population is eliminated, and the asymmetric population N_asym increases significantly, resulting in the increase of the effective conductivity in the strong field regime. Therefore, the non-monotonic behavior of the effective conductivity of graphene can be understood through the competition between the depletion of the effective carriers and the generation of additional carriers via the Landau–Zener tunneling.

To develop further microscopic insight into the nonlinear conductivity of graphene under strong fields, the conduction population distribution n_{c, k} in the k -space was evaluated. Figure 2(c) shows the conduction band population n_{c, k} at T_e = 300 K in the absence of a static electric field. Here, the Dirac point is set at the origin of the coordinates. As can be seen from figure 2(c), the thermally excited conduction carriers exhibit a circular distribution as a consequence of the circular symmetry of the Hamiltonian. Once a static electric field is applied, the field breaks the circular symmetry of the population distribution, resulting in the intraband current.

To study the field-induced symmetry breaking, the asymmetric population distribution ${\Delta}{n}_{\boldsymbol{k}}^{\text{asym}}=\left\vert {n}_{c,\boldsymbol{k}}-{n}_{c,-\boldsymbol{k}}\right\vert /2$ was evaluated under a static field. Figures 2(d)–(f) show the field-normalized asymmetric distribution ${\Delta}{n}_{\boldsymbol{k}}^{\text{asym}}/{E}_{0}$ for different field strengths. As can be seen from figure 2(d), the induced distribution ${\Delta}{n}_{b,\boldsymbol{k}}^{\text{asym}}$ is close to the circular symmetry distribution in the weak field regime (E₀ = 10⁻³ MV m⁻¹), reflecting the circular symmetry distribution of the thermal carriers. By contrast, once the field strength increases, the induced asymmetric distribution is elongated along the field direction, as shown in figures 2(d) and (e). The elongation can be understood through the intraband motion of the carriers: since the field is aligned along the x-axis, the carriers move along this axis, and a significant asymmetric distribution is thus formed. Once the field strength becomes even larger, the carriers can move across larger distances in the k -space before being scattered. As a result, a significant elongation along the field direction can be formed.

We schematically summarize the microscopic picture behind the nonlinear conductivity of graphene in figure 3. For simplicity, we consider the electron-doped system here. Figure 3(a) shows the electron distribution in the Dirac cone in the absence of the field. Once a weak field is applied, the electron distribution is displaced in the Brillouin zone, as shown in figure 3(b), resulting in the intraband current. If the applied field strength becomes stronger, a large proportion of the carriers are displaced, as shown in figure 3(c), causing the depletion of carriers and suppression of the conductivity. When the field strength becomes even stronger, the Landau–Zener process is activated, and additional carriers are supplied from the valence bands via the excitation, as described in figure 3(c). As a result, the carrier depletion is overcome, and the nonlinear conductivity starts increasing.

Zoom In Zoom Out Reset image size

Having clarified the microscopic physics behind the nonlinear conductivity of graphene, a method is here proposed to control the charge transport in graphene via few-cycle THz pulses. For a THz pulse, the following form of the electric field was considered in the domain −T_d/2 < t < T_d/2,

$$\begin{equation}E\left(t\right)=\frac{{E}_{0}}{{\omega }_{\text{THz}}}\frac{\mathrm{d}}{\mathrm{d}t}\left[-\mathrm{sin}\left({\omega }_{\text{TH}z}t+{\phi }_{\text{CEP}}\right){\mathrm{cos}}^{8}\left(\pi \frac{t}{{T}_{\mathrm{d}}}\right)\right],\end{equation} \tag{ 5 }$$

whereas the field was set to zero outside this domain. Here, ω_THz is the mean frequency of the THz pulse, and ϕ_CEP is the carrier-envelope phase (CEP). The pulse duration T_d is set to 5.38π/ω_THz so that the full width at half maximum of the pulse becomes half of the cycle, i.e., π/ω_THz. Figure 4(a) shows the time profile of the applied electric fields with different CEP values (ϕ_CEP = 0, π/2). As can be seen from the figure, the few-cycle pulses can give rise to inversion symmetry breaking, depending on the value of ϕ_CEP.

Zoom In Zoom Out Reset image size

Based on the inversion symmetry breaking of few-cycle light pulses, light-induced charge transport can be realized [16, 31, 32]. To evaluate the transported charge, the THz fields were assumed to vary sufficiently slowly, so that the induced current could be determined via the instantaneous field strength as $J\left(t\right)=J\left(E\left(t\right)\right)$. The transported charge is then defined by means of the time integration of the current ${Q}_{\mathrm{c}}={L}_{\text{eff}}{\int }_{-\infty }^{\infty }\mathrm{d}t\enspace J\left(E\left(t\right)\right)$, where L_eff is the cross-length that corresponds to the cross-section for bulk systems [31, 32]. It should be noted that the transported charge is proportional to the inverse frequency, 1/ω_THz, because of the instantaneous field dependence of the current, $J\left(t\right)=J\left(E\left(t\right)\right)$. Although the amount of the transported charge Q_c depends on the effective cross-length L_eff and frequency ω_THz, one can introduce a normalized quantity as Q_c ω_THz/L_eff, which is invariant to L_eff and ω_THz.

Figure 4(b) shows the normalized transported-charge Q_c ω_THz/L_eff as a function of the peak field strength E₀ of the THz pulses for ϕ_CEP = 0 and different chemical potential values. By increasing the peak field strength, the charge is first transported in the direction opposite to the peak field direction for all the investigated chemical potentials. A larger amount of charge is transported for a larger chemical potential shift μ. These features of the light-induced charge transport can be understood in terms of the suppression of the nonlinear conductivity σ(E₀) shown in figure 1. This figure shows that the nonlinear conductivity first decreases upon increasing the field strength. Thus, the induced current around the pulse peak timing is suppressed. As a result, the charge transport toward the field peak direction is also suppressed, and the total charge is effectively transported toward the opposite direction. Furthermore, since a larger reduction in the conductivity is observed for the largest chemical potential shift investigated, the resulting charge transport becomes also larger for higher μ values.

The inset of figure 4(b) displays the transported charge for μ = 0. Here, the charge is first transported in the direction opposite to the peak field direction. The transport direction then switches once the peak field strength exceeds 1.25 MV m⁻¹. The change in direction of the transported charge originates from the change in the trend of the nonlinear conductivity: in the relatively weak field regime, the effective conductivity for μ = 0 decreases upon increasing the field strength due to the depletion of the thermal carriers, resulting in the observed negative charge transport driven by the few-cycle THz pulses. By contrast, in the strong field regime, the conductivity increases with the increase of E₀ due to the carrier generation though the Landau–Zener tunneling, thus resulting in the positive charge transport. Hence, the direction of the charge transport can be controlled by tuning the field strength. It should be noted that a similar sign change in charge transport was observed in the infrared field regime and was interpreted in terms of the interference of photo-injected carriers [16–18]. By contrast, the results of the present work are based on the quasi-static picture, and thus the microscopic mechanism behind the present results is different from the one proposed in previous works.

The CEP dependence of the transported charge was investigated further. Figure 4(c) shows the transported charge as a function of ϕ_CEP. Here, μ is set to zero, and the results for the different peak field strength are shown. As can be seen from the figure, the magnitude and sign of the transported charge can be controlled by manipulating the CEP in both the depletion (E₀ = 0.7 MV m⁻¹) and tunneling (E₀ = 5 MV m⁻¹) regimes. The results indicate that ultrafast control of the electric current by means of CEP tuning can be realized also in the THz regime, similar to what has been demonstrated for the infrared regime [16, 31].

In conclusion, the nonlinear conductivity of graphene was investigated using a fully-quantum mechanical model, and it was shown that the effective conductivity exhibits a non-monotonic behavior as a function of the field strength. In the relatively weak field regime, the conductivity decreases upon increasing the field strength, reflecting the depletion of the effective carriers. By contrast, in the strong field regime, the conductivity increases as the field strength due to the carrier generation though the Landau–Zener process. Based on these findings, the possibility of controlling the electric current with few-cycle THz pulses was explored. It was demonstrated that the magnitude and sign of the transported charge can be controlled by manipulating the peak field strength E₀, the CEP ϕ_CEP, and the chemical potential shift μ. These findings may open paths toward achieving ultrafast control of charge transport by light through nonequilibrium and nonlinear electron dynamics in matter.

This work was supported by JSPS KAKENHI Grant Number JP20K14382, the European Research Council (ERC-2015-AdG694097), the Cluster of Excellence 'Advanced Imaging of Matter' (AIM), Grupos Consolidados (IT1249-19) and SFB925 'Light induced dynamics and control of correlated quantum systems'. We thank Enago for the English language review. The Flatiron Institute is a division of the Simons Foundation.

The data generated and/or analysed during the current study are not publicly available for legal/ethical reasons but are available from the corresponding author on reasonable request.