Nonlinear electric conductivity and THz-induced charge transport in graphene

Based on the quantum master equation approach, the nonlinear electric conductivity of graphene is investigated under static electric fields for various chemical potential shifts. The simulation results show that, as the field strength increases, the effective conductivity is firstly suppressed, reflecting the depletion of effective carriers due to the large displacement in the Brillouin zone caused by the strong field. Then, as the field strength exceeds $1$~MV/m, the effective conductivity increases, overcoming the carrier depletion via the Landau--Zener tunneling process. Based on the nonlinear behavior of the conductivity, the charge transport induced by few-cycle THz pulses is studied to elucidate the ultrafast control of electric current in matter.

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][20][21][22][23][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 * ssato@ccs.tsukuba.ac.jp 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.

II. METHOD
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], where ρ k (t) is the one-body reduced density matrix at each k-point,D [ρ k (t)] 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 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 E(t) = −Ȧ(t) in the dipole approximation. The Fermi velocity v F was set to 1.12 × 10 6 m/s in accordance with a previous ab initio calculation [29]. The relaxation operatorD [ρ k (t)] was constructed using the relaxation time approximation [30].D [ρ k (t)] depends on the longitudinal relaxation time T 1 , transverse relaxation time T 2 , electron temperature T e , and chemical potential µ. T 1 , T 2 , and T e were arXiv:2103.11969v2 [cond-mat.mes-hall] 2 Dec 2021 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 Figs. 2 (c-f)), we consider the electron dynamics only in a narrow region around the Dirac point. Therefore, the linear band approximation of Eq. (2) is expected to be valid. The electron dynamics in graphene was computed under a static electric field described as where E 0 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 Due to the circular symmetry of the Hamiltonian in Eq. (2), the induced current has only the x-component, i.e., J (E 0 ) = J(E 0 )e x . With this notation, the effective conductivity is here denoted as σ(E 0 ) = J(E 0 )/E 0 . It should be noted that, in the weak field limit, the effective conductivity σ(E 0 ) approaches the linear conductivity i.e., σ 0 = lim E0→0 σ(E 0 ).

III. RESULTS
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 0 ). This approximation becomes accurate when the THz field frequency ω THz is sufficiently smaller than the intrinsic relaxation rates such as 1/T 1 and 1/T 2 . Figure 1 shows the computed conductivity σ(E 0 ) as a function of the field strength E 0 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(E 0 )E 0 = σ(E 0 )E 2 0 , the conductivity reduction is directly related to the fieldinduced transparency of graphene observed in the experiments [19][20][21][22][23][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 fullyquantum model can describe the reduction of the electric conductivity without considering the change of the relaxation times, T 1 and T 2 . 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. The inset of Fig. 1 illustrates the nonlinear conductivity σ(E 0 ) at the charge neutrality point (µ = 0). While the conductivity decreases as the field strength increases up to around E 0 = 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 where the band population n b,k+A(t) is defined as Figure 2 (a) shows the conductivities computed using the full current J (E 0 ) and intraband current J intra (E 0 ). 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.
Due to the circular symmetry of the model, the band velocity exhibits a circular symmetry. Therefore, the intraband current in Eq. (4) originates from the nonsymmetric population distribution in the nonequilibrium steady state. To elucidate the population dis- tribution under field, the total conduction population, N tot = dkn c,k , and corresponding asymmetric population, N asym = dk |n c,k − n c,−k | /2 were evaluated. Figure 2 (b) shows the total and asymmetric populations as a function of the field strength E 0 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 0 = 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 0 (see the inset of Fig. 2 (b)) is characterized by a similar reduction to that of the conductivity shown in Fig. 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 fieldinduced 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 0 > 1 MV/m), as shown in Fig. 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 Fig. 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 ∆n asym k = |n c,k − n c,−k | /2 was evaluated under a static field. Figures 2 (d-f) show the field-normalized asymmetric distribution ∆n asym k /E 0 for different field strengths. As can be seen from Fig. 2 (d), the induced distribution ∆n asym b,k is close to the circular symmetry distribution in the weak field regime (E 0 = 10 −3 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 di-rection, as shown in Figs. 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 Fig. 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 Fig. 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 Fig. 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 Fig. 3 (c). As a result, the carrier depletion is overcome, and the nonlinear conductivity starts increasing. 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, 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 . 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(t) = J (E(t)). The transported charge is then defined by means of the time integration of the current Q c = L eff ∞ −∞ dtJ (E(t)), where L eff is the cross-length that corresponds to the crosssection 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(t) = J (E(t)). 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 noamlized transported-charge Q c ω THz /L eff as a function of the peak field strength E 0 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 0 ) shown in Fig. 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 Fig. 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 0 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][17][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 = 0.7 MV/m) and tunneling (E 0 = 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].

IV. CONCLUSION
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 0 , 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.