Kinetic equation approach to graphene in strong external fields

The report presents the results of using the nonperturbative kinetic approach to describe the excitation of plasma oscillations in a graphene monolayer. As examples the constant electric field as well as an electric field of short high-frequency pulses are considered. The dependence of the induced conduction and polarization currents characteristics on the pulse intensity, pulse duration, and polarization is investigated. The characteristics of secondary electromagnetic radiation resulting from the alternating currents is investigated. The nonlinear response to the external electric field characterizes graphene as an active medium. Qualitative agreement is obtained with the existing experimental result of measurements of currents in constant electric fields and radiation from graphene in the case of excitation by means of the infrared and optical pulses.


I. Introduction
Graphene is a unique real material, promising for microelectronics and described by a fairly simple quantum field model (massless D = 2 + 1 QED). These features of graphene and the possibility of the experimental verification of the theoretical predictions explain the great interest in theoretical studies in the electrodynamics of graphene. A specifics here is the lack of analyticity in the coupling constant. This makes the use of nonperturbative approaches relevant here.
One of these directions is associated with the adaptation of general methods of QED for strong fields based on exact solutions of the main QED equations for some simple models of the external electric field [1][2][3]. In graphene, this approach allows a detailed study of the electrodynamics of graphene in the case of a constant [4] and pulsed (Sauter) electric field [5]. For the case of a constant field, satisfactory agreement with experiment was obtained [6][7][8].
Another recently proposed line of research in graphene is based on the use of methods of nonperturbative kinetic theory, which in the simplest case is valid for spatially homogeneous nonstationary electromagnetic fields. Within the framework of standard QED, such an approach was proposed in a number of papers [9][10][11][12] and is currently successfully applied in various problems of QED in strong external fields (for example, [13,14]) in describing the vacuum production and evolution of an electron positron plasma. A promising field of application of this approach is also in QCD, where problems of the creation and evolution of a quark-gluon plasma in the initial stage of the collision of ultrarelativistic heavy ions can be considered (for a review of the early works in this direction, see [15]). This universal approach was recently adapted to the model of single-layer graphene in [16][17][18].
The work consists of two interdependent parts. Sects. II -V contain account of the nonperturbative kinetic theory and electrodynamics of graphene based on the quasiparticle representation. On this foundation the analysis of currents and plasma oscillations is performed in Sects. VI -VIII for the cases of a constant and pulsed external fields. These results are compared with the existing experimental data. Finally, some actual perspectives on forthcoming investigations are discussed in Sects. IX.

II. Statement of the problem
The aim of this work is the construction and subsequent study of the extended kinetic theory of carriers and radiation in single-layer plane graphene excited by the action of an external semiclassical field A µ ex . Internal currents in turn generate internal fields A µ in . Taking into account the possibilities of the theory, quasiclassical fields are assumed to be spatially homogeneous, arbitrarily depending on time and acting in the graphene plane, so that in the Hamiltonian gauge A 0 = 0 the structure of the acting field is as follows: These limitations already mean that the radiation of graphene at the frequencies of plasma oscillations is an isolated problem related to the extension of D2 → D3 electrodynamics. Quasiclassical fields can be strong and are nonperturbatively taken into account, while possible quantum fields are described in the framework of the relevant perturbation theory and can be spatially inhomogeneous and leave the limits of the graphene plane. Thus, the effective fields in graphene will be equal (k = 1, 2) arXiv:2004.03759v1 [physics.optics] 8 Apr 2020 In the nonperturbative part of the problem, the methods and terminology of the kinetic theory based on the quasiparticle representation are used. Supplementing this system with the Maxwell equation for the internal field, we obtain a closed self-consistent system of equations describing at the kinetic level the dynamics of carriers and fields in graphene.
As basic model is considered below the simplest lowenergy model of graphene in the presence of the external field (1) [19][20][21], which describes the excitations in the vicinity of one of the two Dirac points at the boundaries of the Brillouin zone. A generalization of the formalism to a tight-binding model can be found in [17].
We write the equation of motion and the Hamiltonian of the basic model in an effective semiclassical field A k (t) (2) i Ψ ( x, t) = v Fˆ P σΨ( x, t). (3) is the quasimomentum (k = 1, 2), σ k are the Pauli matrices, and v F = 10 6 m/s is the Fermi velocity. The charge of an electron is −e, e > 0. The wave function in (5), (6) is a two-component spinor. The corresponding current density is equal III. Transition to the quasiparticle representation Standard quantum field theory is based on the transition to the representation of second quantization in momentum space, which allows for the usual physical interpretation, where the key element is the concept of particles, real or virtual, which are considered as excitations of the physical vacuum. As a result, in the case of free fields, all necessary quantities become state-additive in momentum space with a certain population determined by the statistics of the fields.
The introduction of interaction with an external field violates the additivity of the observables in the Fock space and complicates the interpretation of the formalism in terms of particles and antiparticles.
The transition to the quasiparticle representation [9] to a certain extent allows us to solve this problem and restore clarity in the description of processes in strong fields.
In graphene, the entire scheme for constructing the kinetic theory of vacuum particle production in the standard QED in the quasiparticle representation is preserved, but it is much simplified when constructing the dynamics when implementing the canonical Bogolyubov transformation. The reason for the simplification is the absence of a mass and the reduction in the number of spatial dimensions. This leads to the fact that the spin degrees of freedom are hidden and degenerate. Their existence is reflected only in the flavor number N s = 2. The inclusion of real spin degrees of freedom requires a D2 → D3 generalization of the theory.
We suppose that graphene is located in a region bounded by a square with side L and go to the momentum representation (V = L 2 ), In this representation, we write the basic Hamiltonian (4) where P = p − e c A(t) is the quasi-momentum, and the equation of motion (3) is The diagonalization of the Hamiltonian (8) can be done explicitly using the unitary transformation Ψ = U Φ with the matrix [22] The parameter tan κ = P 2 /P 1 is fixed here by the equality where ε( p, t) is the quasienergy According to the selected low-energy model, the dispersion law (12) is valid in the vicinity of the Dirac point p 2 = 0 at the boundaries of the Brillouin zone.
The equation of motion in the quasiparticle representation will have the form where the function describes the transition between states with positive and negative energies and can be found from the equality In the formula (14) E k (t) = − 1 cȦ k (t) is the electric field strength.
If we now write the spinor Φ( p, t) taking into account expansion (7) in the form then the Hamiltonian of the system in the new representation will be equal to: The spinor (16) can be associated with the field operator are the unit spinors. We now turn to the representation of occupation numbers and define the creation and annihilation operators of electrons and holes over the instantaneous vacuum state |t , demanding for these operators the standard anticommutation relations, The remaining elementary anticommutators are equal to zero.
The equations of motion (13) can be written using (16) in terms of amplitudes or in the form of Heisenberg-type equations with respect to the operator ϑ where is the Hamiltonian of the fermion subsystem, including the Hamiltonian H(t) (17) and the Hamiltonian describing the effects of vacuum polarization under the influence of the field.

IV. Basic kinetic equation
Below, in the quasiparticle representation, under kinetic equation (KE) will be understood the closed integro-differential equations for the distribution functions of electrons and holes f e ( p, t) = in|a + ( p, t)a( p, t)|in , Subsequently, these functions are considered equal, by virtue of the assumption of electroneutrality of the system at each moment of time. The averaging in (25), (26) is performed over the in-vacuum state.
To get a closed system of KE, we differentiate f ( p, t) with respect to time and use the equations of motion (21) where anomalous averages are introduced by Differentiating them in time, we obtaiṅ As a result of the integration of these equations over time and substitution in (28), we obtain the desired KE of non-Markovian typė where the dynamic phase is introduced by The equation (33) should be supplemented with an initial condition compatible with the requirement of electroneutrality (27). The simplest option f 0 ( p) = f ( p, t = t 0 ) = 0 corresponds to the absence of excitations at the initial moment of time t 0 . For the first time, KE of this type were obtained in independent works [9,10,11] in the framework of standard QED. The specifics of the quantum field system under consideration consists in the ability to perform canonical Bogolyubov transformations in an explicit form (Section 2). Note that the KE (33) is valid in the case of an arbitrary polarization of the external field (1). The integro-differential KE (33) can be written in the form of an equivalent system of ordinary differential equationṡ with appropriate initial conditions The auxiliary functions u( p, t) and v( p, t) in the ODE system (35) describe vacuum polarization effects (Section 4) and can be written in terms of anomalous averages (29), (30) The KE system (35) has an integral of motion compatible with the initial conditions (36). An interesting property of the KE (33) (or its equivalent system (35)) is the preservation of the non-negativity of the distribution function f ( p, t) ≥ 0 in the entire domain of its definition.

V. Macroscopic averages
Since, due to the uncertainty relation, the description of the system in terms of time-dependent quasienergy (16) is conditional and acquires physical meaning only in the asymptotic region t → ∞, the distribution function itself is rather conventional. This feature of the quasiparticle approach was noted in the literature (for example, [23]). On the other hand, it is understood that under certain conditions (for example, in the case of systems located in a limited region of space, V < ∞), the evolution of a macroscopic quantum field system in timedependent external conditions can be controlled at any time in various ways (using radiation in to the external region or various responses to weak external probes). Such, for example, are experiments on graphene samples in the optical excitation range [24,25] or planned experiments to detect an e − e + plasma in the focal spot of computer-propagating powerful laser beams [26][27][28].
In such a situation, macroscopic averages obtained by integrating the dynamic characteristics of quasiparticles over the momentum space, statistically weighted with the distribution function f ( p, t) playing the role of time-dependent functions -intermediaries between the dynamic characteristics of quasiparticles and physical observables.
The density of excited e − h pairs is the simplest quantity of this type.
where N f = N S N D = 4 is the total number of flavors in the model under consideration. From the full Hamiltonian of the fermionic subsystem (23) follows the expression for the total energy density of the electron -hole subsystem, which consists of the sum of the quasiparticle and polarization parts where where v( p, t) is the vacuum polarization function (37). The current density can also be represented as the sum of the conduction and polarization currents where To write these currents in terms of the functions f, u, v, we use the formula (6) and perform a unitary transition to the quasiparticle representation there (Section 2), Using formulas (10), (16) and the definition (37) here, we obtain the expression for the conduction current density as a result where is the propagation velocity of the quasiparticle excitation.
In the derivation of formula (46), the vacuum unit was omitted 2f − 1 → 2f . The polarization current density is where the vector of "conjugate velocity"ṽ k Using this definition, one can write the amplitude (14) in the following form: Thus, the functions u( p, t) and v( p, t) determine the energy and current due to the polarization of the medium. We now obtain the law of conservation of energy in the quasiparticle subsystem. Differentiating in time the density of quasiparticle energy (41) and using the formulas (46)-(49) for the current density and relation (50), we obtainĖ where J(t) is the total current density (43). This relation allows a different formulation. We write the Maxwell equation for the internal plasma fielḋ and write in Eq. (51) the total electric field as a sum As a result, we obtain Similarly, relation (42) implies the relation describing the energy balance of vacuum polarizatioṅ In the next Sections the electrodynamics of graphene will be considered in detail in the cases of a constant electric field (Sect. VII) and an alternating field (Sect. VI) as well as some features of the inner plasma field.

VI. Constant field
We begin the demonstration of the opportunities of the developed approach with the case of a constant external electric field: Constant is both the absolute value of the field strength and its direction. Further, the first axis of the used twodimensional coordinate system is associated with this direction. In this case, it is convenient to compare with the results of using the Landauer approximation (or the Landauer-Datta-Lundstrom (LDL) model) [29][30][31] for calculating the characteristics of the current arising in graphene samples with the results of the presented approach. Some nontriviality is that quite different methods and different spaces of description are used here: the t -representation in the considered kinetic method and the x -representation in the LDL -tunneling approach. A constant field does not have any characteristic time scales. Those can only be sought in the characteristics of the material itself or of the simulated sample. For a material (graphene), the natural unit of the time scale is the ratio of the lattice constant to the Fermi velocity a/v F ≈ 2.46 × 10 −16 s. It makes no sense to consider processes on a smaller time scale. The upper boundary of the time scale is determined by the condition of applicability of the assumption of spatial homogeneity of the system. It depends on the characteristic size of the sample and can be determined by the ratio L/v F .
The system of equations (35) allows only a numerical investigation. For each point of the momentum space {p 1 , p 2 }, it is solved independently. The region in which the distribution function is localized and the required density of its coverage by the computational grid at the nodes of which solutions for (35) will be sought is determined by the necessary accuracy in calculating the integral characteristics. This is realized by a sequential iterative procedure with stepwise control of the accuracy of the results obtained. Fig. 1 shows the form of the distribution function for two consecutive points of time at a field strength of 0.1 V/µm. The "natural" value /a is used hereinafter as a unit value for the p 1 , p 2 (a is the lattice constant). Such kind of accumulative behavior of the distribution function is in agreement with the results of [32] and [33]. The latter work based on the Greens function method takes into account the dissipative mechanism of inelastic scattering of optical phonons.
The distribution function of the excitation on the initial stage is formed and reaches the maximum value f max = 1 very rapidly. at invariable Gauss-like distribution on the transversal momentum p 2 . The halfwidth of the p 2 -distribution is defined by the field strength E (Fig. 2). As a result, the number density of carriers grows proportionally to the action time T of the field, n(T ) ∝ T , as in the case of standard QED (see, e.g., [9]).
In the LDL approach, the problem of calculating the current density through the sample of a material of finite width L a = 0.246 nm limited by two parallel electrodes is solved by calculating the transmission probability T ( , p 2 , V ) in the presence of a given potential difference V (for dispersion of graphene = v F (p 1 ) 2 + (p 2 ) 2 ) [6]: Therefore, the momentum component p 2 does not change its value during tunneling. In spite of the fact that graphene is a gapless semiconductor, the presence of finite conserved values p 2 leads to the appearance of an energy gap with the width ∆ = 2v F p 2 . In the considered case of a vacuum the initial state temperature and chemical potential are equal to zero, the Fermi energy F = 0 and there are no free carriers in the interelectrode space. Under these conditions the process of carrier transmission can proceed only by Zehner-Klein tunneling. The probability of this tunneling in the WKB approximation is (taking into account the relationship of and p 2 ): In this case from (57) and (58) one obtains [6]: First of all, let us note that the distribution (58) of carriers over the transverse momentum p 2 accurately re-produces the result (Fig. 2, right panel) obtained in the t-representation.
The basic problem is that the conduction current density (46) in the absence of dissipation and spatial boundaries depends on the time after switching on the field and increases continuously, which, obviously, is not observed in the experiment, where each value of the potential difference corresponds to its steady-state current density. So, under the conditions of the real experiment [6] it is necessary to take into account the presence of the electrodes, which limit the lifetime of the carriers. From this point of view the process of carrier generation must permanently continue throughout the entire measurement process and be uniform in the area of the sample. Knowing the strong anisotropy of the carrier spectrum, it is possible to assume in the first approximation that all of them move towards the electrode of the corresponding polarity with the velocity v F . In this case, the average lifetime of carriers is τ = L/2v F 0.5 × 10 −12 s for L = 1.0 µm. At the end of this time after switching on the field, the rate of generation of carriers becomes equal to the rate of their escape through the electrodes. The steady-state current values will be constant and can be calculated in both approaches (Tab.1).

V j in x-representation j in t-representation [V]
[   The results of these calculations coincide with high degree of accuracy in the considered range of parameters. Indirectly this shows a good coincidence with experiment [6]. The given list of values strictly corresponds to the law I ∼ V 3/2 ∼ E 3/2 . This dependence is corroborated by other sources [4,5].
The kinetic method presented here can be a valid outside the framework of the applicability of the xrepresentation, for example, in alternating electric fields in the region of sufficiently high frequencies. An indirect confirmation of this is the difference in the results for a thinner (tenfold) sample with a correspondingly shorter carrier lifetime (Fig. 3). In this case the difference reaches almost 10%. Both approaches also predict superlinear growth of the current density. But the exponent is already 1.70 for both tand x -representations.
Above, we limited ourselves to considering only the conduction current in the situation of a stationary process, when the number of carriers born compensates for their departure from the sample. In the kinetic approach, the behavior of the conduction current (46) and polarization current is fundamentally different from conduction current behavior and has the form of damped oscillations asymptotically reaching a constant negative value. Nevertheless, the contribution of the conduction current becomes dominant quite quickly.

VII. Short high frequency pulses
Next, we consider the field model of pulses with a cyclic carrier frequency ω and a Gaussian envelope of width τ of the form (60) The phase shift ϕ sets the position of the absolute maximum of the field relative to the maximum value of the envelope (which will be associated with the time t = 0). The spatial direction of the field is constant. Since no restrictions are imposed on the carrier frequency and pulse width, this model can be used in a very wide range up to optical radiation.
We choose the characteristics of the field so that they correspond to the parameters of the experiments described in [34]. These are very short pulses with ω = 2π × 2 THz carrier frequency and duration parameter τ = 3/ω ≈ 2.4 × 10 −13 s, which almost coincides with the constant-field action time considered above. But since in this case the field turns on and off relatively slowly, to accurately reproduce the behavior of the model, we will consider a several times longer interval. The electric field amplitude is E a = 3.0 V/µm, which is an order of magnitude greater than that considered above for a constant field. To accurately reproduce the parameters, we used a nonzero value of the carrier phase shift ϕ = 0.85 π. An explicit form of the dependence of the field on time is shown in Fig. 5. The field is formed by a linearly polarized electromagnetic wave and is assumed to be directed along the x-axis. The field strengths are expressed as The distribution function is localized in the vicinity of the Dirac point. The distribution of populated states is much more complicated than in the case of a constant field. Thus reproduction of integral characteristics in this case may be more time consuming. The figure is based on the results of solving kinetic equations on a maximally simplified 222 × 46 grid. The carrier density (39) evolution reproduced from these data for a time interval from −6.0 × 10 −13 s to 6.0 × 10 −13 s is shown in the following Starting from zero, the carrier density is quite difficult to evolve and reaches a final constant value, which is the residual carrier density. A short field pulse generates a number of carriers, which, in the absence of dissipation, continue to exist even after the field is turned off.
We consider one more variant of a field of the form (60). This will be a relatively long infrared pulse, similar in parameters to those considered in [25]. We define its cyclic frequency as ω = 2π×96.7 THz, which corresponds to a wavelength of 3.1 µm and a photon energy of 0.4 eV. The duration parameter is τ ≈ 26/ω ≈ 4.28 × 10 −14 s. The phase shift was chosen equal to zero. The maximum electric field strength was determined based on the declared energy flux density in the focal spot of 7 GW/cm 2 and reaches ≈ 2.2 × 10 2 V/µm and exceeds by almost two orders of magnitude the corresponding value from the previous example. The impulse with these characteristics has the form shown in Figure 8. The form of the distribution function at the final stage of the action of the external field is shown in Figure 9. At the stage of increasing field amplitude, an increase in the density of carriers with small oscillations is observed. At the stage of turning off the field, a drop to a certain residual value takes place.
Conduction and polarization currents for a pulsed alternating field with a cyclic frequency ω = 2π × 2 THz (Fig. 5) are shown in Fig. 11. For the conduction cur- demonstrates a very good reproduction of the general nature of the acting field: it is alternating, the period of polarity change obviously coincides with the period of the carrier of the simulated pulse. The amplitude of the conduction current increases and decreases along with an increase and decrease in the amplitude of the acting field. Obviously, there is a phase shift between it and the external field, because at the time t = 0, when the field reaches its maximum value (and this is the absolute maximum), the conductivity current is close to zero. The behavior of the polarization current is fundamentally more complex. Over the entire time interval presented, in which the conduction current has time to increase by several of orders of magnitude and then also to decrease its amplitude, the polarization current does not show obvious signs of a regular change in its amplitude, but only irregular changes in the current value. In the initial section of the time interval presented, one can try to associate its behavior with the current values of the external field, but then a complex interference pattern is observed, which is characteristic for superposition of oscillations with different frequencies. We can draw an analogy with the slow relaxation of the polarization current already after the completion of the action of a short unipolar pulse, which was noted in [35].

VIII. Plasma field
The appearance of currents, by virtue of Maxwell's equation (52), leads to the generation of a plasma field, which will add to the external field and in the general case will affect the process of further evolution of the system. We use the well-known solution to the problem of the field of a time-dependent current on an infinite plane (see, for example, Ref. [36]): Here µ 0 is the magnetic susceptibility of vacuum, c is the speed of light and z is the distance of the observation point from the graphene plane. It follows from this expression that, in the case under consideration, the plasma field is determined by the current density. From the expression (61) it follows that the alternating current will generate radiation that carries information about its characteristics. When considering a real sample with a characteristic finite size L, all of the above will be true only for a distance z from its surface that satisfies the condition z L. One of the reasons for the continued interest in graphene is the observation of a nonlinear response in it under the action of pulsed fields of the form (61) [25,34]. Does the presented model based on the quantum kinetic equation reflect these properties of graphene? To answer this question, we study the spectral composition of internal currents shown in Fig. 11 and 12.
As a tool for studying time series, we use the discrete Fourier transform. The figures (periodograms) below show the squares of the moduli of the coefficients of such a transformation and reflect the contributions of various frequencies to the total field energy. The horizontal axis of the frequencies is calibrated in T Hz, the vertical axis has a logarithmic scale. Since the absolute values of the periodogram depend, inter alia, on the length of the sample that we have, only relative values within one periodogram have a well-defined meaning. The spectrum of the acting field plays the role of a reference point.
But in order to have an initial reference point, we first present in Fig. 13 the result obtained from a number of discrete values of the external electric field shown in Fig. 5. Since the pulse is very short, its spectrum is greatly broadened. Nevertheless, the maximum for the value of 2 THz is quite visible, and outside it there is a smooth and fairly monotonic decrease in the intensity of the spectral components. On this background are shows the results of exactly the same processing of a sequence of conduction and polarization currents. The behavior of the spectrum of the conduction current is similar to the behavior of the spectrum of the acting field. But the curve for the polarization current clearly contains local maxima for frequencies of 2, 6, 10, and 14 THz, which should be interpreted as the appearance of odd harmonics in the spectrum of the internal field with multiplicities of 3, 5, and 7 with respect to the fundamental frequency. If at the fundamental frequency the contribution to the field energy from the polarization current is about three orders of magnitude less then from the conduction current, at the third harmonic they are comparable, and at the fifth  5) and generated by it currents (Fig. 11).
and seventh contributions of the polarization current are almost an order of magnitude larger.
Let us now see how this looks in the case of a longer pulse, the spectrum of which is more localized. The spectrum of the pulse itself is shown in Fig. 14 FIG. 14: Periodograms for a pulse of a periodic electric field (Fig. 8) and generated by it currents (Fig. 12).
currents are also shown in this figure. In this case, the spectrum of both currents clearly shows the odd harmonics with numbers 3, 5, 7, 9, ... In the above figure, the upper frequency limit is 1000 THz, but if we raise it, we will see harmonics of a higher order. At the carrier frequency, the contribution of the conduction current also dominates. But already at the next, third, the contribution of the polarization current is about one and a half orders of magnitude greater. With increasing frequency, the contributions of the conduction and polarization currents become almost the same. The appearance of odd harmonics is confirmed in experiments [25,34]. The obtained results are in qualitative agreement with the existing experimental data, so that the suggested kinetic theory is surely verified.

IX. Conclusion
In this work we have considered the application of the nonperturbative kinetic equation approach to describe the excitation of plasma oscillations in a graphene monolayer. As examples the constant electric field as well as an electric field of short high-frequency pulses have been considered. The dependence of the induced conduction and polarization current characteristics on the pulse intensity, pulse duration, and polarization have been investigated numerically for these examples. The characteristics of secondary electromagnetic radiation resulting from the alternating currents has been studied and a nonlinear response to the external electric field has been found which characterizes graphene as an active medium.
A perspective direction of development of the kinetic theory of graphene is a generalization to the case of additionally accounting for the interaction with the quantized electromagnetic field. One can proceed from the analogy with the works [37,38], where such a generalization was performed in the standard QED on the basis of the Bogolyubov-Born-Green-Kirkwood-Yvon chain of equations in the single photon approximation. Preliminary investigations of radiation on this basis met large difficulties [39]. Analogous research in graphene would allow to compare the quasiclassical radiation (Sect. VIII) with the quantum one and to understand deeper the situation in standard QED. Considerable interest represents a study of cascade processes in graphene (e.g., [40,41] in standard QED) and also of spin phenomena.