**Nonperturbative Kinetic Description of Electron-Hole Excitations in Graphene in a Time Dependent Electric Field of Arbitrary Polarization**

**Stanislav A. Smolyansky 1,2 , Anatolii D. Panferov <sup>1</sup> , David B. Blaschke 3,4,5,***<sup>∗</sup>* **and Narine T. Gevorgyan <sup>6</sup>**


Received: 16 December 2018; Accepted: 9 April 2019; Published: 16 April 2019

**Abstract:** On the basis of the well-known kinetic description of *e*−*e*<sup>+</sup> vacuum pair creation in strong electromagnetic fields in *D* = 3 + 1 QED we construct a nonperturbative kinetic approach to electron-hole excitations in graphene under the action of strong, time-dependent electric fields. We start from the simplest model of low-energy excitations around the Dirac points in the Brillouin zone. The corresponding kinetic equations are analyzed by nonperturbative analytical and numerical methods that allow to avoid difficulties characteristic for the perturbation theory. We consider different models for external fields acting in both, one and two dimensions. In the latter case we discuss the nonlinear interaction of the orthogonal currents in graphene which plays the role of an active nonlinear medium. In particular, this allows to govern the current in one direction by means of the electric field acting in the orthogonal direction. Investigating the polarization current we detected the existence of high frequency damped oscillations in a constant external electric field. When the electric field is abruptly turned off residual inertial oscillations of the polarization current are obtained. Further nonlinear effects are discussed.

**Keywords:** graphene; dynamic critical phenomena; high-field and nonlinear effects

**PACS:** 81.05.Uw, 64.60.Ht, 73.50.Fq

#### **1. Introduction**

In recent years considerable interest has developed in a nonperturbative, dynamical description of transport phenomena in condensed matter physics inspired by the physics of strong electromagnetic fields [1]. Particular attention was devoted to graphene (see, e.g., [2,3]). In this case there is an obvious similarity with the dynamical Schwinger effect in QED, the creation of electron-positron pairs from the vacuum in strong electromagnetic fields [4–6]. In this context the nonperturbative kinetic approach has proven successful. It is based on the transition to a quasiparticle representation in the presence of an external, quasiclassical electric field facilitated by a time dependent Bogoliubov transformation [7–10]. It would be natural to adopt these methods to specific problems in condensed matter physics and, in particular, to the physics of graphene. Such an adaptation is performed in the present work. The application of these methods allows for advancement to nonperturbative investigations of nonlinear effects in graphene in the presence of strong external electric fields.

We want to give a detailed outline of the contents of the present work. In Section 2, Section 2.1 the basic kinetic equation (KE) for the simplest dynamical model of graphene [2,3,11–14] (a single layer graphene sheet with two Dirac points of the Brillouin zone and absence of the standard scattering mechanism of carriers) is obtained using nonperturbative techniques for the case of a spatially homogeneous, time-dependent external electric field of arbitrary polarization in the graphene plane. The transition to the quasiparticle representation is obtained with the help of a unitary transformation expressed in explicit form [15]. All subsequent consideration is essentially nonperturbative.

The process of electron-hole (*e*-*h*) pair creation in a strong electric field can be considered as a specific field-induced phase transition in a system with broken symmetry [1,10,16]. In Section 2.2 some features of this process are considered in graphene. Section 3 is devoted to the connection of observables such as the quasiparticle number and current densities with the kinetic theory. Here we discuss also the energy conservation law for a system in an external electric field. Here, in particular, an order parameter is introduced which describes the polarization properties of graphene. It is shown that after switching off the external field pulse the order parameter survives and oscillates with momentum dependent amplitude. In other words, the evolution of the order parameter is defined by the entire prehistory of the graphene evolution during the application of the external field. In particular, this effect becomes apparent in the damped oscillations of the residual polarization current on the background of a constant residual conduction current (Section 4). Here it is also shown that the polarization current dominates over the conduction current. This dominance turns out also in calculations of the currents in the framework of the standard perturbation theory. For example, in the Appendix A, we reproduce the well-known results for the polarization and conduction currents in the leading orders of the expansion with respect to *E*/*E*<sup>0</sup> 1, where *E*<sup>0</sup> is the characteristic field (1).

Section 5 contains results of the numerical calculations of the distribution functions of the carriers for electric fields of different magnitude and spectral composition models both for linear and elliptic polarizations. This fact is an important hint that the similar situation is valid also in *D* = 3 + 1 QED, where analogous calculations can be very complicated [17].

In Section 6 we outline the effect of manipulating a weak signal with a current by means of generating active properties of graphene with the help of another (basic) field.

Finally, in Section 7, by analogy with Section 2, we derive the KE in the *D* = 2 + 1 tight binding model of the nearest neighbor interaction [3,11,12]. Also in this case the conduction and polarization currents are obtained. Their detailed investigation will be performed in a separate work.

The conclusions are drawn in Section 8.

We use the metric *<sup>g</sup>μν* <sup>=</sup> diag (1, <sup>−</sup>1, <sup>−</sup>1) and the coordinates *<sup>x</sup><sup>μ</sup>* <sup>=</sup> % *vFt*, *x*1, *x*<sup>2</sup> & . We will proceed from the basic parameters of the model: *a* = 2.46 Å is the lattice spacing, *γ* = 2.7 eV is the hopping energy, and *vF* = 106 m/s is the Fermi velocity. We define a set of scale factors for the physical quantities time (*t*0), momentum (*p*0), and field strength (*E*0) according to

$$t\_0 = \frac{a}{\upsilon\_F}, \ p\_0 = \frac{\hbar}{a'} \quad E\_0 = \frac{\hbar \upsilon\_F}{ea^2}. \tag{1}$$

#### **2. Kinetic Equation**

In this section the basic KE for the description of electron-hole excitations in external, time-dependent electric fields will be derived for the *D* = 2 + 1 QED model of graphene in the framework of a low-energy model (for a tight-binding model, see Section 7 below). Some necessary prerequisites for such a derivation have already been obtained earlier [15] by means of the diagonalization of the initial Hamiltonian of the model. Our approach is based on the consistent usage of the occupation number representation and the adaptation of a method that is well known in *D* = 3 + 1 QED for the description of the creation of an electron-positron plasma from the vacuum in strong fields [7,8,18].

Let us assume the graphene layer is located in the plane % *x*<sup>1</sup> = *x*, *x*<sup>2</sup> = *y* & . A time dependent spatially homogeneous electric field acts in this plane, i.e., the corresponding vector potential in the Hamiltonian gauge is *A<sup>k</sup>* (*t*) = % 0, *A*1(*t*), *A*2(*t*) & . The spatial homogeneity of the electric field can be provided, for example, in the focal spot of two coherent laser beams counter propagating along the axis perpendicular to the graphene layer. It is assumed that the field model is finite, i.e., that the field strength *<sup>E</sup>*(*t*) = <sup>−</sup><sup>1</sup> *c* ˙ *<sup>A</sup>* (*t*) vanishes before switching on and after switching off the laser fields, lim *<sup>t</sup>*→±<sup>∞</sup> *<sup>E</sup>*(*t*) = 0 (the dot above the symbol denotes its time derivative). This is necessary for the correct definition of the in- and out- states of the vacuum with *A*in = *A*(*t* → −∞) and *A*out = *A*(*t* → ∞).

#### *2.1. The Low-Energy Approximation*

The Dirac-type equation for the low-energy excitations in graphene in a time dependent electric field described above is

$$i\hbar \Psi \left( \vec{x}, t \right) = \upsilon\_F \hat{\vec{P}} \vec{\sigma} \Psi \left( \vec{x}, t \right), \tag{2}$$

where *P*ˆ *<sup>k</sup>* = −*ih*¯ ∇*<sup>k</sup>* − (*e*/*c*)*Ak*(*t*) is the quasi-momentum (*k* = 1, 2) and *σ<sup>k</sup>* are the Pauli matrices corresponding to the pseudospin structure of graphene.

The Hamiltonian of the theory,

$$H(t) = \frac{i\hbar}{2} \int d^2 \mathbf{x} \left[ \Psi^\dagger(\vec{x}, t) \Psi(\vec{x}, t) - \Psi^\dagger(\vec{x}, t) \Psi(\vec{x}, t) \right],\tag{3}$$

is the 00 component of the corresponding energy-momentum tensor and it can be transformed with help of the equation of motion (2) to the form

$$H(t) = \upsilon\_F \int d^2 \mathbf{x} \Psi^\dagger(\vec{x}, t) \hat{\vec{P}} \vec{\sigma} \Psi(\vec{x}, t). \tag{4}$$

Here we dropped the spin indices.

The wave function here is a two-component spinor permitting the decomposition

$$\Psi^T(\vec{x},t) = \frac{1}{\left(2\pi\hbar\right)^2} \int d^2p \left(\Psi^{(1)}\_{\vec{p}}\left(t\right), \Psi^{(2)}\_{-\vec{p}}\left(t\right)\right) e^{i\vec{p}\vec{x}/\hbar},\tag{5}$$

which translates the Hamiltonian function (4) to the momentum representation.

For the physical interpretation of the model it is appropriate to go over to the quasiparticle representation, where the Hamiltonian of the theory is diagonal. As it was shown in the work [15], this is achieved with the unitary transformation

$$
\mathcal{U}^\dagger(t)\upsilon\_{\overline{\mathbb{P}}}\overline{\mathcal{P}}\overline{\mathcal{U}}I(t) = \varepsilon(\vec{\mathbb{P}},t)\sigma\_{\overline{\mathbb{S}}} = H\_{\overline{\mathbb{P}}}(t),\tag{6}
$$

and Φ = *U*†Ψ with the unitary matrix [15]

$$
\mathcal{U}(t) = \frac{1}{\sqrt{2}} \begin{pmatrix}
\exp(-i\varkappa/2) & \exp(-i\varkappa/2) \\
\exp(i\varkappa/2) & -\exp(i\varkappa/2)
\end{pmatrix}.
\tag{7}
$$

The function κ is defined by the condition (6) [15], corresponding to tanκ = *P*2/*P*1, where *P<sup>k</sup>* = *<sup>p</sup><sup>k</sup>* <sup>−</sup> (*e*/*c*)*Ak*(*t*). The quasienergy *<sup>ε</sup>*(*p*, *<sup>t</sup>*) in (6) is determined by the dispersion relation in the vicinity of the Dirac points

$$
\varepsilon(\vec{p},t) = \upsilon\_F \sqrt{P^2} = \upsilon\_F \sqrt{(P^1)^2 + (P^2)^2}.\tag{8}
$$

Equation (2) transforms then to the form

$$i\hbar \dot{\Phi} = H\_{\overline{\rho}}(t)\Phi + \frac{1}{2}\lambda \hbar \sigma\_1 \Phi\_\prime \tag{9}$$

where *<sup>H</sup>p*(*t*) is defined by Equation (6) and

$$
\lambda\left(\vec{p},t\right) = \dot{\varkappa} = \frac{\varepsilon v\_F^2 \left[E\_1 P\_2 - E\_2 P\_1\right]}{\varepsilon^2 \left(\vec{p},t\right)}.\tag{10}
$$

Introducing the notation

$$\Phi(\vec{p},t) = \begin{bmatrix} a(\vec{p},t) \\ b^\dagger(-\vec{p},t) \end{bmatrix} \tag{11}$$

the Hamiltonian function (4) can be rewritten in the quasiparticle form

$$\begin{split} H(t) &= \int [dp] \varepsilon(\vec{p}, t) \Phi^{\dagger} \left( \vec{p}, t \right) \sigma\_3 \Phi \left( \vec{p}, t \right) \\ &= \int [dp] \varepsilon(\vec{p}, t) \left[ a^{\dagger} (\vec{p}, t) a(\vec{p}, t) - b(-\vec{p}, t) b^{\dagger} (-\vec{p}, t) \right] \,, \end{split} \tag{12}$$

where the abbreviation [*dp*] = *d*<sup>2</sup> *p*(2*πh*¯)−<sup>2</sup> has been used.

Apparently, the realization of the unitary transformation in the explicit form in both the low-energy and the tight-binding (see below Section 7) models is a result of the fact that these models belong to the class of conformal-invariant field theories (see, e.g., Ref. [6]).

At this stage one can go over to the occupation number representation and replace the amplitudes *a*†(*t*),*a*(*t*) and *b*†(*t*), *b*(*t*) by the corresponding creation and annihilation operators for electrons and holes considered as quasiparticles. These operators are defined on the in-vacuum state |in with vector potential *<sup>A</sup>* in and satisfy the canonical anti-commutation relations

$$\left\{ a(\vec{p},t), a^\dagger(\vec{p}',t) \right\}\_+ = \left\{ b(\vec{p},t), b^\dagger(\vec{p}',t) \right\}\_+ = (2\pi)^2 \delta(\vec{p} - \vec{p}').\tag{13}$$

Other elementary anti-commutators are equal to zero.

From Equations (2), (6) and (11) it follows the equations of motion of the Heisenberg type for the description of the unitary evolution of the creation and annihilation operators, e.g.,

$$\dot{a}(\vec{p},t) = \underbrace{\dot{\hbar}}\_{\cdot} \left[ H(t), a(\vec{p},t) \right] - \frac{\dot{i}}{2} \lambda \left( \vec{p},t \right) b^{+}(-\vec{p},t) = \frac{\dot{i}}{\hbar} \left[ H\_{ltl}(t), a(\vec{p},t) \right],\tag{14}$$

$$\dot{\vec{p}}(\vec{p},t) = \frac{\dot{i}}{\hbar} \left[ H(t), b(-\vec{p}, t) \right] + \frac{\dot{i}}{2} \lambda \left( \vec{p}, t \right) a^+(\vec{p}, t) = \frac{\dot{i}}{\hbar} \left[ H\_{\text{tot}}(t), b(-\vec{p}, t) \right],\tag{15}$$

where the amplitude of the transitions between states with the positive and negative energies of the quasiparticles is defined by Equation (10). From Equations (9), (14) and (15) it follows that evolution of the system is unitary. The Fock space is constructed on the time dependent vacuum state. In Equations (14) and (15) *Htot* = *H* + *Hpol*, where

$$H\_{pol}(t) \quad = \frac{\hbar}{2} \int [dp] \lambda(\vec{p}, t) [a^\dagger(\vec{p}, t) b^\dagger(-\vec{p}, t) - b(-\vec{p}, t) a(\vec{p}, t)] \tag{16}$$

describes the dynamics of vacuum polarization.

Now one can obtain the KE. Let us introduce the distribution functions for the electrons and the holes,

$$f^{\varepsilon}(\vec{p},t) \quad = \quad \langle \text{in} | a^{+}(\vec{p},t)a(\vec{p},t) | \text{in} \rangle,\tag{17}$$

$$f^h(\vec{p}, t) \quad = \quad \langle \text{in} | b^+( -\vec{p}, t) b(-\vec{p}, t) | \text{in} \rangle. \tag{18}$$

The averaging procedure here is carried out under the in-vacuum state |in. Differentiation with respect to time and taking into account Equations (14) and (15) results in

$$\dot{f}^{\varepsilon}(\vec{p},t) = \frac{i\lambda}{2} \left(\vec{p},t\right) \left\{ f^{(+)} (\vec{p},t) - f^{(-)} (\vec{p},t) \right\},\tag{19}$$

where anomalous averages have been introduced

$$f^{(+)}(\vec{p},t) \quad = \quad \langle \text{in} | a^+(\vec{p},t)b^+(-\vec{p},t) | \text{in} \rangle,\tag{20}$$

$$f^{(-)}(\vec{p},t) \quad = \quad \langle \text{in} | b(-\vec{p},t)a(\vec{p},t) | \text{in} \rangle \,. \tag{21}$$

The equations of motion for these functions have the form

$$\dot{f}^{(+)}\left(\vec{p},t\right) = \frac{2i}{\hbar} \varepsilon(\vec{p},t)f^{(+)}\left(\vec{p},t\right) - \frac{i\lambda\left(\vec{p},t\right)}{\ldots\frac{2}{\ddots\cdots\cdots}}[1-2f(\vec{p},t)],\tag{22}$$

$$f^{(-)}\left(\vec{p},t\right) = \frac{-2i}{\hbar}\varepsilon(\vec{p},t)f^{(-)}\left(\vec{p},t\right) + \frac{i\lambda\left(\vec{p},t\right)}{2}[1-2f(\vec{p},t)].\tag{23}$$

Here it was assumed that *f <sup>e</sup>* = *f <sup>h</sup>* = *f* holds as a consequence of the electroneutrality condition.

Let us rewrite Equations (22) and (23) in integral form. Substitution of this result in Equation (19) leads to a KE of non-Markovian type

$$f(\vec{p},t) = \frac{1}{2}\lambda\left(\vec{p},t\right)\int\_{t\_0}^{t}dt'\lambda\left(\vec{p},t'\right)\left[1 - 2f(\vec{p},t')\right]\cos\theta(t,t'),\tag{24}$$

where

$$\theta(t, t') = \frac{2}{\hbar} \int\_{t'}^{t} dt'' \varepsilon(\vec{p}, t'') \tag{25}$$

is the dynamical phase.

In the present work the KE (24) and its reformulation in the form of an equivalent system of ordinary differential equations (ODE), shown below in Equation (27), are considered only for zero initial conditions, *f*<sup>0</sup> = *f*(*t*0) = 0. For the first time a KE of such type was obtained in the works [6,7,19] in *D* = 3 + 1 QED for the description of vacuum creation of electron-positron pairs under the action of a time dependent spatially homogeneous linearly polarized electric field. This method is based on the usage of unitary nonequivalent canonical transformations for the transition to the quasiparticle representation [6]. In the considered situation this approach is applicable and leads to the KE (24) that has the same mathematical structure as in the *D* = 3 + 1 QED case [7,8,18]. However, in the massless *D* = 2 + 1 QED case the transition to the quasiparticle representation is possible in the framework of a unitary transformation [15] (see, e.g., Equation (6)).

An advantage of the unitary approach is also the possibility of a generalization of this method [15] to the case of a two-dimensional electric field with the vector potential *Ak*(*t*)(*k* = 1, 2). Let us remark that the transition from the one-dimensional electric field (linear polarization) to two or three field dimensions (arbitrary polarization) in *D* = 3 + 1 QED is connected with the necessity to take into account a larger number of spin degrees of freedom and is accompanied with a significant increase in the number of necessary KE's [20–22].

The main feature of the KE (24) is the absence of an energy gap in the quasienergy (8). Such kind of models were considered long ago [23] (see also [6]) and have been investigated sufficiently well. In the following this feature will be investigated in the situation when the e-h-system in graphene

is exposed to a time dependent electric field. In the presence of the external field the Dirac points *ε*0(*p*) = 0 are transformed to a family of Dirac lines *LD* which depend parametrically on time,

$$P\_i^D = p\_i^D - \frac{\mathcal{e}}{\mathcal{c}} A\_i(t) = 0, \ i = 1, 2. \tag{26}$$

For the numerical analysis of the KE (24) for different field models it is appropriate to rewrite it in the form of an equivalent system of ODEs [6,8],

$$f = \frac{1}{2}\lambda u, \quad \dot{u} = \lambda \left(1 - 2f\right) - \frac{2\varepsilon}{\hbar}v, \quad \dot{v} = \frac{2\varepsilon}{\hbar}u,\tag{27}$$

with the corresponding initial conditions *<sup>f</sup>*(*t*0) = *<sup>u</sup>*(*t*0) = *<sup>v</sup>*(*t*0) = 0. The auxiliary functions *<sup>u</sup>*(*p*, *<sup>t</sup>*) and *<sup>v</sup>*(*p*, *<sup>t</sup>*) describe polarization effects (Section 3) and can be expressed via the anomalous averages (20) and (21)

$$u = \frac{i}{2} \left[ f^{(+)} - f^{(-)} \right], \\ v = \frac{1}{2} \left[ f^{(+)} + f^{(-)} \right]. \tag{28}$$

A concrete physical interpretation of these functions will be given in Section 3.

For the system of Equation (27) one readily obtains the integral of motion

$$(1 - 2f)^2 + u^2 + v^2 = 1,\tag{29}$$

which is compatible with the zero initial conditions.

There is an approximate nonperturbative solution [24] of the KE (24) which is valid for small occupation numbers, 2 *f* 1 (low density approximation),

$$f\_{\rm LD}(\vec{p},t) \quad = \; f(\vec{p},t) = \frac{1}{2} \int\_{t\_0}^{t} dt' \lambda(\vec{p},t') \int\_{t\_0}^{t'} dt'' \lambda(\vec{p},t'') \cos\theta(t',t''). \tag{30}$$

This integral plays an important role in the formulation of the other nonperturbative approach based on the Markovian approximation (see below).

The polarization function

$$u(\vec{p},t) = \int\_{t\_0}^{t} dt' \lambda(\vec{p},t') [1 - 2f(\vec{p},t')] \cos\theta(t,t')\tag{31}$$

is transformed in the low density approximation to the quadrature formula

$$\mu\_{LD}(\vec{p},t) = \int\_{t\_0}^{t} dt' \lambda(\vec{p},t') \cos\theta(t,t'). \tag{32}$$

From the low density approximation formula (30) and Equation (10) it follows that the distribution function tends to infinity when approaching the Dirac line, *<sup>p</sup>* <sup>→</sup> *pD*(*t*). This indicates also the non applicability of the standard perturbation theory. Thus, close by the lines *LD* an essentially nonperturbative analysis of the KE (24) is required. One such nonlinear approximate solution is obtained in the Markovian approximation based on the neglect of the retardation on the r.h.s. of the KE (24), *<sup>f</sup>*(*p*, *<sup>t</sup>* ) <sup>→</sup> *<sup>f</sup>*(*p*, *<sup>t</sup>*). This results in the quadrature formula

$$f\_M(\vec{p}, t) = \frac{1}{2} \left\{ 1 - \exp\left[ -2J(\vec{p}, t) \right] \right\},\tag{33}$$

which has its analog in the case of D=3+1 QED [24]. From here one can see that the distribution function tends to saturation, *fM*(*p*, *<sup>t</sup>*) <sup>→</sup> 1/2, at *<sup>p</sup>* <sup>→</sup> *pD*(*t*). The polarization function *<sup>u</sup>*(*p*, *<sup>t</sup>*) can be obtained in this approximation on the basis of the first equation of the system (27) and Equation (33)

$$u\_M(\vec{p}, t) = \exp\left[-2f(\vec{p}, t)\right] \int\_0^t dt' \lambda(\vec{p}, t') \cos\theta(t, t'). \tag{34}$$

where *<sup>J</sup>*(*p*, *<sup>t</sup>*) is defined by Equation (30).

#### *2.2. Order Parameter*

By analogy with the standard QED [10], let us introduce the function Φ(*t*) = *u*(*t*) + *iv*(*t*) as an order parameter of the system that describes polarization effects in graphene by means of the anomalous averages (20), (21) and (28) which are characteristic for systems with broken symmetry (e.g. [5,6,16]). We write the corresponding equation of motion

$$
\Phi - \frac{2i\varepsilon}{\hbar} \Phi = \lambda (1 - 2f),
\tag{35}
$$

which follows from Equation (27). The formal solution of this equation with the zero initial condition is

$$\Phi(t) = \int\_{t\_0}^{t} dt' \lambda(t') \left[1 - 2f(t')\right] \exp\left[\frac{2i}{\hbar} \int\_{t'}^{t} d\tau \varepsilon(\tau)\right]. \tag{36}$$

Let us consider now a finite electric field which is switched off at the point of time *t*off, i.e., *<sup>E</sup>*(*<sup>t</sup>* <sup>&</sup>gt; *<sup>t</sup>*off) = 0 and hence according to Equation (10) *<sup>λ</sup>*(*<sup>t</sup>* <sup>&</sup>gt; *<sup>t</sup>*off) = 0. Then, for *<sup>t</sup>* <sup>&</sup>gt; *<sup>t</sup>*off it follows from Equation (36) that the order parameter is different from zero and oscillates with the frequency 2*ε*out/¯*h*, i.e.,

$$\Phi(t > t\_{\rm off}) = \Phi\_{\rm out}(\vec{p}) \exp\left[\frac{2i\varepsilon\_{\rm out}}{\hbar}(t - t\_{\rm off})\right],\tag{37}$$

where the asymptotical value of the quasienergy (8) is equal to

$$
\varepsilon\_{\rm out} = \varepsilon(t \to \infty) = v\_F \sqrt{(\vec{p} - \frac{\mathcal{E}}{c} \vec{A}\_{\rm out})^2} \,\,\,\,\tag{38}
$$

*Ak* out = lim *<sup>t</sup>*→<sup>∞</sup> *<sup>A</sup>k*(*t*). In Equation (37) the momentum dependent amplitude

$$\Phi\_{\rm out}(\vec{p}) = \int\_{t\_0}^{t\_{\rm off}} dt' \lambda(t') \left[1 - 2f(t')\right] \exp\left[\frac{2i}{\hbar} \int\_{\vec{\nu}}^{t\_{\rm off}} d\tau \varepsilon(\tau)\right] \tag{39}$$

is defined by the entire prehistory of the system evolution in a given external field.

The presence of such residual oscillations of the order parameter is a prerequisite for the analogous behavior of the polarization current (see Section 5).

Thus,

$$|\Phi(t > t\_{\rm off})|^2 = |\Phi\_{\rm out}(\vec{p})|^2 = \text{const} \tag{40}$$

after switching off the external field, i.e., the long-lived order is formed.

The amplitude <sup>Φ</sup>out(*p*) of oscillations of the order parameter in the residual state can be defined from the integral of motion (29) by rewriting it in the form

$$(1 - 2f\_{\rm out})^2 + |\Phi\_{\rm out}|^2 = 1.\tag{41}$$

The order parameter Φ(*t*) reflects the role of anomalous averages in the kinetics of the excitation process in graphene, that can be considered as a peculiar field induced phase transition [1,16]. Some other features of this process in graphene will be considered below.

#### **3. Observables**

It is straightforward to write expressions for the pair number density

$$m(t) \quad = \quad \text{N} \int [dp] f(\vec{p}, t). \tag{42}$$

The factor *N* corresponds to number of species (or flavors) of quasiparticles in graphene [3,14,25]: *N* = 4 in the low energy model and *N* = 2 in the tight binding model.

For exact solutions of the ODE system (27) and correct nonperturbative solutions of the type (33) and (34) it follows from the normalization integral (42) that the distribution function is limited everywhere, *<sup>f</sup>*(*p*, *<sup>t</sup>*) <sup>≤</sup> 1. Then both polarization functions *<sup>u</sup>*(*p*, *<sup>t</sup>*) and *<sup>v</sup>*(*p*, *<sup>t</sup>*) are limited also everywhere under the integral of motion (29). This conclusion relates to the neighborhood of the Dirac lines (26) and to the ultraviolet behavior of these functions as well.

The current density consists of two components, the conduction and polarization current densities,

$$j\_k(t) = j\_k^{\text{cond}}(t) + j\_k^{\text{pol}}(t). \tag{43}$$

These currents are defined by the distribution function *<sup>f</sup>*(*p*, *<sup>t</sup>*) and the polarization function *<sup>u</sup>*(*p*, *<sup>t</sup>*), correspondingly [26].

Firstly we consider the currents in the low-energy model. On the basis of the standard definition of the current density [27] (*k* = 1, 2)

$$\begin{array}{rcl}j\_k(t) &=& -\mathfrak{e}\frac{\delta H(t)}{\delta A\_k(t)}\end{array} \tag{44}$$

one can obtain for the theory with the Hamiltonian (4) taking into account the flavor number

$$\left(j\_k(t)\right) = -4ev\_F \int d^2 \mathbf{x} \Psi^\* \left(\vec{x}, t\right) \sigma\_k \Psi \left(\vec{x}, t\right) \,. \tag{45}$$

Going over to the quasiparticle representation with the help of the unitary operator (7) we obtain

$$\dot{j}\_k(t) = 4\epsilon v\_F \int [dp] \Phi^\dagger(\vec{p}, t) \, \mathcal{U}^\dagger(t) \sigma\_k \mathcal{U}(t) \Phi \left(\vec{p}, t\right) \,. \tag{46}$$

Taking into account the spinor (11) and the definition (43), one can separate the conduction and polarization currents,

$$j\_l^{\text{cond}}(t) \quad = \quad 8 \int [dp] v\_q^i(\vec{p}, t) f(\vec{p}, t), \tag{47}$$

$$j\_l^{\text{pol}}(t) \quad = \quad 4 \int [dp] \varepsilon(\vec{p}, t) l\_i(\vec{p}, t) u(\vec{p}, t),$$

where *v<sup>i</sup> <sup>q</sup>*(*p*, *<sup>t</sup>*) = *Pi*/*ε*(*p*, *<sup>t</sup>*) and the vector *li*(*p*, *<sup>t</sup>*) = *δλ*(*p*, *<sup>t</sup>*)/*δE<sup>i</sup>* (*t*) is defined by the components

$$l\_1(\vec{p}, t) = \frac{\varepsilon \upsilon\_F^2 P\_2}{\varepsilon^2}, \quad l\_2(\vec{p}, t) = -\frac{\varepsilon \upsilon\_F^2 P\_1}{\varepsilon^2}. \tag{48}$$

One can see from the system (27) and its nonperturbative solutions (33) and (34), that the polarization effects dominate in the leading approximation for weak fields, *α* = *E*/*E*<sup>0</sup> 1, i.e., *<sup>f</sup>* <sup>∼</sup> *<sup>α</sup>*2, *<sup>u</sup>* <sup>∼</sup> *<sup>α</sup>* and so it follows that

$$|j^{\text{pol}}(t)| \gg |j^{\text{cond}}(t)|\,. \tag{49}$$

This conclusion is supported also by direct numerical calculations.

Let us note, that the conduction and polarization currents are not collinear in the general case. In the case of the linearly polarized electric field collinearity of the currents (47) rebuilds. In order to ascertain this fact, let us consider the situation when the electric field acts along the axis *<sup>x</sup>*1, *E*(*E*1(*t*), 0). Then *<sup>P</sup>*<sup>2</sup> <sup>→</sup> *<sup>p</sup>*<sup>2</sup> and the functions *<sup>f</sup>*(*p*, *<sup>t</sup>*) and *<sup>u</sup>*(*p*, *<sup>t</sup>*) are even and odd under reflection *<sup>p</sup>*<sup>2</sup> → −*p*2, respectively, as it can be seen from the structure of the amplitude (10) and Equation (27). This makes the integrals for *j* cond <sup>2</sup> (*t*) and *j* pol <sup>2</sup> (*t*) in Equation (47) vanish. In order to investigate the theory we calculate the currents (47) in the framework of the perturbation theory in the minimal leading approximation for relatively small external field, see Appendix A.

From Equation (47) it follows that the function *<sup>u</sup>*(*p*, *<sup>t</sup>*) determines the vacuum polarization current. The physical meaning of the other polarization function *<sup>v</sup>*(*p*, *<sup>t</sup>*) is revealed if one considers the total energy density of the quasiparticles including the polarization energy. From Equations (12) and (16) one can obtain *E*tot = *Eq* + *E*pol, where

$$E\_{\vec{q}}(t) \quad = \ \ \ \ \ \delta \int [dp] \epsilon(\vec{p}, t) f(\vec{p}, t) \,, \tag{50}$$

$$E\_{\rm pol}(t) \quad = \ \ \ \ \delta \int [dp] \hbar \lambda(\vec{p}, t) v(\vec{p}, t) \,. \tag{51}$$

Taking the time derivative of the quasiparticle energy *Eq*(*t*) (50) one obtains

$$E\_q(t) = \vec{E}(t)[\vec{j}^{\text{cond}}(t) + \vec{j}^{\text{pol}}(t)] = \vec{E}(t)\vec{j}\_{\text{tot}}(t),\tag{52}$$

where the currents*<sup>j</sup>* cond(*t*) and*<sup>j</sup>* pol(*t*) are defined by Equation (47).

On the other hand, let us write the Maxwell equation for the internal electric field *Ein*(*t*) generated by the motion of the *eh*-plasma,

$$
\dot{E}\_{\text{in}}(t) = -\vec{j}\_{\text{tot}}(t). \tag{53}
$$

On this stage we will imply that the total electric field *Etot*(*t*) is formed by an external field *E*(*t*) and an internal field *<sup>E</sup>*in(*t*), i.e., *<sup>E</sup>*tot <sup>=</sup> *<sup>E</sup>*(*t*) + *E*in. Let us substitute now in Equation (52) the external field *<sup>E</sup>*(*t*) by the total field *E*tot(*t*). Using here Equation (53), we obtain the conservation law of the energy

$$\frac{d}{dt}\left[E\_q(t) + \frac{1}{2}E\_{\rm in}^2\right] = \vec{\mathcal{E}}\,\vec{j}\_{\rm tot}(t). \tag{54}$$

So, the work of external electric fields (r.h.s. of Equation (54)) is distributed between the energy of e-h excitation and the internal electric field.

#### **4. Residual Currents**

Here we consider some nonperturbative effects in graphene which are not sufficiently studied in the standard QED or possess some specific features. We restrict ourselves here to the case of a linearly polarized electric field directed along the axis *x*1.

Let us begin by investigating the residual currents that persist in graphene after the passage of a strong electric field pulse. In the nondissipative model considered here the conduction current discontinues its evolution and remains constant while the polarization current performs damped oscillations. The character of these oscillations and their damping depends on the form of the electric field pulse, as it follows from Equations (37) and (39).

Thus, some oscillating and damped component will be present in the total residual current. In order to calculate it, we will use the formulas for the polarization currents (47) in the low energy model and their analogues in the tight binding model (below in Section 7) with the corresponding polarization function *<sup>u</sup>*out(*t*) for *<sup>t</sup>* > *<sup>t</sup>*off,

$$u^{\rm out}(t) = \text{Re}\,\Phi(t > t\_{\rm out}),\tag{55}$$

where the order parameter <sup>Φ</sup>(*<sup>t</sup>* > *<sup>t</sup>*out) in the out-state is defined by Equations (37)–(39).

According to Equation (37) the frequency of the order parameter is defined by the doubled quasienergy (38). However, while these oscillations are smoothed out upon integration over the momentum space in the polarization current (47), their influence remains quite appreciable, see Figure 1.

**Figure 1. Upper panel:** Supergaussian electric field (56) where *tmax* <sup>=</sup> <sup>1000</sup> *<sup>t</sup>*<sup>0</sup> (2.46 <sup>×</sup> <sup>10</sup>−<sup>13</sup> s) and *Ea* <sup>=</sup> 0.00001 *<sup>E</sup>*<sup>0</sup> (1.088 <sup>×</sup> 103 V/cm). **Lower panel**: The density of the polarization current.

We select the supergaussian model of the electric field

$$E(t) = -\dot{A}(t) = E\_d \exp[-(t - t\_{\text{max}})^4/(2\tau^4)],\tag{56}$$

where *t*max determines the position of the maximum amplitude *Ea* of the field. This choice of the pulse waveform allows to realize abrupt fronts of switching on and off, see the upper panel of Figure 1, and to clearly identify the presence of a alternating polarization current, see the lower panel of Figure 1. This picture demonstrates also dominance of the polarization current.

Another feature of the polarization current becomes apparent in presence of a constant electric field

$$E(t) = E\_{\text{el}} = \text{const}, \ A(t) = -E\_{\text{0}}t. \tag{57}$$

Here the oscillations of the polarization function (Figure 2) transform to damped oscillations of the polarization current (Figure 3). This damping is caused by the monotonic growth of the quasienergy (8) with time at *t* ≥ 0 and, as a consequence, by the decrease of the oscillation amplitude of the polarization current. This mechanism can be traced visually in the Markovian approximation (33).

**Figure 2.** The polarization function *<sup>u</sup>*(*p*1, *<sup>p</sup>*<sup>2</sup> <sup>=</sup> <sup>0</sup>) for the times: *<sup>t</sup>*<sup>1</sup> <sup>=</sup> <sup>50</sup>*t*<sup>0</sup> (1.23 <sup>×</sup> <sup>10</sup>−<sup>14</sup> s, **upper graph**), *<sup>t</sup>*<sup>2</sup> <sup>=</sup> <sup>2050</sup>*t*<sup>0</sup> (5.043 <sup>×</sup> <sup>10</sup>−<sup>13</sup> s, **middle graph**), *<sup>t</sup>*<sup>3</sup> <sup>=</sup> <sup>4050</sup>*t*<sup>0</sup> (9.963 <sup>×</sup> <sup>10</sup>−<sup>13</sup> s, **lower graph**) in a constant electric field (57) *Ea* <sup>=</sup> <sup>5</sup> <sup>×</sup> <sup>10</sup>−<sup>6</sup> *<sup>E</sup>*<sup>0</sup> (5.44 <sup>×</sup> 102 V/cm).

**Figure 3.** The density of the polarization current in a constant field (57) with parameters of Figure 2.

#### **5. Numerical Analysis**

The numerical analysis will be based on the system of ordinary differential Equation (27) rewritten in terms of the corresponding dimensionless values.

We will investigate the response of the system to four electric field models: the constant electric field (57), the Eckart - Sauter field model

$$E(t) = E\_a \cosh^{-2}(t/T), \ A(t) = -E\_a \tanh(t/T), \tag{58}$$

the harmonic function with a constant amplitude

$$E(t) = E\_a \sin(\omega t), \ A(t) = \frac{E\_a}{\omega} - \frac{E\_a}{\omega} \cos(\omega t), \tag{59}$$

where *ω* is the angular frequency, and the "laser field" model [28]

$$\begin{aligned} E(t) &= -E\_d \cos(\omega t) \exp(-t^2/2\tau^2), \\ A(t) &= -\sqrt{\frac{\pi}{8}} E\_d \tau \exp\left(-\sigma^2/2\right) \text{erf}(\frac{t}{\sqrt{2}\tau} - i\frac{\sigma}{\sqrt{2}}) + c.c., \end{aligned} \tag{60}$$

where *σ* = *ωτ*. In all the cases in this section we assume that the electric field is directed along the first coordinate axis.

We start with the most convenient model of the field (58). From Equation (27) follows that the speed of the filling process of the conduction band is determined by the amplitude of the transitions (10). In the denominator of (10) the quasienergy (8) takes zero values on the Dirac line (26). This feature of the amplitude of transitions should be reflected in the behavior and properties of the distribution function. From the form of the evolution of the vector potential (58) it follows that the Dirac line in this case should be represented in the momentum space by a segment with the endpoints determined by *A*(*t* → −∞) and *A*(*t* → ∞) in accordance with the conditions (26).

In Figure 4 we demonstrate the presence of such characteristic features of the distribution function. On the left panel the Dirac line has the end point coordinates *p*<sup>1</sup> ∓ 0.1, *p*<sup>2</sup> = 0.0 while on the right panel the pulse duration is five times larger so that the coordinates of the end points are *p*<sup>1</sup> ∓ 0.5, *p*<sup>2</sup> = 0.0. The Dirac line itself cuts in the distribution function a very thin canyon that is not visible in this figure owing to the selected scale of the numerical calculations.

**Figure 4.** The distribution function in the planar momentum space after the action of the Eckart-Sauter pulse (58). **Left panel**: *Ea* <sup>=</sup> 0.01 *<sup>E</sup>*<sup>0</sup> (1.088 <sup>×</sup> <sup>10</sup><sup>8</sup> V/m) and *<sup>T</sup>* <sup>=</sup> <sup>10</sup> *<sup>t</sup>*<sup>0</sup> (2.46 <sup>×</sup> <sup>10</sup>−<sup>15</sup> s), **Right panel**: *Ea* <sup>=</sup> 0.01 *<sup>E</sup>*<sup>0</sup> (1.088 <sup>×</sup> <sup>10</sup><sup>8</sup> V/m) and *<sup>T</sup>* <sup>=</sup> <sup>50</sup> *<sup>t</sup>*<sup>0</sup> (1.23 <sup>×</sup> <sup>10</sup>−<sup>14</sup> s).

In the next step, we consider the constant field (57) at *t* ≥ 0. The distribution function at the time *t* = 10.0 *t*<sup>0</sup> of the field action is presented on the left panel of Figure 5. Results of the field action with five times longer duration are represented on the right panel of Figure 5.

Another frequently used model is the harmonic electric field (59). The procedure of switching on at *t* = 0 and off at *tm* = 2*πm*/*ω* can be realized with sufficient accuracy in the numerical calculations. The shape of the distribution function and its change in time (*m* = 1, 2, 4 and 10) for the field (59) are presented in Figure 6.

*Particles* **2019**, *2*

**Figure 5.** The distribution function in the constant electric field (57) with *Ea* = 0.01 *E*<sup>0</sup> (1.088 × <sup>10</sup><sup>6</sup> V/cm) at *<sup>t</sup>* <sup>=</sup> <sup>10</sup> *<sup>t</sup>*<sup>0</sup> (2.46 <sup>×</sup> <sup>10</sup>−<sup>15</sup> s) after switching on the field (**left panel**) and at *<sup>t</sup>* <sup>=</sup> <sup>50</sup> *<sup>t</sup>*<sup>0</sup> (1.23 <sup>×</sup> 10−<sup>14</sup> s) after switching on the field (**right panel**).

**Figure 6.** The stages of evolution of the distribution function under the action of the field of type (59) with the number of periods *m* = 1, 2, 4 and 10, respectively.

Finally, let us consider the more realistic model (60) of the "laser pulse". In this case the vector potential and the field strength are changed smoothly at any moment of time and do not bear any problems for the numerical calculations. The shape of the distribution function and its dependence on the pulse width determined by parameter *σ* are presented in Figure 7.

The above results correspond to very short time intervals from *<sup>T</sup>* <sup>=</sup> <sup>10</sup> *<sup>t</sup>*<sup>0</sup> (2.46 <sup>×</sup> <sup>10</sup>−<sup>15</sup> s) to *<sup>T</sup>* <sup>=</sup> <sup>50</sup> *<sup>t</sup>*<sup>0</sup> (1.13 <sup>×</sup> <sup>10</sup>−<sup>14</sup> s) of the electric field action for the models (57) and (58) and for a very high frequency of oscillations *t*0/*T* = 0.1(≈400 THz) for the models (59) and (60). Figure 8 demonstrates the distribution function for the field model (57) and its change for the large time intervals *T* = 406,500 *<sup>t</sup>*<sup>0</sup> (1.0 <sup>×</sup> <sup>10</sup>−<sup>10</sup> s) and *<sup>T</sup>* <sup>=</sup> 1,219,500 *<sup>t</sup>*<sup>0</sup> (3.0 <sup>×</sup> <sup>10</sup>−<sup>10</sup> s) at the field strength *Ea* <sup>=</sup> 9.19 <sup>×</sup> 10−<sup>8</sup> *E*0(10 V/cm). The top row shows images with a linear scale for the color code of the distribution function. This allows to demonstrate that the generated carriers are concentrated in momentum space in a very narrow area close to zero values of the momentum in the direction perpendicular to the direction of the field. The bottom row shows the same distribution function with a logarithmic scale of the color coding. In this case the complicated structure of the distribution function outside the main area of the carrier generation becomes apparent. The main area, however, is absolutely dominant. Other parameters of the electric field can change this picture. This issue requires further research.

**Figure 7.** The residual distribution function in momentum space after the action of the electric field (60) with *Ea* <sup>=</sup> 0.01 *<sup>E</sup>*<sup>0</sup> (1.088 <sup>×</sup> 106 V/cm), *<sup>ω</sup>* <sup>=</sup> <sup>2</sup>*π*0.1. **Left**: *<sup>σ</sup>* <sup>=</sup> 5. **Right**: *<sup>σ</sup>* <sup>=</sup> 10.

**Figure 8.** The residual distribution function in momentum space after the action of the field of type (57) with *Ea* <sup>=</sup> 10 V/cm for a duration of 1.0 <sup>×</sup> <sup>10</sup>−<sup>10</sup> s (**left column**) and 3.0 <sup>×</sup> <sup>10</sup>−<sup>10</sup> s (**right column**). A logarithmic color scale used for the bottom panel.

The set of parameters used here is quite realizable in the experiment. We note that the behavior of the distribution function in momentum space has not undergone fundamental changes in comparison to the ones in Figure 5. Figure 8 demonstrates agreement with the results of the work [29] where the same parameters have been used for the numerical calculations in the framework of another formalism.

The behavior of a quasiparticle plasma under the action of periodic fields (59) and (60) with increasing pulse duration is not trivial. Figure 9 shows the distribution function for the field model (60) with the parameters *<sup>ω</sup>* <sup>=</sup> <sup>2</sup>*<sup>π</sup>* <sup>×</sup> 2.46 <sup>×</sup> <sup>10</sup>−<sup>4</sup> *<sup>t</sup>* −1 <sup>0</sup> (corresponding to 2*<sup>π</sup>* <sup>×</sup> 1.0 THz), *Ea* <sup>=</sup> 9.19 <sup>×</sup> <sup>10</sup>−<sup>6</sup> *<sup>E</sup>*<sup>0</sup> (corresponding to 1000 V/cm) for *σ* = 3, 10, 25, and 50.

**Figure 9.** The residual distribution function under the action of the field of type (60) with the frequency 1.0 THz, amplitude *Ea* = 1000 V/cm and *σ* = 3, 10, 25, 50.

Let us now discuss the analysis of some observable values. We have studied the behavior of the density of carriers (42) in dependence on the amplitude of the electric field for the three models (57), (58) and (60). A summary of the results is presented below.

It should be noted that these results are determined solely by the filling of the conduction band due to the quasiparticle excitation by an external electric field. The presence of thermally excited carriers has not been taken into account as well as the relaxation processes since the considered times are much shorter than the relaxation time. Figure 10 demonstrates that the number of quasiparticles created during the action of the field increases quadratically with increasing electric field strength. The quadratic dependence can be traced quite rigorously in relatively weak fields *Ea* ≤ 0.001 *E*<sup>0</sup> (-10<sup>5</sup> V/cm). The increment of the pair number density is slowing down somewhat with further increase of the electric field. This can be explained by a saturation effect.

The presented values correspond to the pulse of the constant field (57) with duration 20*t*0, the Eckart-Sauter pulse (58) with duration parameter *T* = 10 *t*<sup>0</sup> and a "laser" pulse (60) with a period of the carrier frequency equal to 2*π*/*ω* = 10 *t*0. Such a proximity of the characteristics of the compared field pulses provides very similar values for the surface density of the charge carriers. However, it should be noted that above it has been demonstrated that there is a strong difference between the quasiparticle spectrum in the field of type (60) and the quasiparticle spectrum produced by fields of the type (57) or (58). Nevertheless, the density of carriers and their dependence on the amplitude of the electric field are very similar, see Figure 10.

The dotted line in Figure 10 indicates the approximate level of the thermal carrier density at room temperature. For short pulses their contribution to the total number of carriers will be noticeable only at high electric fields. On the other hand, the spectrum of thermal quasiparticles and quasiparticles generated by the field pulse are different. These differences appear at any electric field strength.

**Figure 10.** The dependence of the carrier density for the electric field models (57), (58) and (60).

Now we come back to the constant field and look at the dynamics of the process of creation of the quasiparticles in the period of the field action. We consider a weak field strength of about 1 V/cm and a field action time of 5 <sup>×</sup> 105*t*<sup>0</sup> (corresponding to 1.23 <sup>×</sup> <sup>10</sup>−<sup>10</sup> s). The left panel of Figure <sup>11</sup> shows the evolution of the distribution function along the direction of the electric field (for *p*<sup>2</sup> 0). The sections of the distribution function along the *p*<sup>1</sup> axis are presented for six time points from *t*<sup>1</sup> = 25,000 *t*<sup>0</sup> to *t*<sup>6</sup> = 500,000 *t*0. This figure shows in more detail the dynamics that we have already seen in Figures 5 and 8. The complete picture is presented in the right panel of Figure 11. It shows the evolution of a slice of the distribution function for the value of *p*<sup>1</sup> = −0.002002 *p*<sup>0</sup> in which the distribution function at the initial period (*t*<sup>1</sup> <sup>=</sup> 2.5 <sup>×</sup> <sup>10</sup>4*t*0, *<sup>t</sup>*<sup>2</sup> <sup>=</sup> <sup>1</sup> <sup>×</sup> 105*t*0) is not large. At the time *<sup>t</sup>*<sup>3</sup> <sup>=</sup> <sup>2</sup> <sup>×</sup> <sup>10</sup>5*t*<sup>0</sup> there is a transition in the state of saturation and then the picture becomes almost stationary.

Figure 12 shows the evolution of three observables under the action of a constant electric field with the same parameters as in Figure 11. The left panel shows the time dependence of the density of the charge carriers (42). The dashed line shows the linear extrapolation of the initial values. The middle panel shows the time dependence of the density of the conduction current in the direction of the field. The dashed line also shows a linear extrapolation. It can be concluded that the creation of charge carriers in a weak constant field proceeds at a constant rate. The energy of the carriers is proportional to their momentum. The number of carriers and the average value of their momentum increase under the influence of the field. As a result, the energy density of the carriers in graphene increases quadratically. This is shown in the right panel of Figure 11 (the dashed line shows the quadratic extrapolation of the initial values).

**Figure 11.** The distribution function in momentum space for the range of time: *<sup>t</sup>*<sup>1</sup> <sup>=</sup> 2.5 <sup>×</sup> <sup>10</sup><sup>4</sup> *<sup>t</sup>*<sup>0</sup> (6.15 <sup>×</sup> <sup>10</sup>−<sup>12</sup> s), *<sup>t</sup>*<sup>2</sup> <sup>=</sup> <sup>1</sup> <sup>×</sup> 105 *<sup>t</sup>*<sup>0</sup> (2.46 <sup>×</sup> <sup>10</sup>−<sup>11</sup> s), *<sup>t</sup>*<sup>3</sup> <sup>=</sup> <sup>2</sup> <sup>×</sup> <sup>10</sup><sup>5</sup> *<sup>t</sup>*<sup>0</sup> (4.92 <sup>×</sup> <sup>10</sup>−<sup>11</sup> s), *<sup>t</sup>*<sup>4</sup> <sup>=</sup> <sup>3</sup> <sup>×</sup> 105 *<sup>t</sup>*<sup>0</sup> (7.38 <sup>×</sup> <sup>10</sup>−<sup>11</sup> s), *<sup>t</sup>*<sup>5</sup> <sup>=</sup> <sup>4</sup> <sup>×</sup> <sup>10</sup><sup>5</sup> *<sup>t</sup>*<sup>0</sup> (9.84 <sup>×</sup> <sup>10</sup>−<sup>11</sup> s), *<sup>t</sup>*<sup>6</sup> <sup>=</sup> <sup>5</sup> <sup>×</sup> <sup>10</sup><sup>5</sup> *<sup>t</sup>*<sup>0</sup> (1.23 <sup>×</sup> <sup>10</sup>−<sup>10</sup> s). The electric field strength *Ea* <sup>=</sup> <sup>1</sup> <sup>×</sup> <sup>10</sup>−<sup>8</sup> *<sup>E</sup>*<sup>0</sup> (1.088 V/cm). The dependence of the distribution function on *<sup>p</sup>*<sup>1</sup> for *<sup>p</sup>*<sup>2</sup> <sup>=</sup> 0 is shown on the left panel while the dependence on *p*<sup>2</sup> for *p*<sup>1</sup> = −0.002002 *p*<sup>0</sup> is shown on the right panel.

**Figure 12.** From left to right: density of charge carriers (42), density of conduction current and energy density the carriers for range of the time 25,000 *<sup>t</sup>*0<sup>−</sup> 500,000 *<sup>t</sup>*<sup>0</sup> (6.15 <sup>×</sup> <sup>10</sup>−12*<sup>s</sup>* <sup>−</sup> 1.23 <sup>×</sup> <sup>10</sup>−10*s*). The constant electric field strength *Ea* <sup>=</sup> <sup>1</sup> <sup>×</sup> <sup>10</sup>−<sup>8</sup> *<sup>E</sup>*<sup>0</sup> (1.088 V/cm).

#### **6. Graphene as Active Medium**

The nonlinear properties of graphene allow to activate it by some basic electric field for driving by the current of another weak signal. Below we will consider an example when a rather strong basic field is aligned with the *x*1—axis while the probe field aligned with the *x*2—axis allows the application of perturbation theory.

The corresponding perturbation theory can be constructed both on the system of Equation (27) and in the framework of the Markovian approximation (33). The latter variant allows to proceed with analytical calculations.

In the framework of such an approximation we can limit ourselves in the weak field approximation to the Markovian solution (33) and (34). The distribution function in this approximation is *fM* ≈ *f* (0) + *f* (1),

$$f^{(0)} = \frac{1}{2}(1 - e^{-f^{(0)}}), \quad f^{(1)} = f^{(1)}e^{-f^{(0)}},\tag{61}$$

where labels (0) and (1) correspond to the "basic" field *A*1(*t*) (nonperturbative solutions) and to the perturbing field *<sup>A</sup>*2(*t*), respectively. The function *<sup>J</sup>*(*p*, *<sup>t</sup>*) is defined by Equation (33), so that

$$f^{(1)}(\vec{p},t) = 8 \int\_{t\_0}^{t} dt' \lambda^{(1)}(\vec{p},t') \int\_{t\_0}^{t'} dt'' \lambda^{(0)}(\vec{p},t'') \cos\theta^{(0)}(t',t''),\tag{62}$$

where according to Equation (10)

$$\lambda^{(0)}(\vec{p},t) = \frac{\varepsilon v\_{\rm F}^2 E\_1(t) p\_2}{2[\varepsilon^{(0)}(\vec{p},t)]^2}, \quad \lambda^{(1)}(\vec{p},t) = \frac{\varepsilon v\_{\rm F}^2 E\_2(t) P\_1}{2[\varepsilon^{(0)}(\vec{p},t)]^2}. \tag{63}$$

In Equation (62) we have neglected in the phase (25) the frequency shift under the influence of the weak field *<sup>A</sup>*2(*t*), i.e., *<sup>θ</sup>* <sup>→</sup> *<sup>θ</sup>*(0). An analogous decomposition of the polarization function (34) leads to the result: *uM <sup>u</sup>*(0) <sup>+</sup> *<sup>u</sup>*(1), where

$$u^{(1)}(\vec{p},t) = J^{(0)}(\vec{p},t) \left\{ \int\_0^t dt' \lambda^{(1)}(\vec{p},t') \cos\theta^{(0)}(t,t') - J^{(1)}(\vec{p},t) \int\_0^t dt' \lambda^{(0)}(\vec{p},t') \cos\theta^{(0)}(t,t') \right\}. \tag{64}$$

Substitution of Equations (61)–(64) into Equation (47) results in the perturbed current calculated in the first order of perturbation theory under the weak field,

$$j\_1^{(1)}(t) = 0, \; j\_2^{(1)}(t) = \int\_{t\_0}^t dt' \sigma(t, t') E\_2(t'), \tag{65}$$

where *σ*(*t*, *t* ) is the linear induced conductivity of graphene controlled by the external field *A*2(*t*),

$$
\sigma(t, t') = -8\varepsilon v\_F^2 \int \frac{[dp]}{\varepsilon^{(0)}(\vec{p}, t)} \left\{ \frac{\delta f^{(1)}(\vec{p}, t)}{\delta E\_2(t')} + f^{(2)}(\vec{p}, t) \frac{ev\_F^2 P\_1^2(t)}{[\varepsilon^{(0)}(\vec{p}, t)]^2} \right. \\
$$

$$
+ P\_1(t) \frac{\delta u^{(1)}(\vec{p}, t)}{\delta E\_2(t')} - u^{(0)}(\vec{p}, t) \frac{ev\_F^2 p\_2}{[\varepsilon^{(0)}(\vec{p}, t)]^2} \}. \tag{66}
$$

Here the first and second groups of terms correspond to the contributions of the conduction and polarization currents.

The dependence of the conductivity (66) on the magnitude and spectral decomposition of the basic field will be considered separately.

#### **7. Tight Binding Model**

It is not difficult to obtain now the analogous KE in the *D* = 2 + 1 tight binding model of the nearest neighbor interaction [3,11,12]. The Hamiltonian function in the momentum representation in this case is

$$H\_{\vec{p}}(t) = \begin{pmatrix} 0 & h\_{\vec{p}}(t) \\ h\_{\vec{p}}^{\*}(t) & 0 \end{pmatrix} = h\_{\vec{p}}^{\prime}(t)\sigma\_1 - h\_{\vec{p}}^{\prime\prime}(t)\sigma\_2. \tag{67}$$

where

$$h\_{\overrightarrow{p}}(t) = h\_{\overrightarrow{p}}'(t) + ih\_{\overrightarrow{p}}'(t) = -\gamma \sum\_{a} \exp\left(\frac{i}{\hbar}P\delta\_{a}'\right) \tag{68}$$

with *γ* ≈ 2.7 eV being the hopping energy, and

$$
\vec{\delta}\_1 = \frac{a}{3} (0, \sqrt{3}), \,\vec{\delta}\_2 = \frac{a}{3} (\pm 3/2, -\sqrt{3}/2) \tag{69}
$$

are the locations of the nearest neighbors, *a* ≈ 3.

An external electric field is introduced here according to the rule *<sup>p</sup>* <sup>→</sup> *<sup>P</sup>* <sup>=</sup> *<sup>p</sup>* <sup>−</sup> *<sup>e</sup>*/*cA* (*t*). Such a method was used in the work [3] in the case of a constant electric field *A*2(*t*) = −*eEt* and resulted immediately in a nonlinear interaction with external field. Such a theory belongs to the class of theories with the highest derivatives.

The Hamiltonian functions (4) and (67) have the same pseudospin structure. Therefore one can follow the way of derivation of KE (24) in the theory with the Hamiltonian function (67). The quasienergy (8) and the amplitude (10) are changed only by the following formal substitutions

$$
\sigma\_F P\_1 \to h\_{\vec{p}}'(t), \; \upsilon\_F P\_2 \to -h\_{\vec{p}}^{\prime\prime}(t). \tag{70}
$$

This results in

$$\varepsilon(\vec{p},t) = \sqrt{h\_{\vec{p}}^\*(t)h\_{\vec{p}}(t)} = |h\_{\vec{p}}(t)|,\tag{71}$$

$$\lambda\left(\vec{p},t\right) = \frac{1}{|h\_{\vec{p}}(t)|^2} \left\{ \dot{h}\_{\vec{p}}''(t)h\_{\vec{p}}'(t) - \dot{h}\_{\vec{p}}'(t)h\_{\vec{p}}''(t) \right\}.\tag{72}$$

An analogous KE can be obtained also for the case of the multilayer graphene model [30].

The conduction and polarization currents have the following form (the flavor number in the given model is equal to 2 [3])

$$j\_k^{\text{cond}}(t) = -4c\gamma \int [dp] f(\vec{p}, t) [F\_k^{(1)}(\vec{p}, t) \cos \chi + F\_k^{(2)}(\vec{p}, t) \sin \chi] \,\tag{73}$$

$$\vec{\eta}\_{k}^{\text{pol}}(t) = -4c\gamma \int [dp] f(\vec{p}, t) \left[ -F\_{k}^{(1)}(\vec{p}, t) \sin \chi + F\_{k}^{(2)}(\vec{p}, t) \cos \chi \right] \,\tag{74}$$

with the vector formfactors

$${}^{1}F\_{k}^{(1)}(\vec{p},t) = \sum\_{a} \delta\_{a}^{(k)} \sin\left(\frac{1}{\hbar} \mathcal{P} \delta\_{a}^{\prime}\right),\tag{75}$$

$$F\_k^{(2)}(\vec{p},t) = \sum\_a \delta\_a^{(k)} \cos\left(\frac{1}{\hbar}\vec{P}\vec{\delta}\_a\right) \tag{76}$$

and *χ* being the angle of the unitary rotation in the matrix of the type (7),

$$
\chi = -h\_{\vec{p}}^{\prime\prime}(t) / h\_{\vec{p}}^{\prime}(t). \tag{77}
$$

Let us rewrite Equations (73) and (74) for the currents to obtain

$$\mathcal{G}\_{k}^{\text{cond}}(t) = -4c\gamma \int [dp] f(\vec{p}, t) \sum\_{a} \delta\_{a}^{(k)} \sin \left( \chi + \frac{1}{\hbar} \vec{P} \vec{\delta}\_{a} \right), \tag{78}$$

$$\delta f\_k^{\text{Col}}(t) = -4c\gamma \int [dp] \mu(\vec{p}, t) \sum\_a \delta\_a^{(k)} \cos \left( \chi + \frac{1}{\hbar} P \delta\_a' \right) \tag{79}$$

as the final result of this section. The numerical evaluation of the currents for particular external field models is delegated to future work.

#### **8. Conclusions**

We have obtained on a nonperturbative basis the KE for describing electron-hole excitations in graphene under the action of a spatially homogeneous time dependent electric field. To this end the analogy with the well-developed case of kinetic theory of vacuum *e*+*e*<sup>−</sup> plasma generation in strong fields [8,10] in *D* = 3 + 1 QED has been used. As a rule, we used the simplest low energy model. However, the used method admits a straightforward generalization to other realistic models of the carrier dynamics as, e.g., the tight binding model of nearest neighbour interaction. The derivation of the KE is based on the transition to the quasiparticle representation [6]. As shown in Section 2, in the *D* = 2 + 1 QED model of graphene such a derivation can be given in an explicit form with the help of a unitary transformation first introduced in the work [15] for the linearly polarized electric field. It is important that the final KE is valid in the general case of an arbitrarily polarized electric field. Some features of the obtained KE are discussed in that Section. In particular, the non applicability of

the standard perturbation theory in vicinity of the Dirac lines has been demonstrated. However, the corresponding approximate analytical and numerical nonperturbative solutions (e.g., the Markovian approximation) of the KE provide a correct description in physical terms in the entire momentum space. In Section 2.2 we consider also some general properties of the evolution of the excited electron-hole plasma that allow to interpret this phenomenon as a specific field induced phase transition [10,16]. An important characteristics of this process is an order parameter that continues to oscillate in the out-state after the external field pulse ceases. The connection of the observables with both, the distribution function and the polarization functions has been discussed in Section 3. The damped oscillations of the residual polarization current on the background of a conduction current were considered in Section 4. The character of these oscillations is related to features of the external field pulse. The damped oscillations of the polarization current in a constant electric field have demonstrated a similar nature. Apparently, these effects are accessible to experimental observation. It can be assumed that similar phenomena occur in *D* = 3 + 1 QED. Here too it was shown that the polarization current dominates over the conduction current. We have performed a systematic numerical investigation based on the KE for the distribution function of quasiparticle excitations and the corresponding observable values for various models of the external electric field.

We have discussed the possibility of using graphene as an active medium excited by the basic electric field to be probed by another signal current ("pump-and-probe").

Finally, we have derived an analogous KE for the tight binding model that is substantially nonperturbative. In the framework of this model we have obtained and discussed the conduction and polarization currents.

A verification of the developed theory was obtained in the work [31], where good agreement with experiment has been shown for the case of a constant electric field [32].

Let us note that both models considered here led KE's of identical form. Moreover, this form invariance is conserved also in the standard QED in the case of a linearly polarized electric field when the spin degrees of freedom are frozen.

**Author Contributions:** Conceptualization, S.A.S. and A.D.P.; Validation, all authors; Data Curation, A.D.P. and N.T.G.; Writing—Original Draft Preparation, S.A.S. and A.D.P.; Writing—Review & Editing, D.B.B. and N.T.G.; Visualization, A.D.P. and N.T.G.; Funding Acquisition, S.A.S. and D.B.B.

**Funding:** This research was supported in part by the "RUDN University Program 5-100", by RFBR according to the research project No 17-02-00375 A, and by the Polish NCN under grant number UMO-2014/15/B/ST2/03752.

**Acknowledgments:** The authors thank V.V. Dmitriev, B. Dora, D.M. Gitman and R. Moessner for useful discussions. D.B. is grateful to Hayk Sarkisyan for inspiring discussions on low-dimensional quantum systems and for the hospitality extended to him at the Russian-Armenian University.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **Appendix A. Perturbation Theory**

In order to demonstrate the effectiveness of the introduced approach, we will reproduce some well known results in the framework of perturbation theory for relatively small external fields, *<sup>E</sup>* < *<sup>E</sup>*0.

We begin with the analysis of currents in the low density approximation (30) and (32) which corresponds to the one-photon excitation mechanism [33]. In the minimal leading approximation *<sup>ε</sup>*(*p*, *<sup>t</sup>*) <sup>→</sup> *<sup>ε</sup>*0(*p*) = *vF*|*p*<sup>|</sup> we have

$$f\_{LD}^{(2)}(\vec{p},t) \quad = \frac{1}{2} \int^t dt' \lambda^{(1)}(\vec{p},t') \text{ ,} u\_{LD}^{(1)}(\vec{p},t') \tag{A1}$$

$$
\mu\_{LD}^{(1)}(\vec{p},t) \quad = \int dt'' \lambda^{(1)}(\vec{p},t'') \cos[\eta \, p(t-t'')],\tag{A2}
$$

with

$$
\lambda^{(1)}(\vec{p},t) = \lambda\_0(\vec{p}) \left[ E\_1(t)P\_2 - E\_2(t)P\_1 \right] \tag{A3}
$$

where *<sup>λ</sup>*0(*p*) = *ev*<sup>2</sup> *F*/2*ε*<sup>2</sup> <sup>0</sup>(*p*) = *<sup>e</sup>*/2*p*<sup>2</sup> and *<sup>η</sup>* <sup>=</sup> <sup>2</sup>*vF*/¯*h*. The upper indices at the functions *<sup>f</sup>* (2) and *<sup>u</sup>*(1) indicate the order of perturbation theory.

The relations (A1) and (A2) indicate the dominant role of polarization effects in the considered approximation, |*j* (1)pol(*t*)||*<sup>j</sup>* (2)cond(*t*)|.

Let us consider the case of linear polarization *E*(*t*) = (*E*(*t*), 0) with arbitrary time dependence of the electric field.

The polarization current in lowest order perturbation theory according to Equation (47) is

$$j\_1^{(1)\text{pol}}(t) \quad = \quad 8\epsilon v\_F \int [dp] u \sin \varkappa\_\prime \tag{A4}$$

$$j\_2^{(1)\text{pol}}(t) \quad = \quad -8\epsilon v\_F \int [dp] u \cos \varkappa\_\prime$$

where κ is defined in explanation to Equation (7),

$$\varkappa = \arctan(P\_2/P\_1) \approx \arctan(p\_2/p\_1),\tag{A5}$$

where the last step corresponds to the leading approximation. From Equation (A5) follows

$$p\sin\varkappa \approx p\_2/p = \sin\Phi, \quad \cos\varkappa \approx p\_1/p = \cos\Phi,\tag{A6}$$

where Φ is the polar angle in the polar representation of the momentum space. Integration over the momentum *p* in the neighborhood of the Dirac points is limited by the cutoff parameter Λ. It is implied that it can be defined by the limits of the validity range of the linear dispersion law *ε*0(*p*) = *vF p*. However, the results obtained below are universal and do not depend on the choice of Λ.

Taking these remarks into account, one can thus write the polarization current after integration over the angle (here t0 → −∞),

$$\begin{split} \left. f\_{1}^{(1)\text{pol}}(t) \right|\_{} &= \left. \frac{\varepsilon^{2} \upsilon\_{F}}{\pi \hbar^{2}} \int\_{-\infty}^{t} dt' E(t') \int\_{0}^{\Lambda} dp \cos[\eta p (t - t')]\_{}, \\\left. f\_{2}^{(1)\text{pol}}(t) \right|\_{} &= \left. 0 \right. \end{split} \tag{A7}$$

Let us now perform a Fourier transformation of the function *E*(*t*) and after that integrate over the momentum *p*,

$$\begin{split} \left(f\_{1}^{(1)\text{pol}}(t)\right) &= \underbrace{\frac{e^{2}}{\pi\hbar}}\_{-\infty} \int d\omega \, E(\omega) \int\_{-\infty}^{t} dt' \frac{\sin[\Lambda \eta(t-t')]}{t-t'} e^{j\omega t'} \\ &= \underbrace{\frac{e^{2}}{\pi\hbar}}\_{-\infty} \int d\omega \, E(\omega) e^{j\omega t} \int\_{-\infty}^{0} \frac{dx}{x} \sin(\gamma x) \cos x. \end{split} \tag{A8}$$

The last integral does not depend on the parameter *γ* = 2*vF*Λ/¯*hω*,

$$\int\_{-\infty}^{0} \frac{dx}{x} \sin \gamma x \cos x = \int\_{-\infty}^{0} \frac{dx}{x} = \frac{\pi}{2} \text{ .} \tag{A9}$$

so that

$$j\_1^{(1)\text{pol}}(t) = \frac{\varepsilon^2}{4\hbar}E(t). \tag{A10}$$

This result does not depend on the choice of the field model.

In order to calculate the conduction current, it is necessary to find the distribution function. To this end we use perturbation theory as a first step and consider the case of a constant electric field (57) switched on at the time *t*<sup>0</sup> = 0. In the leading approximation from Equation (30) follows the known result [15,34]

$$f\_{LD}^{(2)}(\vec{p},t) = \frac{e^2 \hbar^2 E^2 p\_2^2}{4v\_F^2 p^6} \sin^2 \Omega t \,\,\,\,\,\tag{A11}$$

where Ω = *vF p*/¯*h* is the frequency of the vacuum oscillations.

The anisotropic distribution (A11) and the corresponding nonperturbative Markovian distribution (33) have the center symmetry relative to the Dirac point *pi* → −*pi* whereby the conductivity current (see Equation (47)) vanishes.

In order to break this symmetry, it is necessary to go beyond the leading approximation. However, the next correction leads to secular terms that indicate a problem with perturbation theory. For further details on the transport properties of graphene see, e.g., Refs. [35–39].

#### **References**


c 2019 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).
