Radiation Problems Accompanying Carrier Production by an Electric Field in the Graphene

Abstract

A number of physical processes that occur in a flat one-dimensional graphene structure under the action of strong time-dependent electric fields are considered. It is assumed that the Dirac model can be applied to the graphene as a subsystem of the general system under consideration, which includes an interaction with quantized electromagnetic field. The Dirac model itself in the external electromagnetic field (in particular, the behavior of charged carriers) is treated nonperturbatively with respect to this field within the framework of strong-field QED with unstable vacuum. This treatment is combined with a kinetic description of the radiation of photons from the electron-hole plasma created from the vacuum under the action of the electric field. An interaction with quantized electromagnetic field is described perturbatively. A significant development of the kinetic equation formalism is presented. A number of specific results are derived in the course of analytical and numerical study of the equations. We believe that some of predicted effects and properties of considered processes may be verified experimentally. Among these effects, it should be noted a characteristic spectral composition anisotropy of the quantum radiation and a possible presence of even harmonics of the external field in the latter radiation.

1. Introduction

Particle creation from the vacuum by strong electromagnetic and gravitational fields is a remarkable quantum effect predicted first in References [1,2,3,4,5] and studied then in References [6,7,8,9,10] in the framework of the relativistic quantum mechanics. Its exhaustive explanations and consistent nonperturbative methods of investigation became possible in the framework of quantum field theory (QFT). QFT with external backgrounds is, to a certain extent, an appropriate model for this purpose. In the framework of such a model, the particle creation is related to a violation of the vacuum stability. Backgrounds (external fields) that violate the vacuum stability are electric-like fields that are able to produce nonzero work when interacting with charged particles. Different approaches for calculating the effect were proposed and realized, depending on the structure of such backgrounds. From a quantum mechanical point of view, the most clear formulation of the problem of particle production from the vacuum by external fields is possible for time-dependent homogeneous external electric fields that are switched on and off at infinitely remote times $t\to \pm \infty$, respectively. Such a kind of external fields are called t-electric potential steps (t-steps). Scattering, particle creation from the vacuum and particle annihilation by the t-steps were considered in the framework of the relativistic quantum mechanics, see References [6,7,8,9,10], a more complete list of relevant publications can be found in References [11,12]. A general nonperturbative with respect to the external background formulation of QED was developed in Refs. [13,14,15]. In contrast to the t-electric potential steps, there are many physically interesting situations where the external backgrounds are constant (time-independent) but spatially inhomogeneous, for example, concentrated in restricted space areas. The simplest type of such backgrounds are so-called x-electric potential steps (x-steps), in which the field is only inhomogeneous in one space direction and it represents a spatial-like potential step for charged particles. The x-steps can also create particles from the vacuum, the Klein paradox is closely related to this process [1,2,3,4]. Important calculations of the particle creation by x-steps in the framework of the relativistic quantum mechanics were presented in References [6,7,8] and later developed in Ref. [16]. A general nonperturbative with respect to such external background formulation of QED was formulated in Ref. [17]. It is based on the existence of exact solutions of the Dirac or Klein-Gordon equation (wave equations, in what follows) with corresponding external fields. When such solutions can be found and all of the calculations can be analytically done, we refer these examples as exactly solvable cases. Until now, there are known only few exactly solvable cases that are related to t-steps and to x-steps. In the case of t-steps, these are particle creations in a constant uniform electric field [5,6,7,8], in an adiabatic electric field [9,10], in the so-called T-constant electric field [18,19], in a periodic alternating in time electric field [9,10,20], in an exponentially decaying electric field [21], in an exponentially growing and decaying electric fields [22,23] (see Ref. [24] for the review), in a composite electric field [25], and in an inverse-square electric field [26]. In the case of x-steps, these are particle creation in the Sauter electric field [17], in the so-called L-constant electric field [27], and in the inhomogeneous exponential peak field [28] and inverse-square electric field [26,29].

Until recently, problems that were related to particle creation from the vacuum had mostly theoretical interest. This is related to the fact that the vacuum instability can be observed only in extremely strong external electric fields of the magnitude of $E_c=m^2/e\simeq {10}^{16}$ V/cm ($E_c$ is the Schwinger’s critical field). However, recent technological advances in laser physics suggest that lasers, such as those planned for the Extreme Light Infrastructure project (ELI), may be able to reach the nonperturbative regime of pair production in the near future (see review [30]). Moreover, the situation has changed completely in recent years regarding applications to condensed matter physics: particle creation became an observable effect in graphene physics, an area that is currently under intense development [31,32,33,34]. Briefly, this is explained by two facts: first, the low-energy electronic excitations in the graphene monolayer in the presence of an external electromagnetic field can be described by the Dirac model [35], namely, by a $2+1$ quantized Dirac field in such a background (here, dispersion surfaces are the so-called Dirac cones); and, second, the gap between the upper and lower branches in the corresponding Dirac particle spectra is very small, so that the particle creation effect turns out to be dominant (under certain conditions) as a response to the applied external electric-like field to the graphene. In particular, such an effect is crucial for understanding the conductivity in the graphene, especially in the so-called non-linear regime [36,37,38,39,40,41]. The first experimental observation of non-linear current-voltage characteristics ($I-V$) of graphene devices and its interpretation in terms of the pair-creation has been recently reported in [42]. In the work [43], the quantum electronic and energy transport in the graphene at low carrier density and low temperatures when quantum interference effects are important were studied in the framework of strong field QED. A formulation of the Dirac model in the Fock space that includes an interaction of photons with fermions was considered in Ref. [44].

Because of a limited number of exactly solvable cases, approaches have been developed in parallel that make it possible, in the absence of appropriate exact solutions of the Dirac equation, to use certain approximate methods, including semiclassical and numerical, for nonperturbative calculations of quantum effects that are related to the vacuum instability. In this regard, it should be noted here the method of quantum kinetic equations (KE). In particular, this method is well adopted to using numerical calculations. In the framework of strong field QED, such an approach was also considered in Refs. [45,46,47,48,49,50] (equivalence of the KE method and other exact methods was demonstrated in some cases in References [51,52]) and then applied to problems of QED with strong external fields (see, for example, [53]). This approach was recently adapted to the model of single-layer graphene in [54,55,56].

In the present work, we propose a generalization of the KE method while taking into account an interaction of carriers with a photon reservoir in the graphene excited by a uniform, time-dependent electric field. It is assumed that the interaction of carriers with the external quasiclassical field is taken into account nonperturbatively, in contrast to their interaction with the quantized electromagnetic field. Thus, there appear collision integrals (CI), which take into account interaction with the quantized electromagnetic field in the single-photon approximation, corresponding to two channels: carrier redistribution by momenta as a result of a stimulated absorption or emission of photons and annihilation or creation of pairs. The kinetic description of these processes required the introduction of an appropriate KE for the photon subsystem. The Maxwell equations that describe the generation of an internal plasma field close the system of equations.

In what follows, we restrict ourselves by the study of quantum radiation processes within the framework of the KE approach [57,58]. We consider two different radiation mechanisms. One of them is a collective mechanism for the generation of plasma currents and waves (in the spatially uniform case, plasma oscillations [49,59]). The corresponding electromagnetic fields can leave the region of active action of the external field and, as a result, can be detected far from it outside the graphene plane. In the graphene the excitation region is limited by a simply connected surface and it is not difficult to find the electromagnetic field in the whole space from surface currents [60]. A theoretical and experimental study of this kind of collective radiation was completed in References [61,62] within the framework of an alternative dynamic approach [63,64] and, in [54], to the KE approach [65]. In addition to the above mentioned quasiclassical radiation, a quantum component of the radiation also exists due to elementary acts of interaction of the carriers with photons. In graphene, these mechanisms can be taken into account by analogy with the standard $D=3+1$ QED [57,58,66]. Outside of the KE approach, this radiation mechanism was studied in Ref. [67]. In the general case, these two mechanisms act self-consistently: the generation of the photon field leads to a back reaction problem of the second level, influencing plasma currents and oscillations of electromagnetic fields.

The approach that is considered in the present work is based on an adaptation of KE methods for a description of nonrelativistic plasma-like media in a nonequilibrium state (see, for example, [68,69,70]). The peculiarity of the system under consideration is that in the presence of external fields, massless quasiparticle states in the graphene are only observed indirectly and appear only through macroscopic characteristics, such as currents and the radiation, which are available for an observation during strongly nonequilibrium evolution at any time. A more detailed discussion of the KE methods that were used in the present work and, in particular, the role of the polarization phenomena is given below in Section 5.

The work is organized, as follows. In Section 2.1, we briefly describe (following in main References [54,55,56,65]) a nonperturbative KE approach for studying the evolution of carrier excitations in monolayer graphene under the action of an external quasiclassical time dependent electric field. The set of KE obtained here is an analog of the corresponding equations, which was used in strong field QED (see, e.g., [53]). In Section 2.2, the KE are generalized to include the interaction of carriers with the quantized electromagnetic field. They allow for us to consider a set of new problems: the radiation of the quantized field from the graphene, the photoproduction of carriers under the action of the quantized electromagnetic field, cascade processes, and so on. The obtained closed self-consistent system of KE appears as a result of an application of the truncation procedure of the Bogoliubov–Born–Green–Kirkwood–Yvon (BBGKY) chains of equations in the lower order of the perturbation theory with respect to the interaction of carriers with the photon reservoir. In the following Section 3, we set the problem of how to consistently treat the radiation of the quantized electromagnetic field in the physical system under consideration. We argue that for our purposes it is enough to analyze only photon KE. The corresponding CI are quite complex functionals containing distribution functions of the carrier and photon subsystems. Significant simplifications appear in the case of a small photon density. This case we call the low density approximation. We note that, in addition, everywhere is used the long wave approximation. Some of the results obtained in the framework of such approximations are analyzed in Section 4. They are then used in Section 5 in order to analyze the CI in the annihilation and the momentum redistribution channels. We also present parallel numerical estimates of some of the analytical results, and we are convinced of their qualitative coincidence. In Section 6, we present a summary of all the obtained results, point out open problems, and outlook possible perspectives.

2. Kinetic Equations Describing Quantum Excitations in Graphene Placed in an Electric Field

2.1. Kinetic Equations Describing Zero Order Processes

As was already mentioned in the Introduction, problems that are related to the vacuum instability can be studied in the framework of the KE approach. Here, we use such an approach to study the system of electronic excitations (charged carriers) and their interaction with electromagnetic fields in the graphene placed in an external electromagnetic field. Conditionally, the consideration is divided in two steps. On the first step, we consider the production of the charged carriers from the quantum vacuum in the graphene under the action of a spatially uniform time-dependent external electric field ${\mathbf{E}}_{\mathrm{ext}}\left(t\right)=\left(E_{\mathrm{ext}}^1\left(t\right),E_{\mathrm{ext}}^2\left(t\right)\right)$ situated in the graphene plane $(x^1,x^2)$, distracting from a possible interaction of the carriers with the quantized electromagnetic field and neglecting possible modification of the external field due to a back reaction. According to terminology that is accepted in studying the vacuum instability, in this first step we consider the KE approach for describing only zero order processes with respect to radiative corrections, see Refs. [13,14,15]. On the second step, we take into account the interaction of the carriers with the quantized electromagnetic field, as well as with the modified due to the back reaction external field and, then, on this base, we study the resulting electromagnetic field in the graphene. Partially, such a consideration takes into account the first order processes with respect to the radiative corrections. We stress that the first step consideration is based on the KE approach describing zero order processes in the graphene placed in an external electric field and it is nonperturbative with respect to the interaction with this field.

The external field is described by electromagnetic potentials $A_{\mathrm{ext}}^{\alpha }\left(t\right)$, $\alpha =1,2$. In the general case, the carriers in the graphene are subjected to an effective electric field,

$$E^{\alpha }\left(t\right)=E_{\mathrm{ext}}^{\alpha }\left(t\right)+E_{\mathrm{int}}^{\alpha }\left(t\right),\alpha =1,2,$$

(1)

given by potentials

$$A^{\alpha }\left(t\right)=A_{\mathrm{ext}}^{\alpha }\left(t\right)+A_{\mathrm{int}}^{\alpha }\left(t\right),$$

(2)

where the internal field $E_{\mathrm{int}}^{\alpha }\left(t\right)$ is induced by a back-reaction mechanism (as will be shown, $E_{\mathrm{int}}\left(t\right)$ can be neglected under some suppositions). We let the electric field be switched on at $t_{\mathrm{in}}$ and switched off at $t_{\mathrm{out}}$, so that the interaction between the Dirac field and the electric field vanishes at all time instants outside the interval $t\in \left(t_{\mathrm{in}},t_{\mathrm{out}}\right)$.

In order to describe the carrier quantum motion, we use the Dirac model of the graphene that describes the carries in a vicinity of one of the two Dirac points at boundaries of the Brillouin zone, see References [31,71,72] for a review. This model considers non interacting between themselves carriers placed in the external field. A wave function of a carrier is a two-component spinor $\psi (\mathbf{r},t)$. The latter satisfies the corresponding massless Dirac equation1;

$$\begin{array}{ccc}& & i\hslash \dot{\psi }(\mathbf{r},t)=h\left(t\right)\psi (\mathbf{r},t),h\left(t\right)={\mathrm{v}}_F\widehat{\mathbf{P}}\left(t\right)\sigma ,\mathbf{r}=\left(x^1,x^2\right),\hfill \\ & & \widehat{\mathbf{P}}\left(t\right)=\widehat{\mathbf{p}}+\frac{e}{c}\mathbf{A}\left(t\right),e>0,\widehat{\mathbf{p}}=-i\hslash \nabla .\hfill \end{array}$$

(3)

The wave function $\psi (\mathbf{p},t)$ in the momentum representation is defined by the decomposition

$$\psi (\mathbf{r},t)=\frac{1}{\sqrt{S}}\sum _{\mathbf{p}}\psi (\mathbf{p},t)e^{i\mathbf{pr}/\hslash },$$

(4)

where S is the area of the standard box regularization, and it satisfies the following equation:

$$i\hslash \dot{\psi }\left(\mathbf{p},t\right)=h_{\mathbf{p}}\left(t\right)\psi \left(\mathbf{p},t\right),h_{\mathbf{p}}\left(t\right)={\mathrm{v}}_F\mathbf{P}\left(t\right)\sigma ,$$

(5)

whereas, $\mathbf{P}$ is the kinetic momentum,

$$\mathbf{P}\left(t\right)=\mathbf{p}+\frac{e}{c}\mathbf{A}\left(t\right)=\left(P^1,P^2\right).$$

(6)

Let us perform an unitary transformation $\psi \left(\mathbf{p},t\right)=U\left(t\right)\phi \left(\mathbf{p},t\right)$, where the matrix U has the form [37]:

$$U\left(t\right)=\frac{1}{\sqrt{2}}\left(\begin{array}{cc}\exp (-i\varkappa /2)& \exp (-i\varkappa /2)\\ \exp (i\varkappa /2)& -\exp (i\varkappa /2)\end{array}\right).$$

(7)

Subsequently, we come to an auxiliary quasienergy eigenvalue problem for the transformed Hamiltonian ${\tilde{h}}_{\mathbf{p}}\left(t\right)=U^{†}\left(t\right)h_{\mathbf{p}}\left(t\right)U\left(t\right).$ We fix the parameter $\varkappa$ by the condition $tan\varkappa =P^2/P^1$, such that the corresponding quasienergy (the excitation quasienergy or a dispersion law) is $\epsilon (\mathbf{p},t)={\mathrm{v}}_F\sqrt{{\mathbf{P}}^2}$ and:

$${\tilde{h}}_{\mathbf{p}}\left(t\right)u_{\pm 1}=\pm \epsilon (\mathbf{p},t)u_{\pm 1},u_{+1}=\left(\begin{array}{c}1\\ 0\end{array}\right),u_{-1}=\left(\begin{array}{c}0\\ 1\end{array}\right).$$

(8)

The spinor $\phi \left(\mathbf{p},t\right)$ satisfies the equation

$$i\hslash \dot{\phi }(\mathbf{p},t)={\tilde{h}}_{\mathbf{p}}\left(t\right)\phi (\mathbf{p},t)+\frac{1}{2}\lambda \hslash {\sigma }_1\phi (\mathbf{p},t),$$

(9)

where the excitation function $\lambda (\mathbf{p},t)$ is determined from the equation $2i{U}^{†}\dot{U}=\lambda {\sigma }_1$ and it has the form:

$$\lambda (\mathbf{p},t)=\frac{e{v}_F^2[E_{\mathrm{ext}}^1\left(t\right)P^2-E_{\mathrm{ext}}^2\left(t\right)P^1]}{{\epsilon }^2(\mathbf{p},t)}.$$

(10)

Note that the Dirac Hamiltonians $h\left(t\right)$ do not commute at distinct time instants, $\left[{\tilde{h}}_{\mathbf{p}}\left(t\right),{\tilde{h}}_{\mathbf{p}}\left(t^{\prime }\right)\right]\ne 0$ if $t\ne t^{\prime }$. In the model under consideration, the dispersion law holds true in a vicinity of the Dirac point ${\mathbf{P}}^2=0$ at the boundaries of the Brillouin zone. This model corresponds to low-energy excitations of the carries. However, the approach under consideration allows for a generalization to a tight-binding model of the nearest neighbor interaction (see References [31,40,73]), as demonstrated in Ref. [54].

In what follows, we are going to consider the so-called adiabatic ansatz (alternatively called the quasiparticle representation), which is widely used in semiclassical approximations and numeric calculations, see, e.g., References [45,47,48,49,51,52,74].

We recall that there are two species of fermions in the model, corresponding to excitations about two distinct Dirac points in the Brillouin zone, i.e., each of species belongs to a distinct valley. The algebra of $\gamma$-matrices has two inequivalent representations in $(2+1)$-dimensions and a distinct (pseudospin) representation is associated with each Dirac point. Another doubling of fields is due to the (real) spin of the electron. Consequently, there are four species of fermions in the model. Thus, one should to take into account the degeneracy factor $N_f=4$ in order to find real mean values of physical quantities. We note that a transition to the quasiparticle representation in the model was used, for example, in Refs. [37,54,65,75]. Below, we follow the works by [54,65].

A quantum Dirac field $\mathsf{\Psi }(\mathbf{r},t)$ that is associated with the function $\psi (\mathbf{r},t)$ satisfies Equation (3) and the standard equal-time canonical anticommutation relations. It describes a fermion species of the model. The field operator $\mathsf{\Psi }(\mathbf{p},t)$ in the momentum representation is defined by the decomposition:

$$\mathsf{\Psi }(\mathbf{r},t)=\frac{1}{\sqrt{S}}\sum _{\mathbf{p}}\mathsf{\Psi }(\mathbf{p},t)e^{i\mathbf{pr}/\hslash }.$$

(11)

Subsequently, it is convenient to introduce the field operator $\mathsf{\Phi }(\mathbf{p},t)=U^{-1}\left(t\right)\mathsf{\Psi }\left(\mathbf{p},t\right)$, which satisfies Equation (9). According to the adiabatic ansatz, we define two kinds of creation and annihilation operators ($a^{†}(\mathbf{p},t),a(\mathbf{p},t)$ and $b^{†}(\mathbf{p},t),b(\mathbf{p},t)$) decomposing the operator $\mathsf{\Phi }\left(\mathbf{p},t\right)$ into solutions $u_{\pm 1}$ (see Equation (8)),

$$\mathsf{\Phi }(\mathbf{p},t)=a(\mathbf{p},t)u_{+1}+b^{†}(-\mathbf{p},t)u_{-1}.$$

(12)

Their nonzero anticommutation relations read:

$${\left[a(\mathbf{p},t),a^{†}({\mathbf{p}}^{\prime },t)\right]}_{+}={\left[b(\mathbf{p},t),b^{†}({\mathbf{p}}^{\prime },t)\right]}_{+}={\delta }_{\mathbf{p},{\mathbf{p}}^{\prime }}.$$

(13)

In each time instant t, one can formally introduce a Fock space that is equipped by instantaneous vacuum vectors $|0,t\rangle ,$

$$a(\mathbf{p},t)|0,t\rangle =b(\mathbf{p},t\left)\right|0,t\rangle =0,\forall \mathbf{p},$$

(14)

and a corresponding basis originated by the action of the creation operators $a^{†}(\mathbf{p},t)$ and $b^{†}(\mathbf{p},t)$ on the corresponding vacuum vectors. Equation (9) for the operator $\mathsf{\Phi }\left(\mathbf{p},t\right)$ implies the following equations for the creation and annihilation operators:

$$\begin{array}{ccc}& & i\hslash \dot{a}(\mathbf{p},t)=\epsilon (\mathbf{p},t)a(\mathbf{p},t)-\frac{1}{2}\hslash \lambda (\mathbf{p},t)b^{†}(-\mathbf{p},t),\hfill \\ & & i\hslash \dot{b}(-\mathbf{p},t)=\epsilon (\mathbf{p},t)b(-\mathbf{p},t)+\frac{1}{2}\hslash \lambda (\mathbf{p},t)a^{†}(-\mathbf{p},t).\hfill \end{array}$$

(15)

The QFT Hamiltonian and corresponding charge operator Q in the model have the form:

$$\begin{array}{ccc}& & H\left(t\right)=\int d\mathbf{r}{\mathsf{\Psi }}^{†}(\mathbf{r},t)h\left(t\right)\mathsf{\Psi }(\mathbf{r},t),\hfill \\ & & Q\left(t\right)=-\frac{e}{2}\int d\mathbf{r}\left[{\mathsf{\Psi }}^{†}(\mathbf{r},t),\mathsf{\Psi }(\mathbf{r},t)\right],\hfill \end{array}$$

(16)

where the integration is over the finite area S. They can be diagonalized at any time instant t while using decomposition (12),

$$\begin{array}{ccc}& & H\left(t\right)=\sum _{\mathbf{p}}\epsilon (\mathbf{p},t)\left[a^{†}(\mathbf{p},t)a(\mathbf{p},t)-b(-\mathbf{p},t)b^{†}(-\mathbf{p},t)\right],\hfill \\ & & Q=-e\sum _{\mathbf{p}}\left[a^{†}(\mathbf{p},t)a(\mathbf{p},t)-b^{†}(-\mathbf{p},t)b(-\mathbf{p},t)\right].\hfill \end{array}$$

(17)

The introduced above creation and annihilation operators become in- and out-operators of real quasiparticle at time instants $t_{\mathrm{in}}$ and $t_{\mathrm{out}}$, because the external electric field vanishes for $t\in \left(-\infty ,t_{\mathrm{in}}\right)\cup \left(t_{\mathrm{out}},+\infty \right).$ They act in the corresponding Fock spaces with initial and final vacua $|0,\mathrm{in}\rangle =|0,t_{\mathrm{in}}\rangle$ and $|0,\mathrm{out}\rangle =|0,t_{\mathrm{out}}\rangle$, respectively. One interprets $a^{†}(\mathbf{p},t_{\mathrm{in}})$ and $a(\mathbf{p},t_{\mathrm{in}})$ as creation and annihilation operators of initial electrons, $b^{†}(\mathbf{p},t_{\mathrm{in}})$ and $b(\mathbf{p},t_{\mathrm{in}})$ as the creation and annihilation operators of initial holes, whereas $a^{†}(\mathbf{p},t_{\mathrm{in}})a(\mathbf{p},t_{\mathrm{in}})$ and $b^{†}(\mathbf{p},t_{\mathrm{in}})b(\mathbf{p},t_{\mathrm{in}})$ are operators of initial electron and hole numbers, respectively. Operators $a^{†}(\mathbf{p},t_{\mathrm{out}})$ and $a(\mathbf{p},t_{\mathrm{out}})$ are interpreted as creation and annihilation operators of final electrons, $b^{†}(\mathbf{p},t_{\mathrm{out}})$ and $b(\mathbf{p},t_{\mathrm{out}})$ as creation and annihilation operators of final holes, whereas $a^{†}(\mathbf{p},t_{\mathrm{out}})a(\mathbf{p},t_{\mathrm{out}})$ and $b^{†}(\mathbf{p},t_{\mathrm{out}})b(\mathbf{p},t_{\mathrm{out}})$ are interpreted as operators of final electron and hole numbers, respectively.

Because the particles are massless in the model under consideration, and the spin degrees of freedom are hidden and manifest themselves only in the population of states, the charge is the only characteristic of the quasiparticles. It can be shown that Equation $\left[Q\left(t\right),H\left(t\right)\right]=0$ holds true at any time moment and, therefore, the electroneutrality of the system is preserved over time.

It is useful to introduce the following auxiliary distribution functions of quasiparticles (time-evolving adiabatic particle numbers)

$$\begin{array}{ccc}\hfill f^e(\mathbf{p},t)& =& \langle 0,\mathrm{in}\left|a^{+}(\mathbf{p},t)a(\mathbf{p},t)\right|0,\mathrm{in}\rangle ,\hfill \end{array}$$

(18)

$$\begin{array}{ccc}\hfill f^h(\mathbf{p},t)& =& \langle 0,\mathrm{in}\left|b^{+}(-\mathbf{p},t)b(-\mathbf{p},t)\right|0,\mathrm{in}\rangle .\hfill \end{array}$$

(19)

Averaging charge operator (16) over the in-vacuum, and taking into account the charge conservation low, we obtain

$$f^e(\mathbf{p},t)=f^h(\mathbf{p},t)=f(\mathbf{p},t).$$

(20)

Because $a(\mathbf{p},t_{\mathrm{in}})|0,\mathrm{in}\rangle =b(\mathbf{p},t_{\mathrm{in}})|0,\mathrm{in}\rangle =0,$ the initial value of the function $f(\mathbf{p},t)$ is zero, $f(\mathbf{p},t_{\mathrm{in}})=0$. When the electric field is turned off in an asymptotically distant future, the dispersion laws of quasiparticles goes to the mass surface,

$${\epsilon }_{\mathrm{out}}(\mathbf{p},t)={\mathrm{v}}_F\sqrt{{\left(\mathbf{p}+\frac{e}{c}{\mathbf{A}}_{\mathrm{out}}\right)}^2},{\mathbf{A}}_{\mathrm{out}}=\underset{t\to \infty }{lim}\mathbf{A}\left(t\right).$$

Subsequently, functions (18) and (19) describe the momentum distributions of real (observable) particles,

$$f(\mathbf{p},t_{\mathrm{out}})=\langle 0,\mathrm{in}|a^{+}(\mathbf{p},t_{\mathrm{out}})a(\mathbf{p},t_{\mathrm{out}})|0,\mathrm{in}\rangle =\langle 0,\mathrm{in}|b^{+}(\mathbf{p},t_{\mathrm{out}})b(\mathbf{p},t_{\mathrm{out}})|0,\mathrm{in}\rangle .$$

(21)

To obtain a closed set of KE, we differentiate $f(\mathbf{p},t)$ with respect to the time. Subsequently, using Equation (15), we obtain:

$$\dot{f}(\mathbf{p},t)=\frac{i\lambda \left(\mathbf{p},t\right)}{2}\left[f^{(+)}(\mathbf{p},t)-f^{(-)}(\mathbf{p},t)\right],$$

(22)

where the following anomalous expectation values

$$\begin{array}{ccc}& & f^{(+)}(\mathbf{p},t)=\langle 0,\mathrm{in}|a^{+}(\mathbf{p},t)b^{+}(-\mathbf{p},t)|0,\mathrm{in}\rangle ,\hfill \\ & & f^{(-)}(\mathbf{p},t)=\langle 0,\mathrm{in}\left|b(-\mathbf{p},t)a(\mathbf{p},t)\right|0,\mathrm{in}\rangle \hfill \end{array}$$

(23)

are introduced and $\lambda \left(\mathbf{p},t\right)$ is given by Equation (10).

The time derivatives of the functions $f^{(\pm )}(\mathbf{p},t)$ have the form:

$$\begin{array}{ccc}\hfill {\dot{f}}^{(+)}(\mathbf{p},t)& =& \frac{2i}{\hslash }\epsilon (\mathbf{p},t)f^{(+)}(\mathbf{p},t)-\frac{i\lambda (\mathbf{p},t)}{2}[1-2f(\mathbf{p},t)],\hfill \\ \hfill {\dot{f}}^{(-)}(\mathbf{p},t)& =& \frac{-2i}{\hslash }\epsilon (\mathbf{p},t)f^{(-)}(\mathbf{p},t)+\frac{i\lambda (\mathbf{p},t)}{2}[1-2f(\mathbf{p},t)].\hfill \end{array}$$

(24)

As a result of the integration of these equations over time and substitution into Equation (22), we obtain a KE of non-Markovian type in the following form:

$$\begin{array}{ccc}& & \dot{f}(\mathbf{p},t)=I(\mathbf{p},t),\hfill \\ & & I(\mathbf{p},t)=\frac{1}{2}\lambda (\mathbf{p},t){\int }_{t_0}^t\lambda (\mathbf{p},t^{\prime })[1-2f(\mathbf{p},t^{\prime })]\cos \theta (\mathbf{p};t,t^{\prime })d{t}^{\prime },\hfill \end{array}$$

(25)

where the phase $\theta$ in the source function $I(\mathbf{p},t)$ reads:

$$\theta (\mathbf{p};t,t^{\prime })=\frac{2}{\hslash }{\int }_{t^{\prime }}^{t}d{t}^{″}\epsilon (\mathbf{p},t^{″}).$$

(26)

The main task of the first stage is to find the distribution function $f\left(\mathbf{p},t\right)$ of created carriers. To this end, it is convenient to introduce new functions $u(\mathbf{p},t)$ and $\nu (\mathbf{p},t)$ and reduce integro-differential Equation (25) to an equivalent set of ordinary differential equations:

$$\begin{array}{ccc}& & \dot{f}(\mathbf{p},t)=\frac{1}{2}\lambda (\mathbf{p},t)u(\mathbf{p},t),\hfill \\ & & \dot{u}(\mathbf{p},t)=\lambda (\mathbf{p},t)\left[1-2f(\mathbf{p},t)\right]-\frac{2\epsilon (\mathbf{p},t)}{\hslash }\nu (\mathbf{p},t),\dot{\nu }(\mathbf{p},t)=\frac{2\epsilon (\mathbf{p},t)}{\hslash }u(\mathbf{p},t).\hfill \end{array}$$

(27)

The functions $u(\mathbf{p},t)$ and $\nu (\mathbf{p},t)$ describe the polarization effects and they are expressed in terms of anomalous means (23) as:

$$u(\mathbf{p},t)=f^{(+)}(\mathbf{p},t)+f^{(-)}(\mathbf{p},t),\nu (\mathbf{p},t)=i\left[f^{(+)}(\mathbf{p},t)-f^{(-)}(\mathbf{p},t)\right].$$

(28)

The total Hamiltonian of the fermion subsystem $H_{\mathrm{tot}}\left(t\right)=H\left(t\right)+H_{\mathrm{pol}}\left(t\right)$ contains two parts: the Hamiltonian $H\left(t\right)$ of the quasiparticle excitations (17) and a polarization Hamiltonian $H_{\mathrm{pol}}\left(t\right)$, which corresponds to the second term in the RHS of Equation (9),

$$H_{\mathrm{pol}}\left(t\right)=-i\frac{\hslash }{2}\sum _{\mathbf{p}}\lambda (\mathbf{p},t)\left[a^{†}(\mathbf{p},t)b^{†}(-\mathbf{p},t)-b(-\mathbf{p},t)a(\mathbf{p},t)\right].$$

(29)

The corresponding vacuum polarization energy density can be represented as

$$E_{\mathrm{pol}}\left(t\right)=-\frac{\hslash }{2S}\sum _{\mathbf{p}}\lambda (\mathbf{p},t)v(\mathbf{p},t),$$

(30)

where $v(\mathbf{p},t)$ is given by Equation (28).

The vacuum mean value of the current density

$$\mathbf{j}(\mathbf{r},t)=\langle 0,\mathrm{in}\left|e{\mathrm{v}}_F{\mathsf{\Psi }}^{†}(\mathbf{r},t)\sigma \mathsf{\Psi }(\mathbf{r},t)\right|0,\mathrm{in}\rangle$$

can be written in the following form:

$$\mathbf{j}\left(t\right)={\mathbf{j}}_{\mathrm{cond}}\left(t\right)+{\mathbf{j}}_{\mathrm{pol}}\left(t\right),$$

(31)

see [54,65]. Here, the conductivity ${\mathbf{j}}_{\mathrm{cond}}\left(t\right)$ and polarization current ${\mathbf{j}}_{\mathrm{pol}}\left(t\right)$ densities are:

$$\begin{array}{ccc}\hfill j_{\mathrm{cond}}^{\alpha }\left(t\right)& =& \frac{2e}{S}\sum _{\mathbf{p}}{v}^{\alpha }(\mathbf{p},t)f(\mathbf{p},t),\hfill \end{array}$$

(32)

$$\begin{array}{ccc}\hfill j_{\mathrm{pol}}^{\alpha }\left(t\right)& =& -\frac{e}{S}\sum _{\mathbf{p}}{\tilde{v}}^{\alpha }(\mathbf{p},t)u(\mathbf{p},t),\hfill \end{array}$$

(33)

where

$$v^{\alpha }(\mathbf{p},t)=\frac{\partial \epsilon \left(\mathbf{p},t\right)}{\partial P^{\alpha }}=\frac{{\mathrm{v}}_F^2{P}^{\alpha }}{{\epsilon }^2\left(\mathbf{p},t\right)}$$

(34)

is a group velocity and ${\tilde{v}}^{\alpha }(\mathbf{p},t)$ is the vector of the so-called conjugate velocity that is determined by components of vector (34), as follows: ${\tilde{v}}^{\alpha }(\mathbf{p},t)=\left(v^2,-v^1\right)$. Note that the total current of all fermion species is $N_f\mathbf{j}\left(t\right)$. Section 5 discusses the role of current densities (32) and (33) in the electromagnetic emission in the graphene.

The plasma classical electric field ${\mathbf{E}}_{\mathrm{int}}\left(t\right)$ is generated by the internal current $N_f\mathbf{j}\left(t\right)$ and it satisfies the Maxwell equation,

$${\dot{\mathbf{E}}}_{\mathrm{int}}\left(t\right)=-N_f\mathbf{j}\left(t\right).$$

(35)

This field contributes to the effective quasiclassical electric field (1), which, in turn, has an effect on the dynamics of carriers according to Equation (27).

The above formulation of the quantum kinetic theory describing excitations in the graphene [54,65] is constructed by analogy with the quantum kinetic theory describing the vacuum production of $e^{-}{e}^{+}$ plasma (see, e.g., [46,76]).

In the thermodynamical limit $S\to \infty$, replacing the sum over the momenta in Equations (32) and (33) by an integral the conductivity and polarization current densities take the form

$$\begin{array}{ccc}& & {\mathbf{j}}_{\mathrm{cond}}\left(t\right)=\frac{2e}{{(2\pi \hslash )}^2}\int \mathbf{v}(\mathbf{p},t)f(\mathbf{p},t)d\mathbf{p},\hfill \end{array}$$

(36)

$$\begin{array}{ccc}& & {\mathbf{j}}_{\mathrm{pol}}\left(t\right)=-\frac{e}{{(2\pi \hslash )}^2}\int \epsilon (\mathbf{p},t)\mathbf{l}(\mathbf{p},t)u(\mathbf{p},t)d\mathbf{p},\hfill \end{array}$$

(37)

where $\mathbf{v}(\mathbf{p},t)$ given by Equation (34) is a propagation velocity of quasiparticle excitations, $\mathbf{l}(\mathbf{p},t)$ is a polarization function,

$$\begin{array}{ccc}& & e\mathbf{l}(\mathbf{p},t)=(\partial \lambda (\mathbf{p},t)/\partial E_{\mathrm{ext}}^1\left(t\right),\partial \lambda (\mathbf{p},t)/\partial E_{\mathrm{ext}}^2\left(t\right)),\hfill \\ & & l_1(\mathbf{p},t)=v_F^2{P}_2/{\epsilon }^2(\mathbf{p},t),l_2(\mathbf{p},t)=-v_F^2{P}_1/{\epsilon }^2(\mathbf{p},t).\hfill \end{array}$$

(38)

The function $\lambda (\mathbf{p},t)$ (10) can be expressed via $\mathbf{l}(\mathbf{p},t)$, as follows:

$$\lambda (\mathbf{p},t)=e{\mathbf{E}}_{\mathrm{ext}}\left(t\right)\mathbf{l}(\mathbf{p},t).$$

(39)

Thus, we see that the polarization energy (30) and polarization current (37) are expressed in terms of polarization functions (28). We note that the separation of the total current in the sum of the conduction current and the polarization current is, in a certain sense, conditional, but useful in approximate calculations.

Among the macroscopic averages, it is necessary to also refer the number density of pairs of $N_f$ fermion species

$$n\left(t\right)=\frac{N_f}{{(2\pi \hslash )}^2}\int f(\mathbf{p},t)d\mathbf{p},$$

(40)

Following the works [46,50], we can discover that, for a finite t the asymptotic behavior of solutions of system (27) for $\epsilon (p,t)\to \infty$ is described by the following leading terms

$$f(\mathbf{p},t)\approx \frac{1}{16}{\left(\frac{\lambda (\mathbf{p},t)}{\epsilon (\mathbf{p},t)}\right)}^2,u(\mathbf{p},t)\approx \frac{1}{4\epsilon (\mathbf{p},t)}\frac{d}{dt}\frac{\lambda (\mathbf{p},t)}{\epsilon (\mathbf{p},t)},\nu (\mathbf{p},t)\approx \frac{\lambda (\mathbf{p},t)}{2\epsilon (\mathbf{p},t)}.$$

(41)

Therefore, integrals (36), (37), and (40) converge. Note that, for $t\to \infty$ the term (37) containing $u(p,t)$ vanishes because it is represented via an integral over momenta of a rapidly oscillating function. Assuming the external field being switched off for $t\to \infty$ the term $f(p,t)$ exponentially disappears with the growth of $\epsilon (p,t)$.

2.2. Inclusion an Interaction with Quantized Electromagnetic Field

Here, we assume that excitations in the graphene interact with both effective classical electromagnetic field (2) and the quantized electromagnetic field situated in the graphene plane. Below, we generalize the KE for such an extended system. The total QFT Hamiltonian $H_{\mathrm{tot}}\left(t\right)$ that corresponds to the system consists of a part that is originated from the Dirac model of excitations in the graphene interacting with both an external classical electromagnetic field and with the quantized electromagnetic field and of a Hamiltonian of the quantized electromagnetic field. Note that we explicitly consider an interaction with only one of a fermion species of the model. Thus, to find real mean values of physical quantities, one should to take the degeneracy factor into account $N_f=4$.

In the quasiparticle representation, as introduced in Section 2.1, such a Hamiltonian reads:

$$H_{\mathrm{tot}}\left(t\right)=H\left(t\right)+H_{\mathrm{pol}}\left(t\right)+H_{\mathrm{int}}\left(t\right)+H_{\mathrm{A}}\left(t\right),$$

(42)

where Hamiltonians $H\left(t\right)$, as given by Equation (17), and $H_{\mathrm{pol}}\left(t\right)$ given by Equation (28) describe excitations in the graphene and their interaction with existing in the graphene quasiclassical electromagnetic field ${\mathbf{E}}_{\mathrm{ext}}\left(t\right)$ whereas $H_{\mathrm{A}}\left(t\right)$ is a Hamiltonian of the free quantized electromagnetic field in the Coulomb gauge.

The electromagnetic field is not confined to the graphene surface, $z=0$, but it rather propagates in the ambient $3+1$ dimensional space-time, where z is the coordinate of axis normal to the graphene plane. We allow for the graphene sheet to have a global momentum $p_z$ along the z axis, in order to account for the possibility of a momentum transfer in this direction to some external system. However, only the projection of the electromagnetic field operator on the graphene plane $\mathbf{A}\left(x\right)=\left(A^1\left(x\right),A^2\left(x\right)\right)$ interacts with electrons and holes. Thus, we can exclude the noninteracting component, $A^3\left(x\right)$ of the electromagnetic field operator from the consideration and which corresponds to the case $z=0$. If the wave vector of a photon field, $\left(k^1,k^2,k_z\right)$, has component $k_z\ne 0$ then the photon immediately leaves the interaction area. Such a free photon field can be represented with standard annihilation and creation operators and formulate, in such terms, a theory of emission in the first-order approximation with respect to electron-photon interaction [44]. However, in the present model, we are interested in a different emission mechanism, which allows for one to consider the photon field propagating along the graphene plane and having a sufficiently small momenta $k_z\to 0$. We chose the range of $k_z$ so as to emitted photon did not leave a graphene film of thickness d.Such a condition determines a limiting angle between the emission direction and graphene plane,

$$\left|k_z\right|/k<\phi ,\phi \sim d/\sqrt{S}\ll 1,$$

(43)

where $k=\left|\mathbf{k}\right|$ and $\mathbf{k}=\left(k^1,k^2\right)$ is a two-component wave vector.

The operator of such a field can be written in the following form:

$$\begin{array}{ccc}& & {\widehat{A}}^{\alpha }(\mathbf{r},t)={\widehat{A}}^{(+)\alpha }(\mathbf{r},t)+{\widehat{A}}^{(-)\alpha }(\mathbf{r},t),\alpha =1,2,\hfill \\ & & {\widehat{A}}^{(\pm )\alpha }(\mathbf{r},t)=\sqrt{\frac{\hslash c}{V}}\sum _{\mathbf{k}}\sum _{k_z<\phi k}\frac{1}{\sqrt{2k}}{\widehat{A}}^{(\pm )\alpha }(\pm \mathbf{k},k_z,t)e^{-i\mathbf{kr}},\hfill \\ & & k=\left|\mathbf{k}\right|,{\left[{\widehat{A}}^{(\pm )\alpha }(\mathbf{k},k_z,t)\right]}^{†}=A^{(\mp )\alpha }(-\mathbf{k},-k_z,t),\hfill \end{array}$$

(44)

where $V=SL$ is the volume of the regularization box (L is the length of the edge of the box normal to the graphene plane). This field is polarized and its polarization vector is transversal to the wavevector $\mathbf{k}$. The operators ${\widehat{A}}^{(\pm )\alpha }(\mathbf{k},k_z,t)$ satisfy the standard equal time commutation relations:

$$\begin{array}{ccc}\hfill \left[A^{(-)\alpha }(\mathbf{k},k_z,t),A^{(+)\beta }({\mathbf{k}}^{\prime },k_z^{\prime },t)\right]& =& {\delta }_{\alpha \beta }{\delta }_{{\mathbf{kk}}^{\prime }}{\delta }_{k_z,k_z^{\prime }},\hfill \\ \hfill \left[A^{(\pm )\alpha }(\mathbf{k},k_z,t),A^{(\pm )\beta }({\mathbf{k}}^{\prime },k_z^{\prime },t)\right]& =& 0.\hfill \end{array}$$

(45)

Then, the Hamiltonian $H_{\mathrm{A}}\left(t\right)$ of the free quantized electromagnetic field can be written as

$$H_{\mathrm{A}}\left(t\right)=c\hslash \sum _{\mathbf{k}}\sum _{k_z<\phi k}k{\widehat{A}}^{(+)\alpha }(\mathbf{k},k_z,t){\widehat{A}}^{(-)\alpha }(\mathbf{k},k_z,t),$$

(46)

The effective interaction of the carriers in graphene with the quantized electromagnetic field can be defined as:

$$H_{\mathrm{int}}\left(t\right)=-e\frac{v_F}{c}\int d\mathbf{r}{\mathsf{\Psi }}^{†}(\mathbf{r},t)\sigma \mathsf{\Psi }(\mathbf{r},t)\widehat{\mathbf{A}}(\mathbf{r},t).$$

(47)

Substituting decompositions (11), (12), and (44) in Hamiltonian (47), we obtain:

$$\begin{array}{ccc}& & H_{\mathrm{int}}\left(t\right)=-e{v}_F\sum _{\mathbf{p},\mathbf{k}}\sum _{k_z<\phi k}\sqrt{\frac{\hslash }{2cVk}}\left\{{\mathsf{\Gamma }}_{uu}^{\alpha }(\mathbf{p},\mathbf{p}-\hslash \mathbf{k};t)a^{†}(\mathbf{p},t)a(\mathbf{p}-\hslash \mathbf{k},t)\right.\hfill \\ & & +{\mathsf{\Gamma }}_{uv}^{\alpha }(\mathbf{p},\mathbf{p}-\hslash \mathbf{k};t)a^{†}(\mathbf{p},t)b^{†}(-\mathbf{p}+\hslash \mathbf{k},t)+{\mathsf{\Gamma }}_{vu}^{\alpha }(\mathbf{p},\mathbf{p}-\hslash \mathbf{k};t)b(-\mathbf{p},t)a(\mathbf{p}-\hslash \mathbf{k},t)\hfill \\ & & +\left.{\mathsf{\Gamma }}_{vv}^{\alpha }(\mathbf{p},\mathbf{p}-\hslash \mathbf{k};t)b(-\mathbf{p},t)b^{†}(-\mathbf{p}+\hslash \mathbf{k},t)\right\}{\widehat{A}}^{\alpha }(\mathbf{k},k_z,t).\hfill \end{array}$$

(48)

Here, ${\mathsf{\Gamma }}_{\xi \eta }^{\alpha }(\mathbf{p},{\mathbf{p}}^{\prime },t)$ are vertex matrix functions,

$${\mathsf{\Gamma }}_{\xi \eta }^{\alpha }(\mathbf{p},{\mathbf{p}}^{\prime },t)={\xi }^{†}{U}^{†}(\mathbf{p},t){\sigma }^{\alpha }U({\mathbf{p}}^{\prime },t)\eta ,{\mathsf{\Gamma }}_{\xi \eta }^{\alpha *}(\mathbf{p},{\mathbf{p}}^{\prime };t)={\mathsf{\Gamma }}_{\eta \xi }^{\alpha }({\mathbf{p}}^{\prime },\mathbf{p};t),$$

(49)

spinors $\xi$ and $\eta$ are given by Equation (8), and the evolution matrix $U(\mathbf{p},t)$, which describes the influence of the external field on the interaction of the introduced quasiparticles with photons is given by Equation (7).

Hamiltonian (42), being written in the quasiparticle representation, determines the time evolution of creation and annihilation operators of all the particles,

$$\begin{array}{ccc}& & \dot{a}(\mathbf{p},t)=-\frac{i}{\hslash }\epsilon (\mathbf{p},t)a(\mathbf{p},t)+\frac{i}{2}\lambda (\mathbf{p},t)b^{†}(-\mathbf{p},t)+e{v}_F\sum _{\mathbf{k}}\sum _{k_z<\phi k}\frac{i}{\sqrt{2\hslash cVk}}\hfill \\ & & \times \left\{{\mathsf{\Gamma }}_{uu}^{\alpha }(\mathbf{p},{\mathbf{p}}^{\prime };t)a(\mathbf{p}-\hslash K,t)+{\mathsf{\Gamma }}_{uv}^{\alpha }(\mathbf{p},{\mathbf{p}}^{\prime };t)b^{†}(-{\mathbf{p}}^{\prime },t)\right\}{\widehat{A}}^{\alpha }(\mathbf{k},k_z,t),\hfill \\ & & \dot{b}(-\mathbf{p},t)=-\frac{i}{\hslash }b(-\mathbf{p},t)-\frac{i}{2}\lambda (\mathbf{p},t)a^{†}(\mathbf{p},t)-e{v}_F\sum _{\mathbf{k}}\sum _{k_z<\phi k}\frac{i}{\sqrt{2\hslash cVk}}\hfill \\ & & \times \left\{{\mathsf{\Gamma }}_{uv}^{\alpha }(\mathbf{p}+\hslash K,\mathbf{p};t)a^{†}(\mathbf{p}+\hslash \mathbf{k},t)+{\mathsf{\Gamma }}_{vv}^{\alpha }(\mathbf{p}+\hslash K,\mathbf{p};t)b(-{\mathbf{p}}^{\prime },t)\right\}{\widehat{A}}^{\alpha }(\mathbf{k},k_z,t),\hfill \\ & & {\dot{\widehat{A}}}^{(\pm )\alpha }(\pm \mathbf{k},k_z,t)=\pm iK{\mathcal{A}}^{(\pm )\oslash }(\pm \mathbf{k},k_z,t)\mp e{v}_F\sum _{\mathbf{p}}\frac{i}{\sqrt{2\hslash cVk}}\hfill \\ & & \left\{{\mathsf{\Gamma }}_{uu}^{\alpha }(\mathbf{p},{\mathbf{p}}^{\prime };t)a^{†}(\mathbf{p},t)a({\mathbf{p}}^{\prime },t)+{\mathsf{\Gamma }}_{uv}^{\alpha }(\mathbf{p},{\mathbf{p}}^{\prime };t)a^{†}(\mathbf{p},t)b^{†}(-{\mathbf{p}}^{\prime },t)\right.\hfill \\ & & +\left.{\mathsf{\Gamma }}_{vu}^{\alpha }(\mathbf{p},{\mathbf{p}}^{\prime };t)b(-\mathbf{p},t)a({\mathbf{p}}^{\prime },t)+{\mathsf{\Gamma }}_{vv}^{\alpha }(\mathbf{p},{\mathbf{p}}^{\prime };t)b(-\mathbf{p},t)b^{†}(-{\mathbf{p}}^{\prime },t)\right\},\hfill \end{array}$$

(50)

where ${\mathbf{p}}^{\prime }=\mathbf{p}-\hslash \mathbf{k}$.

Now, we introduce the correlation function $f(\mathbf{p},\mathbf{p}{}^{\prime };t)$ of the electron-hole subsystem and the one $F_{\alpha \beta }(\mathbf{k},\mathbf{k}{}^{\prime };t)$ of the photon subsystem,

$$\begin{array}{ccc}& & f(\mathbf{p},\mathbf{p}{}^{\prime };t)=\langle 0,\mathrm{in}|a^{†}(\mathbf{p},t)a(\mathbf{p}{}^{\prime },t)|0,\mathrm{in}\rangle ,\hfill \\ & & F_{\alpha \beta }(\mathbf{k},\mathbf{k}{}^{\prime };t)=\sum _{k_z<\phi k}\langle 0,\mathrm{in}|{\widehat{A}}^{(+)\alpha }(\mathbf{k},k_z,t){\widehat{A}}^{(-)\beta }(\mathbf{k}{}^{\prime },k_z^{\prime },t)|0,\mathrm{in}\rangle .\hfill \end{array}$$

(51)

In the case of spatially-homogeneous systems, the correlation functions, being written in the coordinate representation, only depend on the coordinate difference. Their Fourier transforms depend of the corresponding $\delta$ functions. Note that, under the condition $\phi \to 0$ the momentum distribution of photons does not depend on $k_z$ and the correlation function $F_{\alpha \beta }(k,k{}^{\prime };t)$ is simply proportional to the number of state with $k_z<\phi k$, $n\left(k\right)=\frac{kd}{2\pi }$. Thus, the correlation functions (51) have the form:

$$\begin{array}{ccc}& & f(\mathbf{p},{\mathbf{p}}^{\prime };t)=f(\mathbf{p},t)\frac{S}{{\left(2\pi \right)}^2}{\delta }_{\mathbf{p},{\mathbf{p}}^{\prime }},\hfill \end{array}$$

(52)

$$\begin{array}{ccc}& & F_{\alpha \beta }(\mathbf{k},{\mathbf{k}}^{\prime };t)=F(\mathbf{k},t)n\left(k\right)\frac{S}{{\left(2\pi \right)}^2}{\delta }_{\alpha \beta }{\delta }_{\mathbf{k},{\mathbf{k}}^{\prime }},\hfill \end{array}$$

(53)

where

$$F(\mathbf{k},t)=\frac{1}{2}\sum _{\alpha }{F}_{\alpha \alpha }(\mathbf{k},{\mathbf{k}}^{\prime };t){|}_{{\mathbf{k}}^{\prime }=\mathbf{k}}.$$

(54)

In Equation (52), $f(\mathbf{p},t)$ is a superficial density of electron (hole) numbers with a momentum $\mathbf{p}$ at the time instant $t,$ whereas $F(\mathbf{k},t)$ in Equation (53) is a superficial density of photons with the wave vector $\mathbf{k}$ at the time instant t per unit state with a given $k_z$. In what follows, we use the notation ${\widehat{A}}^{\alpha }(k,t)={\widehat{A}}^{\alpha }(k,0,t)$ for the photon operator with a given $k_z\to 0$.

Equation (50) imply the following equations for the introduced distribution functions:

$$\begin{array}{ccc}& & \dot{f}(\mathbf{p},t)=\frac{1}{2}\lambda (\mathbf{p},t)u(\mathbf{p},t)-ie{v}_F\sum _{\mathbf{k}}\frac{{\left(2\pi \right)}^{2}n\left(k\right)}{\sqrt{2\hslash cVk}}\hfill \\ & & \times \left\{{\mathsf{\Gamma }}_{uu}^{\alpha *}(\mathbf{p},\mathbf{p}-\hslash \mathbf{k};t)\langle 0,\mathrm{in}|a^{†}(\mathbf{p}-\hslash \mathbf{k},t)a(\mathbf{p},t){\widehat{A}}^{\alpha }(\mathbf{k},t)|0,\mathrm{in}\rangle \right.\hfill \\ & & +{\mathsf{\Gamma }}_{uv}^{\alpha *}(\mathbf{p},\mathbf{p}-\hslash \mathbf{k};t)\langle 0,\mathrm{in}|b(-\mathbf{p}+\hslash \mathbf{k},t)a(\mathbf{p},t){\widehat{A}}^{\alpha }(\mathbf{k},t)|0,\mathrm{in}\rangle \hfill \\ & & -{\mathsf{\Gamma }}_{uu}^{\alpha }(\mathbf{p},\mathbf{p}-\hslash \mathbf{k};t)\langle 0,\mathrm{in}|a^{†}(\mathbf{p},t)a(\mathbf{p}-\hslash \mathbf{k},t){\widehat{A}}^{\alpha }(\mathbf{k},t)|0,\mathrm{in}\rangle \hfill \\ & & -\left.{\mathsf{\Gamma }}_{uv}^{\alpha }(\mathbf{p},\mathbf{p}-\hslash \mathbf{k};t)\langle 0,\mathrm{in}|a^{†}(\mathbf{p},t)b^{†}(-\mathbf{p}-\hslash \mathbf{k},t){\widehat{A}}^{\alpha }(\mathbf{k},t)|0,\mathrm{in}\rangle \right\},\hfill \end{array}$$

(55)

and

$$\begin{array}{ccc}& & \dot{F}(\mathbf{k},t)=-ie{v}_F\sum _{\mathbf{p}}\frac{2{\pi }^2}{\sqrt{\hslash cVk}}\left\{{\mathsf{\Gamma }}_{uu}^{\alpha }(\mathbf{p},\mathbf{p}+\hslash \mathbf{k};t)\langle 0,\mathrm{in}|a^{†}(\mathbf{p},t)a(\mathbf{p}+\hslash \mathbf{k},t){\widehat{A}}^{(-)\alpha }(\mathbf{k},t)|0,\mathrm{in}\rangle \right.\hfill \\ & & +{\mathsf{\Gamma }}_{uv}^{\alpha }(\mathbf{p},\mathbf{p}+\hslash \mathbf{k};t)\langle 0,\mathrm{in}|a^{†}(\mathbf{p},t)b^{†}(-\mathbf{p}-\hslash \mathbf{k},t){\widehat{A}}^{(-)\alpha }(\mathbf{k},t)|0,\mathrm{in}\rangle \hfill \\ & & +{\mathsf{\Gamma }}_{vu}^{\alpha }(\mathbf{p},\mathbf{p}+\hslash \mathbf{k};t)\langle 0,\mathrm{in}|b(-\mathbf{p},t)a(\mathbf{p}+\hslash \mathbf{k},t){\widehat{A}}^{(-)\alpha }(\mathbf{k},t)|0,\mathrm{in}\rangle \hfill \\ & & +{\mathsf{\Gamma }}_{vv}^{\alpha }(\mathbf{p},\mathbf{p}+\hslash \mathbf{k};t)\langle 0,\mathrm{in}|b(-\mathbf{p},t)b^{†}(-\mathbf{p}-\hslash \mathbf{k},t){\widehat{A}}^{(-)\alpha }(\mathbf{k},t)|0,\mathrm{in}\rangle \hfill \\ & & -{\mathsf{\Gamma }}_{uu}^{\alpha }(\mathbf{p},\mathbf{p}-\hslash \mathbf{k};t)\langle 0,\mathrm{in}|{\widehat{A}}^{(+)\alpha }(\mathbf{k},t)a^{†}(\mathbf{p},t)a(\mathbf{p}-\hslash \mathbf{k},t)|0,\mathrm{in}\rangle \hfill \\ & & -{\mathsf{\Gamma }}_{uv}^{\alpha }(\mathbf{p},\mathbf{p}-\hslash \mathbf{k};t)\langle 0,\mathrm{in}|{\widehat{A}}^{(+)\alpha }(\mathbf{k},t)a^{†}(\mathbf{p},t)b^{†}(-\mathbf{p}+\hslash \mathbf{k},t)|0,\mathrm{in}\rangle \hfill \\ & & -{\mathsf{\Gamma }}_{vu}^{\alpha }(\mathbf{p},\mathbf{p}-\hslash \mathbf{k};t)\langle 0,\mathrm{in}|{\widehat{A}}^{(+)\alpha }(\mathbf{k},t)b(-\mathbf{p},t)a(p-\hslash \mathbf{k},t)|0,\mathrm{in}\rangle \hfill \\ & & -\left.{\mathsf{\Gamma }}_{vv}^{\alpha }(\mathbf{p},\mathbf{p}-\hslash \mathbf{k};t)\langle 0,\mathrm{in}|{\widehat{A}}^{(+)\alpha }(\mathbf{k},t)b(-\mathbf{p},t)b^{†}(-\mathbf{p}+\hslash \mathbf{k},t)|0,\mathrm{in}\rangle \right\}.\hfill \end{array}$$

(56)

Equation (55) reduces to the first equation of set (27) when the electron-hole system does not interact with photons

The system of Equations (55) and (56) is not closed, their RHS contain higher order correlation functions. The means $\langle 0,\mathrm{in}|a^{†}a{\widehat{A}}^{(\pm )}|0,\mathrm{in}\rangle$ and $\langle 0,\mathrm{in}|b^{†}b{\widehat{A}}^{(\pm )}|0,\mathrm{in}\rangle$ correspond to processes of induced irradiation and absorption, whereas the means $\langle 0,\mathrm{in}|ab{\widehat{A}}^{(+)}|0,\mathrm{in}\rangle$ and $\langle 0,\mathrm{in}|a^{†}{b}^{†}{\widehat{A}}^{(-)}|0,\mathrm{in}\rangle$ correspond to processes with pair annihilation and creation. We do not consider the phantom process with the correlators $\langle 0,\mathrm{in}|a^{†}{b}^{†}{\widehat{A}}^{(+)}|0,\mathrm{in}\rangle$ and $\langle 0,\mathrm{in}|ab{\widehat{A}}^{(-)}|0,\mathrm{in}\rangle$.

Equations (55) and (56) constitute a part of the Bogoliubov–Born–Green–Kirkwood–Yvon (BBGKY) chain of equations in the electron-hole and photon sectors. In order to close the set (55) and (56), it is necessary to obtain equations of the second level for all of the above mentioned correlation functions. These equations will already contain two-particle correlation functions, which, assuming a weak interaction between the subsystems, can again be represented through the single-particle correlator.

Below, we give an example of a truncation procedure:

$$\begin{array}{ccc}& & \langle 0,\mathrm{in}|a^{†}({\mathbf{p}}_1,t)a({\mathbf{p}}_2,t){\widehat{A}}_{\beta }(\mathbf{k}{}^{\prime },t){\widehat{A}}_{\alpha }^{(\pm )}(\mathbf{k},t)|0,\mathrm{in}\rangle =\frac{S^2}{{\left(2\pi \right)}^4}f({\mathbf{p}}_1,t){\delta }_{{\mathbf{p}}_1,{\mathbf{p}}_2}{\delta }_{\alpha \beta }\left\{\begin{array}{c}\left[F(\mathbf{k},t)+1\right]{\delta }_{\mathbf{k},-{\mathbf{k}}^{\prime }}\\ F(\mathbf{k},t){\delta }_{\mathbf{k},{\mathbf{k}}^{\prime }}\end{array}\right.,\hfill \\ & & \langle 0,\mathrm{in}\left|a^{†}({\mathbf{p}}_1,t)a({\mathbf{p}}_2,t)a^{†}(\mathbf{p}{}^{\prime },t)a(\mathbf{p}{}^{″},t)\right|0,\mathrm{in}\rangle \hfill \\ & =& \frac{S^2}{{\left(2\pi \right)}^4}f({\mathbf{p}}_1,t)\left[f({\mathbf{p}}^{\prime },t){\delta }_{{\mathbf{p}}_1,{\mathbf{p}}_2}{\delta }_{{\mathbf{p}}^{\prime },{\mathbf{p}}^{″}}+\left[1-f({\mathbf{p}}^{\prime },t)\right]{\delta }_{{\mathbf{p}}_1,{\mathbf{p}}^{″}}{\delta }_{{\mathbf{p}}_2,{\mathbf{p}}^{\prime }}\right],\hfill \\ & & \langle 0,\mathrm{in}|a^{†}({\mathbf{p}}_1,t)a({\mathbf{p}}_2,t)b(-{\mathbf{p}}^{\prime },t)b^{†}(-{\mathbf{p}}^{″},t)|0,\mathrm{in}\rangle =\frac{S^2}{{\left(2\pi \right)}^4}f({\mathbf{p}}_1,t)\left[1-f({\mathbf{p}}^{\prime },t)\right]{\delta }_{{\mathbf{p}}_1,{\mathbf{p}}_2}{\delta }_{{\mathbf{p}}^{\prime },{\mathbf{p}}^{″}}.\hfill \end{array}$$

(57)

It corresponds to the random-phase-approximation (RPA) [69,70]. We note that, in the truncation procedure, one neglects the polarization effects in the resulting CI and in the evolution equations for the polarization functions $u(\mathbf{p},t)$ and $v(p,t)$ (27). Here, the assumption of spatial homogeneity and its consequence (51) were taken into account.

In such a way, a closed set of equations for the distribution functions $f(\mathbf{p},t)$ and $F(\mathbf{k},t)$ can be obtained. In the thermodynamic limit $V\to \infty$, it reads:

$$\begin{array}{ccc}\hfill \dot{f}(\mathbf{p},t)& =& I(\mathbf{p},t)+C_{\gamma }(\mathbf{p},t)+C_{\mathrm{eh}}(\mathbf{p},t),\hfill \end{array}$$

(58)

$$\begin{array}{ccc}\hfill \dot{F}(\mathbf{k},t)& =& S_{\gamma }(\mathbf{k},t)+S_{\mathrm{eh}}(\mathbf{k},t).\hfill \end{array}$$

(59)

Collision integrals $C_{\gamma }(\mathbf{p},t)$ and $C_{\mathrm{eh}}(\mathbf{p},t)$ that describe the creation and annihilation of $\mathrm{eh}$-pairs or a distribution of carries in momenta in the course of one-photon absorption or emission, respectively, have the form:

$$\begin{array}{ccc}& & C_{\gamma }(\mathbf{p},t)=2\int \frac{d\mathbf{k}n\left(k\right)}{{\left(2\pi \hslash \right)}^2}{\int }_{t_0}^{t}d{t}^{\prime }{K}_{\gamma }(\mathbf{p},\mathbf{p}+\hslash \mathbf{k};t,t^{\prime })\hfill \\ & & \left\{f(\mathbf{p},t^{\prime })f(\mathbf{p}+\hslash \mathbf{k},t^{\prime })[1+F(\mathbf{k},t^{\prime })]-[1-f(\mathbf{p},t^{\prime })][1-f(\mathbf{p}+\hslash \mathbf{k},t^{\prime })]F(\mathbf{k},t^{\prime })\right\},\hfill \end{array}$$

(60)

$$\begin{array}{ccc}& & C_{\mathrm{eh}}(\mathbf{p},t)=2\int \frac{d\mathbf{k}n\left(k\right)}{{\left(2\pi \hslash \right)}^2}{\int }_{t_0}^{t}d{t}^{\prime }{K}_{\mathrm{eh}}(\mathbf{p},\mathbf{p}+\hslash \mathbf{k};t,t^{\prime })\hfill \\ & & \left\{f(\mathbf{p},t^{\prime })[1-f(\mathbf{p}+\hslash \mathbf{k},t^{\prime })][1+F(\mathbf{k},t^{\prime })]-f(\mathbf{p},t^{\prime })[1-f(\mathbf{p}+\hslash \mathbf{k},t^{\prime })]F(\mathbf{k},t^{\prime })\right\}.\hfill \end{array}$$

(61)

Collision integrals in the photon sector read:

$$\begin{array}{ccc}& & S_{\gamma }(\mathbf{k},t)=2\int \frac{d\mathbf{p}}{{\left(2\pi \hslash \right)}^2}{\int }_{t_0}^{t}d{t}^{\prime }{K}_{\gamma }(\mathbf{p},\mathbf{p}+\hslash \mathbf{k};t,t^{\prime })\hfill \\ & & \left\{f(\mathbf{p},t^{\prime })f(\mathbf{p}+\hslash \mathbf{k},t^{\prime })[1+F(\mathbf{k},t^{\prime })]-[1-f(\mathbf{p},t^{\prime })][1-f(\mathbf{p}+\hslash \mathbf{k},t^{\prime })]F(\mathbf{k},t^{\prime })\right\},\hfill \end{array}$$

(62)

$$\begin{array}{ccc}& & S_{\mathrm{eh}}(\mathbf{k},t)=2\int \frac{d\mathbf{p}}{{\left(2\pi \hslash \right)}^2}{\int }_{t_0}^{t}d{t}^{\prime }{K}_{\mathrm{eh}}(\mathbf{p},\mathbf{p}+\hslash \mathbf{k};t,t^{\prime })\hfill \\ & & \left\{f(\mathbf{p},t^{\prime })[1-f(\mathbf{p}+\hslash \mathbf{k},t^{\prime })][1+F(\mathbf{k},t^{\prime })]-f(\mathbf{p},t^{\prime })[1-f(\mathbf{p}+\hslash \mathbf{k},t^{\prime })]F(\mathbf{k},t^{\prime })\right\}.\hfill \end{array}$$

(63)

We note that kernels in collision integrals (60)–(63) in the photon and $\mathrm{eh}$-sectors are the same and have the form:

$$\begin{array}{ccc}& & K_{\gamma }(\mathbf{p},\mathbf{p}+\hslash \mathbf{k};t,t^{\prime })=\frac{{\left(e{v}_F\right)}^2}{2\hslash ck}{\mathsf{\Gamma }}_{uv}^{\alpha }(\mathbf{p},\mathbf{p}+\hslash \mathbf{k};t){\mathsf{\Gamma }}_{uv}^{\alpha *}(\mathbf{p},\mathbf{p}+\hslash \mathbf{k};t^{\prime })\cos {\mathsf{\Theta }}^{(+)}(\mathbf{p},\mathbf{p}+\hslash \mathbf{k};t,t^{\prime }),\hfill \\ & & K_{\mathrm{eh}}(\mathbf{p},\mathbf{p}+\hslash \mathbf{k};t,t^{\prime })=\frac{{\left(e{v}_F\right)}^2}{2\hslash ck}{\mathsf{\Gamma }}_{uu}^{\alpha }(\mathbf{p},\mathbf{p}+\hslash \mathbf{k};t){\mathsf{\Gamma }}_{uu}^{\alpha *}(\mathbf{p},\mathbf{p}+\hslash \mathbf{k};t^{\prime })\cos {\mathsf{\Theta }}^{(-)}(\mathbf{p},\mathbf{p}+\hslash \mathbf{k};t,t^{\prime }),\hfill \\ & & {\mathsf{\Theta }}^{(\pm )}(\mathbf{p},\mathbf{p}+\hslash \mathbf{k};t,t^{\prime })=\frac{1}{\hslash }{\int }_{t^{\prime }}^{t}d\tau \left[\epsilon (\mathbf{p},\tau )\pm \epsilon (\mathbf{p}+\hslash \mathbf{k},\tau )-c\hslash k\right].\hfill \end{array}$$

(64)

This fact allows for one to interpret collision integrals (60), (61) and (62), (63) as a reduction of some unified integrands to the photon and $\mathrm{eh}$-sectors.

3. Setting of the Problem

The set of KE (58) and (59) with CI (60)–(63) and Maxwell Equations (35) describe a self-consistent dynamics of carriers in the graphene and a behavior of the internal electromagnetic fields. The same set of equations describes an electromagnetic radiation from graphene, both classical, being generated by plasma currents, and quantum, due to one-photon processes. Because the classical radiation has already been studied2 on the basis of KE, see [61,62], in the present work we focus our attention on the quantum radiation, in particular, comparing its characteristics with the once of the classical radiation. Below, we study this problem solving photon KE (59). This corresponds to neglecting the back reaction of the quantum radiation on carrier dynamics, which is still being described by Equation (25), but with effective external field (1). Thus, the action of the photon subsystem on the evolution of the $\mathrm{eh}$ subsystem and the cascade processes remains outside the scope of this study, as well as the action of the quantum radiation on effective electric field (1).

In order to solve the problem, we rewrite Equations (62) and (63), neglecting the influence of the photon subsystem,

$$\begin{array}{ccc}& & S_{\gamma }(\mathbf{k},t)=\frac{2}{{\left(2\pi \hslash \right)}^2}\int d\mathbf{p}{\int }_{t_0}^{t}d{t}^{\prime }{K}_{\gamma }\left(\mathbf{p},\mathbf{p}+\hslash \mathbf{k};t,t^{\prime }\right)f(\mathbf{p},t^{\prime })f\left(\mathbf{p}+\hslash \mathbf{k},t^{\prime }\right),\hfill \end{array}$$

(65)

$$\begin{array}{ccc}& & S_{\mathrm{eh}}(\mathbf{k},t)=\frac{2}{{\left(2\pi \hslash \right)}^2}\int d\mathbf{p}{\int }_{t_0}^{t}d{t}^{\prime }{K}_{\mathrm{eh}}\left(\mathbf{p},\mathbf{p}+\hslash \mathbf{k};t,t^{\prime }\right)f(\mathbf{p},t^{\prime })\left[1-f(\mathbf{p}+\hslash \mathbf{k},t^{\prime })\right].\hfill \end{array}$$

(66)

Thus, at this stage, we study Equation (59) with CI (65) and (66), which do not contain the photon distribution function. We only take into account a spontaneous emission from the carrier currents produced by an external electric field. Integrals (65) and (66) contain quadratic combinations of the distribution functions of the carriers. These functions should be found as a result of solving an auxiliary problem on the basis of KE (25) or equivalent Equation (27).

4. Studying Processes in Specific External Fields

4.1. Models of External Fields

Below, we study the KE that describes the behavior of the carriers and the electromagnetic field in the graphene when external field is linearly polarized, $A_{\mathrm{ext}}^1=0$ ($E_{\mathrm{ext}}^1=0$). Two following models of the linear polarized short laser pulse [77] are convenient for numerical modelling of the CI (65) and (66). First, we consider the two following potentials $A_{\mathrm{ext}}^2\left(t\right)$ and corresponding fields $E_{\mathrm{ext}}^2\left(t\right)$:

$$\begin{array}{ccc}& & A_{\mathrm{ext}}^2\left(t\right)=A_c\left(t\right)=-\sqrt{\frac{\pi }{8}}{E}_0\tau \exp (-{\sigma }^2/2)\mathrm{erf}\left(\frac{t}{\sqrt{2\tau }}-i\frac{\sigma }{\sqrt{2}}\right)+c.c.,\hfill \\ & & E_{\mathrm{ext}}^2\left(t\right)=E_c\left(t\right)=E_0\exp (-t^2/2{\tau }^2)\cos \omega t,\sigma =\omega \tau ,\hfill \end{array}$$

(67)

and

$$\begin{array}{ccc}\hfill A_{\mathrm{ext}}^2\left(t\right)& =& A_s\left(t\right)=i\sqrt{\frac{\pi }{8}}{E}_0\tau \exp (-{\sigma }^2/2)\mathrm{erf}\left(\frac{t}{\sqrt{2\tau }}-i\frac{\sigma }{\sqrt{2}}\right)-c.c.,\hfill \\ \hfill E_{\mathrm{ext}}^2\left(t\right)& =& E_s\left(t\right)=E_0\exp (-t^2/2{\tau }^2)\sin \omega t,\hfill \end{array}$$

(68)

where $\mathrm{erf}\left(x\right)$ is the error function, $\omega =2\pi /T$ is the cyclic frequency, and $\tau$ is the envelope pulse length.

We are going also to use a model of a harmonic electric field in some analytical calculations:

$$A_{\mathrm{ext}}^2\left(t\right)=A\left(t\right)=-(E_0/\omega )\sin \omega t,E_{\mathrm{ext}}^2\left(t\right)=E\left(t\right)=E_0\cos \omega t.$$

(69)

4.2. Low-Density Approximation

Below, to analyze CI (65) and (66) by analytical methods, we introduce some approximations for the basic constituents of these CI.

In order to estimate of the distribution function, we will use the low-density approximation $f(\mathbf{p},t)\ll 1$ [52], which immediately allows for one to write a solution with the zero initial date $f\left(t_0\right)=0$ of Equations (25) and (26) as:

$$f\left(t\right)=\frac{1}{2}{\int }_{t_0}^{t}d{t}^{\prime }\lambda \left(t^{\prime }\right){\int }_{t_0}^{t^{\prime }}d{t}^{″}\lambda \left(t^{″}\right)\cos \theta (t^{\prime },t^{″}).$$

(70)

It is convenient to rewrite this solution in a slightly different form using definition (26) of the phase and setting the initial time moment (the time of the switching on the external field) to minus infinity ($t_0=-\infty )$,

$$\begin{array}{ccc}& & f\left(t\right)=\frac{1}{2}{\int }_{-\infty }^{t}d{t}^{\prime }{\lambda }_c\left(t^{\prime }\right){\int }_{-\infty }^{t^{\prime }}d{t}^{″}{\lambda }_c\left(t^{″}\right)+\frac{1}{2}{\int }_{-\infty }^{t}d{t}^{\prime }{\lambda }_s\left(t^{\prime }\right){\int }_{-\infty }^{t^{\prime }}d{t}^{″}{\lambda }_s\left(t^{″}\right),\hfill \\ & & {\lambda }_c\left(t\right)=\lambda \left(t\right)\cos \theta (t,-\infty ),{\lambda }_s\left(t\right)=\lambda \left(t\right)\sin \theta (t,-\infty ).\hfill \end{array}$$

(71)

Subsequently, solution (71) can be represented as:

$$f\left(t\right)=\frac{1}{4}{\left[{\int }_{-\infty }^{t}d{t}^{\prime }{\lambda }_s\left(t^{\prime }\right)\right]}^2+\frac{1}{4}{\left[{\int }_{-\infty }^{t}d{t}^{\prime }{\lambda }_c\left(t^{\prime }\right)\right]}^2,$$

(72)

which immediately implies that $f\left(t\right)\ge 0$.

Subsequently, we introduce the approximation of the effective electromagnetic mass while using more accurate estimates of the regularized energy [20]:

$$\begin{array}{ccc}& & \epsilon (\mathbf{p},t)\to {\epsilon }_{*}\left(\mathbf{p}\right)={\epsilon }_{*}\left(p\right)=v_F\sqrt{m_{*}^2{v}_F^2+p^2},\hfill \\ & & m_{*}^2=\frac{e^2}{c^2{v}_F^2}\frac{1}{2T}{\int }_{-T}^{T}dt{A}^2\left(t\right),p=\left|\mathbf{p}\right|,\hfill \end{array}$$

(73)

where $T=2\pi /\omega$ is the period of field oscillations with cyclic frequency $\omega =2\pi \nu$. In the periodic field model (69), this mass reads:

$$m_{*}=\frac{e{E}_0}{\sqrt{2}{v}_F\omega }.$$

(74)

One can present estimates of $m_{*}$ while using parameters that appeared in two experimental studies: $m_{*}=0.03m_e$ for $E_0=3\times {10}^6\mathrm{V}/\mathrm{m}$, $\nu =2\times {10}^{12}\mathrm{H}z$ [62]; $m_{*}=0.047m_e,E_0=2.3\times {10}^8\mathrm{V}/\mathrm{m}$, $\nu =96.7\times {10}^{12}\mathrm{Hz}$ [61].

In the above mentioned approximation, the phase $\theta (t,t^{\prime })$ and amplitude $\lambda (\mathbf{p},t)$ (26) are:

$$\theta (t,t_0)=\frac{2{\epsilon }_{*}}{\hslash }(t-t_0),{\lambda }_{*}(\mathbf{p},t)=eE\left(t\right)l_{*},l_{*}=-v_F^2{p}_1{\epsilon }_{*}^{-2}.$$

(75)

This leads to the fact that only two main harmonics remain in amplitudes (71), ${\mathsf{\Omega }}_{\pm }=2{\epsilon }_{*}/\hslash \pm \omega$, such that:

$${\lambda }_{*c}\left(t\right)=\frac{1}{2}e{E}_0{l}_{*}\left[\cos {\mathsf{\Omega }}_{+}t+\cos {\mathsf{\Omega }}_{-}t\right],{\lambda }_{*s}\left(t\right)=\frac{1}{2}e{E}_0{l}_{*}\left[\sin {\mathsf{\Omega }}_{+}t+\sin {\mathsf{\Omega }}_{-}t\right].$$

(76)

Substituting Equation (76) into Equation (72), we obtain the distribution function for lengthy impulse $\tau \gg T$:

$$f(\mathbf{p},t)={\left(\frac{e{E}_0{l}_{*}}{4{\mathsf{\Omega }}_{+}{\mathsf{\Omega }}_{-}}\right)}^2({\mathsf{\Omega }}_{+}^2+{\mathsf{\Omega }}_{-}^2+2{\mathsf{\Omega }}_{+}{\mathsf{\Omega }}_{-}\cos \omega t).$$

(77)

While taking into account the definition of ${\mathsf{\Omega }}_{\pm },$ we arrive to the representation:

$$f(\mathbf{p},t)=f^{\left(0\right)}\left(\mathbf{p}\right)+f^{\left(2\right)}(\mathbf{p},t),$$

(78)

where the function

$$f^{\left(0\right)}\left(\mathbf{p}\right)=\frac{{(e\hslash E_0{l}_{*})}^2(4{\epsilon }_{*}^2+{\hslash }^2{\omega }^2)}{8{(4{\epsilon }_{*}^2-{\hslash }^2{\omega }^2)}^2}$$

(79)

corresponds to a stationary background distribution and the function

$$f^{\left(2\right)}(\mathbf{p},t)=\frac{{(e\hslash E_0{l}_{*})}^2}{8(4{\epsilon }_{*}^2-{\hslash }^2{\omega }^2)}\cos 2\omega t$$

(80)

corresponds to the breathing mode on the doubled frequency of the external field.

In the general case, the distribution function $f(\mathbf{p},t)$ only contains the even harmonies of external field. This conclusion follows from the general structure of the basic KE (25), (26), or its alternative form (27).

4.3. Calculating Kernels

Let us estimate at first convolutions of the matrix vertices of functions of type (49) entering in kernels (64). To this end, we calculate the functions themselves componentwise, while using definitions of evolution operator (7) and spinors (8),

$$\begin{array}{ccc}\hfill {\mathsf{\Gamma }}_{uu}^1(\mathbf{p},{\mathbf{p}}_1;t)& =& \cos [\varkappa \left(\mathbf{p}\right)/2+\varkappa \left({\mathbf{p}}_1\right)/2]=-i{\mathsf{\Gamma }}_{uv}^2(\mathbf{p},{\mathbf{p}}_1;t),\hfill \\ \hfill {\mathsf{\Gamma }}_{uu}^2(\mathbf{p},{\mathbf{p}}_1;t)& =& \sin [\varkappa \left(\mathbf{p}\right)/2+\varkappa \left({\mathbf{p}}_1\right)/2]=i{\mathsf{\Gamma }}_{uv}^1(\mathbf{p},{\mathbf{p}}_1;t),\hfill \end{array}$$

(81)

where $\varkappa \left(\mathbf{p}\right)=arctan(P^2/P^1)$. It is convenient to write these functions as some algebraic expressions of the quasimomenta $P^1$ and $P^2.$ It should be noted that functions (81) describe the influence of the external field on elementary acts of interaction of electrons and holes (considered to be quasiparticles) with photons. Each of the functions is periodic with the period of the external field, and their convolutions in kernels of CI (62), (63) depend on the observation time t and the antecedent time $t^{\prime }$ that describes the memory effects in the interaction. The sum rules at coinciding times $t=t^{\prime }$ follow from Equation (81):

$$\sum _{\alpha }|{\mathsf{\Gamma }}_{uu}^{\alpha }(\mathbf{p},{\mathbf{p}}_1;t){|}^2=\sum _{\alpha }|{\mathsf{\Gamma }}_{vv}^{\alpha }(\mathbf{p},{\mathbf{p}}_1;t){|}^2=\sum _{\alpha }{|{\mathsf{\Gamma }}_{uv}^{\alpha }(\mathbf{p},{\mathbf{p}}_1;t)|}^2=1.$$

(82)

Here, the influence of the external field is taken into account, but, in comparison with convolutions in kernels (64), the retardation effect is neglected. Such an approximation corresponds to the neglecting harmonics of the external field. The intensity of these harmonics at the fundamental frequency is small in comparison with (82),

$$\frac{e{v}_F^2{p}_{\Vert }{E}_0}{c\omega {\epsilon }_{*}^2\left(\mathbf{p}\right)}\ll 1,$$

(83)

where $p_{\Vert }=\mathbf{p}\mathbf{k}/k$ is the longitudinal momentum. The properties (82) will be used below for estimates of kernels (64) of the CI.

Electromagnetic mass approximation (73) is also effective in estimating phases in kernels of CI (64). In this approximation

$${\mathsf{\Theta }}_{*}^{(\pm )}(\mathbf{p},{\mathbf{p}}_1;t,t^{\prime })=\left[{\epsilon }_{*}\left(\mathbf{p}\right)\pm {\epsilon }_{*}\left({\mathbf{p}}_1\right)-c\hslash k\right](t-t^{\prime }).$$

(84)

Usually, in the absence of external fields, deriving KE (see e.g., [68,69,70]), after integration over the time in CI, one obtains energy conservation laws in elementary acts of scattering of constituents. In the considered highly nonequilibrium situation, phase (84) can be modified by harmonics of the external field due to the time dependence of the carrier distribution function (an example is function (78). It can significantly influence elementary processes in Equation (84). For example, the spontaneous single-photon annihilation channel (sign “+” in Equation (84)) is forbidden at $m_{*}\ne 0$, but the presence of an external field may lead to its opening (process of the stimulated annihilation). Such a situation is also known in the standard strong field QED [78,79]. It corresponds to a general theory of an external field influence on scattering process in strong nonequilibrium systems, see [66,70,78,79,80].

Now, we can advance in calculating time integrals in CI in the leading harmonic approximation under consideration. Taking into account the structure of CI as functionals of the distribution function $f(\mathbf{p},t)$ (78)–(80), phase (84) can acquire additional contributions $n\omega$,$n=0,\pm 2,\pm 4$ from higher harmonics of the external field. Subsequently, integration in time leads to the appearance of singular functions determined by roots of the equation:

$${\phi }^{(\pm )}(\mathbf{p},\mathbf{k};n)={\epsilon }_{*}\left(\mathbf{p}\right)\pm {\epsilon }_{*}(\mathbf{p}+\hslash \mathbf{k})-c\hslash k+n\hslash \omega =0.$$

(85)

Below, these equations are analyzed in some particular cases.

5. Spectral Composition of Quantum Radiation

According to KE (59), the quantum radiation is formed in the annihilation channel with CI (60) and in the channel of momentum redistribution with CI (61). We consider these channels while using approximate CI (65) and (66) in the case of a linearly polarized electric field, neglecting the retardation in the convolutions of the vertex functions with substitutions of the exact convolutions by single-time expressions (82), also taking into account the additional approximations that were discussed earlier in Section 4:

• the low-density approximation for estimation of carrier distribution functions (78)–(80);
• approximation of effective electromagnetic mass (73) in the estimations of the distribution function and phases (64).

We note that both CI (65) and (66) are functionals of different degrees of nonlinearity under distribution function. They have to be calculated in one way in the biharmonic approximation relative to the external field frequency.

5.1. Annihilation Channel

Let us chose time independent and biharmonic parts of CI (65),

$$\begin{array}{ccc}& & S_{\gamma }(\mathbf{k},t)=S_{\gamma }^{\left(0\right)}\left(\mathbf{k}\right)+S_{\gamma }^{\left(2\right)}(\mathbf{k},t),\hfill \end{array}$$

(86)

$$\begin{array}{ccc}& & S_{\gamma }^{\left(0\right)}\left(\mathbf{k}\right)=\frac{{\left(e{v}_F\right)}^2}{2\hslash cka}\int \frac{d\mathbf{p}}{{(2\pi \hslash )}^2}{\int }^{t}d{t}^{\prime }\cos {\mathsf{\Theta }}_{*}^{\left(+\right)}(\mathbf{p},\mathbf{p}+\hslash \mathbf{k};t,t^{\prime })f^{\left(0\right)}\left(\mathbf{p}\right)f^{\left(0\right)}(\mathbf{p}+\hslash \mathbf{k}),\hfill \end{array}$$

(87)

$$\begin{array}{ccc}& & S_{\gamma }^{\left(2\right)}(\mathbf{k},t)=\frac{{\left(e{v}_F\right)}^2}{2\hslash cka}\int \frac{d\mathbf{p}}{{(2\pi \hslash )}^2}{\int }^{t}d{t}^{\prime }\cos {\mathsf{\Theta }}_{*}^{\left(+\right)}(\mathbf{p},\mathbf{p}+\hslash \mathbf{k};t,t^{\prime })\hfill \end{array}$$

$$\begin{array}{ccc}& & \times \left[f^{\left(0\right)}\left(\mathbf{p}\right)f^{\left(2\right)}(\mathbf{p}+\hslash \mathbf{k})+f^{\left(2\right)}\left(\mathbf{p}\right)f^{\left(0\right)}(\mathbf{p}+\hslash \mathbf{k})\right]\cos 2\omega t^{\prime }.\hfill \end{array}$$

(88)

Here, the functions $f^{\left(0\right)}$ and $f^{\left(2\right)}$ are given by Equations (78)–(80).

The time integral in CI (87) with phase (84) (upper sign “+”) implies the following time independent result:

$${\int }_{-\infty }^{t}d{t}^{\prime }\cos {\mathsf{\Theta }}_{*}^{\left(+\right)}(\mathbf{p},\mathbf{p}+\hslash \mathbf{k};t,t^{\prime })=\pi \hslash \delta \left[{\epsilon }_{*}\left(\mathbf{p}\right)+{\epsilon }_{*}(\mathbf{p}+\hslash \mathbf{k})-\hslash ck\right].$$

(89)

Here, the $\delta$-function reflects the energy conservation law of the single-photon annihilation process that cannot be realized at $m_{*}\ne 0$ (see Section 4). In other words, Equation (85) ${\phi }^{(\pm )}(\mathbf{p},\mathbf{k};n=0)=0$ has no physical roots and, therefore:

$$S_{\gamma }^{\left(0\right)}\left(\mathbf{k}\right)=0.$$

(90)

The breathing mode is described by CI (88). Integration here over the time leads to two types of singularities:

$$\begin{array}{ccc}& & {\int }_{-\infty }^{t}d{t}^{\prime }\cos {\mathsf{\Theta }}_{*}^{\left(+\right)}(\mathbf{p},\mathbf{p}+\hslash \mathbf{k};t,t^{\prime })\cos 2\omega t^{\prime }=\pi \hslash \left\{\delta \left[{\phi }^{(+)}(\mathbf{p},\mathbf{k};n=2)\right]\right.\hfill \\ & & \left.+\delta \left[{\phi }^{(+)}(\mathbf{p},\mathbf{k};n=-2)\right]\right\}\cos 2\omega t+\hslash \mathcal{P}\left[\frac{1}{{\phi }^{(+)}(\mathbf{p},\mathbf{k};n=2)}-\frac{1}{{\phi }^{(+)}(\mathbf{p},\mathbf{k};n=-2)}\right]\sin 2\omega t,\hfill \end{array}$$

(91)

where $\mathcal{P}$ is the symbol of principal value. The equations

$${\phi }^{(+)}(\mathbf{p},\mathbf{k};n=\pm 2)=0$$

(92)

describe two singular energy surfaces in the momentum space, which correspond to the two annihilation processes with the participation of higher harmonics of the external field, namely, an emission of an annihilation photon is accompanied by simultaneous radiation $(n=-2)$ or absorption $(n=2)$ of harmonics from reservoir of the external field. One can speak about emission of “soft” or “hard” photons of the quantized field.

The soft photon annihilation process in integral (91) is described by energy condition (92) with $n=-2$. The long-wave approximation is warranted in this region of wave numbers:

$${\epsilon }_{*}(\mathbf{p}+\hslash \mathbf{k})\simeq {\epsilon }_{*}\left(\mathbf{p}\right)+\hslash \mathbf{k}\partial {\epsilon }_{*}\left(\mathbf{p}\right)/\partial \mathbf{p}.$$

(93)

Together with the condition $v_F\ll c$, it leads to the solutions $p^{(1,2)}=\pm p_{*}$ of Equation (92) with $n=-2$, where

$$p_{*}=\hslash (ck+2\omega )/2v_F.$$

(94)

Subsequently, using the representation

$$\delta \left[\phi \left(p\right)\right]=\sum _i\delta (p-p_i){\left(\frac{d\phi \left(p\right)}{dp}\right)}_{p=p_i},$$

(95)

one can specify the $\delta$-function $\delta \left[{\phi }^{(+)}(\mathbf{p},\mathbf{k};n=-2)\right]$ in integral (91) while using the additional approximation $k\ll 2\omega /c$, in order to obtain:

$$S_{\gamma }^{\left(2\right){(-)}^{\prime }}(\mathbf{k},t)=\frac{{\pi }^3{\alpha }^3{c}^2{v}_F^4{E}_0^4}{16{\hslash }^{2}ka{\omega }^7}\left[\frac{3}{4}+{\left(\frac{v_{F}k}{\omega }\right)}^2{\cos }^2\phi \right]\cos 2\omega t,$$

(96)

where $\phi$ is the angle between the vectors $\mathbf{k}$ and $\mathbf{E}\left(t\right)$ characterizing the direction of the radiation.

In the same approximation, one can calculate the principal value integral that corresponds to the function ${\phi }^{(+)}(\mathbf{p},\mathbf{k};n=-2)$ in Equation (91),

$$S_{\gamma }^{\left(2\right){(-)}^{″}}(\mathbf{k},t)=-\frac{5{\alpha }^3{c}^2{v}_F^4{E}_0^4}{32{\hslash }^{2}ka{\omega }^7}ln\frac{2\hslash \omega }{m_{*}{v}_F^2}·\sin 2\omega t.$$

(97)

Thus, the photon production rate in the annihilation channel in the long-wave approximation is equal to the sum of CI (96) and (97),

$$S_{\gamma }^{\left(2\right)}(\mathbf{k},t)=S_{\gamma }^{\left(2\right){(-)}^{\prime }}(\mathbf{k},t)+S_{\gamma }^{\left(2\right){(-)}^{″}}(\mathbf{k},t).$$

(98)

5.2. Channel of the Momentum Redistribution

The channel of the momentum redistribution is described by CI (66). In the low-density limit $f\ll 1$, this CI is a linear functional with respect to the distribution function $f(\mathbf{p},t)$ and can be represented by a decomposition of type (86). Its stationary background part vanishes by virtue of a violation of the energy conservation law ${\phi }^{(-)}(\mathbf{p},\mathbf{k};n=0)=0$.

Let us consider the breathing mode of CI (66). In the approximation under consideration, it reads:

$$S_{\mathrm{eh}}^{\left(2\right)}(\mathbf{k},t)=\frac{{\left(e{v}_F\right)}^2}{{\left(2\pi \right)}^2{\hslash }^{3}cka}\int d\mathbf{p}{\int }_{t_0\to -\infty }^{t}d{t}^{\prime }{f}^{\left(2\right)}(\mathbf{p},t^{\prime })\cos {\mathsf{\Theta }}_{*}^{\left(-\right)}(\mathbf{p},\mathbf{p}+\hslash \mathbf{k};t,t^{\prime }),$$

(99)

where the phase ${\mathsf{\Theta }}_{*}^{(-)}$ is defined by Equation (84). Integration over the time gives the following result:

$$\begin{array}{ccc}& & {\int }_{-\infty }^{t}d{t}^{\prime }\cos {\mathsf{\Theta }}_{*}^{\left(-\right)}(\mathbf{p},\mathbf{p}+\hslash \mathbf{k};t,t^{\prime })\cos 2\omega t^{\prime }=\frac{\pi \hslash }{2}\left\{\delta \left[{\phi }^{(-)}(\mathbf{p},\mathbf{k};n=2)\right]\right.+\hfill \\ & & \left.+\delta \left[{\phi }^{(-)}(\mathbf{p},\mathbf{k};n=-2)\right]\right\}\cos 2\omega t+\frac{\hslash }{2}\mathcal{P}\left[\frac{1}{{\phi }^{(-)}(\mathbf{p},\mathbf{k};n=2)}-\frac{1}{{\phi }^{(-)}(\mathbf{p},\mathbf{k};n=-2)}\right]\sin 2\omega t.\hfill \end{array}$$

(100)

Two equations

$${\phi }^{(-)}(\mathbf{p},\mathbf{k};n=\pm 2)=0$$

(101)

define the singular energy surfaces in integral (100).

The case $n=-2$ corresponds to synchronous emission by a carrier of one radiation photon with energy $\hslash ck$ and the harmonic with doubled frequency of the external field. Such a process is forbidden.

The case $n=2$ (emission of a photon with capture of the harmonic from the external field reservoir) is open and is considered bellow in long-wave approximation (93). Equation (101) for $n=2$ implies the following condition:

$$\pm \hslash k{p}_{\Vert }\frac{v_F^2}{{\epsilon }_{*}\left(\mathbf{p}\right)}+2\hslash \omega =\hslash ck,$$

(102)

where $p_{\Vert }$ is the longitudinal momentum. It is written in the form where the LHS represents a change of the carrier energy in course of the emission of the radiation photon. In order to analyze relation (102), we note that it strongly breaks down at $\omega =0$,

$$\hslash k{p}_{\Vert }\frac{v_F^2}{{\epsilon }_{*}\left(\mathbf{p}\right)}\ll \hslash ck$$

(103)

such that a rather high frequency is necessary for the restoration of Equality (102). Because the LHS of Inequality (103) is very small in comparison with RHS, one can conclude that the energy of the radiated photon and the energy $2\hslash \omega$ from the reservoir of an external field are correlated, so that $ck\ge 2\omega$. It means that the radiated photon captures most of the energy from the external electromagnetic reservoir.

Now, we can introduce the detuning

$$\mathsf{\Delta }=\frac{c}{v_F}\left(1-\frac{2\omega }{ck}\right)\ge 0$$

(104)

and to write a solution of Equation (102) in the linear approximation relative to the longitudinal momentum $p_{\Vert }$ in the following form:

$$p_{\Vert }=\pm p_{\Vert }^{*},p_{\Vert }^{*}=\frac{{\epsilon }_{*\perp }\left(p_{\perp }\right)}{v_F}\mathsf{\Delta },$$

(105)

where ${\epsilon }_{*\perp }\left(p_{\perp }\right)=v_F\sqrt{m_{*}^2{v}_F^2+p_{\perp }^2}$ is the transversal energy.

Let us rewrite CI (99) keeping in Equation (100) the contribution with the $\delta \left[{\phi }^{(-)}(\mathbf{p},\mathbf{k};n=2)\right]$ only,

$$S_{\mathrm{eh}}^{{\left(2\right)}^{\prime }}(\mathbf{k},t)=\frac{e^4{E}_0^2{v}_F^6}{64\pi cka}\int d^{2}p\frac{p_1^2}{{\epsilon }_{*}^4\left(\mathbf{p}\right)}\frac{\delta \left[{\phi }^{(-)}(\mathbf{p},\mathbf{k};n=2)\right]}{4{\epsilon }_{*}^2\left(\mathbf{p}\right)-{\hslash }^2{\omega }^2}\cos 2\omega t.$$

(106)

Here, the relation

$$\delta \left[{\phi }^{(-)}(\mathbf{p},\mathbf{k};n=2)\right]=\frac{{\epsilon }_{*\perp }\left(p_{\perp }\right)}{k{v}_F^2}\left[\delta (p_{\Vert }-p_{\Vert }^{*})+\delta (p_{\Vert }+p_{\Vert }^{*})\right]$$

(107)

which is valid in long-wave approximation (93) was taken into account. After integration over $p_{\Vert }$, it is convenient to write the residual integral using the dimensionless variable $x={\epsilon }_{*\perp }\left(p_{\perp }\right)/m_{*}{v}_F^2(x\gg 1)$,

$$\begin{array}{ccc}& & S_{\mathrm{eh}}^{{\left(2\right)}^{\prime }}(\mathbf{k},t)=\frac{\pi c\hslash {\alpha }^2{E}_0^2}{16k^{2}a{[1+{\mathsf{\Delta }}^2]}^{5/2}{m}_{*}^2{v}_F^3}{\int }_1^{\infty }\frac{dx}{x^2{(x^2-1)}^{1/2}}\hfill \\ & & \times \frac{{\mathsf{\Delta }}^2{x}^2{\cos }^2\phi +(x^2-1){\sin }^2\phi }{x^2-{\xi }^2}\cos 2\omega t,\xi =\frac{\hslash \omega }{2m_{*}{v}_F^2\sqrt{1+{\mathsf{\Delta }}^2}}.\hfill \end{array}$$

(108)

If we suppose in addition that the inequality $2\omega >ck$ (which that is equivalent to $\left|\mathsf{\Delta }\right|\gg 1)$ holds true, CI (108) will be equal:

$$\begin{array}{ccc}& & S_{\mathrm{eh}}^{{\left(2\right)}^{\prime }}(\mathbf{k},t)=\frac{\pi c\hslash }{16ka\omega }{\left(\frac{\alpha E_0}{m_{*}{v}_F}\right)}^{2}I\left(\tilde{\xi }\right){\cos }^2\phi ·\cos 2\omega t,\hfill \end{array}$$

(109)

$$\begin{array}{ccc}& & I\left(\tilde{\xi }\right)={\int }_1^{\infty }\frac{dx}{\sqrt{x^2-1}}\frac{1}{x^2-{\tilde{\xi }}^2}=\frac{arcsin\left(\tilde{\xi }\right)}{\tilde{\xi }\sqrt{1-{\tilde{\xi }}^2}},\tilde{\xi }=\hslash k/4m_{*}{v}_F.\hfill \end{array}$$

(110)

The residual part of CI (99) is a principal-value integral that corresponds to the contribution $\mathcal{P}{\left[{\phi }^{(-)}(\mathbf{p},\mathbf{k};n=2)\right]}^{-1}$ in Equation (100),

$$S_{\mathrm{eh}}^{{\left(2\right)}^{″}}(\mathbf{k},t)=\frac{{\alpha }^2\hslash c{E}_0^2{v}_F^4}{4k^{2}a}\int \frac{d^{2}p}{{\epsilon }_{*}^3\left(\mathbf{p}\right)}\frac{p_1^2}{4{\epsilon }_{*}^2\left(\mathbf{p}\right)-{\hslash }^2{\omega }^2}\mathcal{P}\left[\frac{1}{p_{\Vert }-p_{\Vert }^{*}}+\frac{1}{p_{\Vert }+p_{\Vert }^{*}}\right]\cos 2\omega t.$$

(111)

Here, $p_1=p_{\Vert }\cos \phi +p_{\perp }\sin \phi$, where $p_{\Vert }$ and $p_{\perp }$ are the longitudinal and transversal components. In deriving Equation (111), it was taken into account splitting (102) of the energy surface ${\phi }^{(-)}(\mathbf{p},\mathbf{k};n=2)=0$. In this form, it is easy to see that CI (111) is equal to zero by virtue of the oddness of the integrand function relatively substitutions $p_{\Vert }\to -p_{\Vert }$ and $p_{\perp }\to -p_{\perp }$, such that:

$$S_{\mathrm{eh}}^{{\left(2\right)}^{″}}(\mathbf{k},t)=0.$$

(112)

5.3. Comments on the Results Obtained

Representations of CI (96)–(98) and (108) allow for us to make the following comments and conclusions:

(1) The channel of radiation based on the mechanism of the momentum redistribution is more effective in comparison with the annihilation channel. In order to explain this conclusion, it is necessary to note that the properties of the eh-plasma are defined by the distribution function $f(\mathbf{p},t)$. In our study, this function is calculated in the low-density approximation, where, according to Equations (78)–(80), $f(\mathbf{p},t)\sim \alpha {E_0}^2$. One can see that CI (65) in the annihilation channel is a quadratic functional with respect to $f(\mathbf{p},t)$, while CI (66) in the momentum redistribution channel is the linear functional (in the low-density approximation) with respect to $f(\mathbf{p},t)$. It is stipulated by distinctions in these mechanisms of radiation: the emission of photons in the case of the momentum redistribution takes place with participation of the quasiparticles from one of the subsystems (subsystem of electrons or holes), whereas radiation in the annihilation channel is in need in two partners from different subsystems. Thus, in the leading approximation, we have: $S_{\gamma }/S_{\mathrm{eh}}\sim \alpha$.

(2) A time independent component is absent in the spectrum of the radiation of the quantized electromagnetic field. In the approximation under consideration, it is stipulated by the absence of the energy feeding of the corresponding single-photon processes from the photon reservoir of the external field.

The basic breathing harmonic of the radiation is the doubled harmonic of the external field. It is basic frequency of oscillations of the eh-plasma (see Equation (80)). The presence of odd harmonics in the spectrum of quantized radiation has principal importance for an experimental identification of the quantum radiation since the competitive quasiclassical radiation on the frequency of the plasma oscillations contain odd harmonics only [61,62,65].

(3) Representations for CI (96)–(98) and (108) were obtained in the long-wave approximation and describe the situation well in the area of small wave numbers. At the same time, CI (96) and (97) of the annihilation channel demonstrate $1/k$-dependence in an explicit form with feebly marked anisotropy. In the standard QED, such $1/k$ behavior was predicted long ago in References [57,67]. The k-dependence in the momentum redistribution channel is more complicated, it corresponds to the emission of electrons and holes that are accelerated by the electric field in the opposite directions (see Equation (105)).

(4) The quantized electromagnetic radiation expands in the graphene plane and it is anisotropic: strong anisotropic effect in the momentum redistribution channel overlaps on isotropic radiation in the annihilation channel. We note that the quasiclassical radiation of the plasma oscillations in the graphene propagates always in the perpendicular direction to the graphene plane.

(5) Figure 1 illustrates some of these features of the quantized field radiation. Here, the spectral composition of the total photon production rate in the long waves region is represented. The central isotropic peak corresponds to the $1/k$-dependence. This feature in the radiation spectrum corresponds to the annihilation mechanism.

The peripheral strong anisotropic distribution is a result of a momentum redistribution mechanism in CI (108) for $1.95\times {10}^{-4}\le \tilde{\xi }\le 0.456$. Such fixation of the parameter guarantees finding the parameters $E_0,k$ and $\omega$ limited by the usability condition (104) in the physical region and it provides sufficient distance from the resonance point $\tilde{\xi }=1$ (see, e.g., Equation (110)). Figure 1 shows that photons are predominantly created in the direction of the acting external field.

$1/k$-dependence in the region of the central peak is reproduced also on the qualitative level by virtue of numerical calculations directly of CI (65) in the annihilation channel.

Specific of the computer calculations forces one to use the field model of a short laser impulse (68) (see Figure 2 for the case of linear polarization, $E_0={10}^6\mathrm{V}/\mathrm{cm},\omega =2\pi \times {10}^{14}\mathrm{Hz},\sigma =3$). The distribution function of the eh-subsystem calculated on the basis of basic KE (27) is depicted for the point of time $t=0$ on Figure 3. The valley $p_1=0$ between two symmetrical fragments of the distribution function is a consequence of the structure of amplitude (10) ($E_{\mathrm{ext}}^1\left(t\right)=0$ by virtue of the chosen polarization of the external field). Finally, Figure 4 shows $1/k$-dependence in projection on the axis of $k_2$ calculated on the basis of CI (65) of the annihilation channel.

(6) As it was shown in Section 5, the considered above CI in the model of periodical field (69) do not contain time independent components. This fact can be explained by the absence of a definite asymptotic limit of field model (69) at $t\to \infty$ [81]. In the general case CI can have non-zero asymptotic at $t\to \infty$ if the back reaction of the system is not taken into account. Figure 5 demonstrates such situation in the model of a short laser pulse (67). We believe that special consideration of such behavior would be quite important.

6. Conclusions

In the present work, we have considered some physical processes that are possible in the graphene under action of a strong time dependent electric field. To this end, nonperturbative methods (see [13,14,15] and relevant References given in the text of the article) of strong field QED with unstable vacuum (in particular, the Dirac model of the graphene) were used in combination with kinetic description of radiation from the electron-positron plasma created from the vacuum under an action of a strong time dependent electric field that was developed in References [57,58,66,67,80]. We would like to emphasize the significant development of KE formalism and study of the structure of these equations presented in the work. This formalism includes a nonperturbative basis oriented on the description of excitations of the eh-plasma in the graphene under the action of strong electric fields and a perturbative part described interaction with the quantized electromagnetic field. Section 5.3 presents the main concrete results obtained in the work on the base of such generalized theory. We stress that some of the predicted properties of the model under consideration may be verified experimentally. First, is the issue of the possible presence of even harmonics of the external field in the quantum radiation spectrum. Another important property that can be tested is the characteristic spectral composition anisotropy of the quantum radiation and its direction in the graphene plane (emission in the external space of the electromagnetic waves of the plasma oscillations is oriented perpendicular to the graphene plane). Some other predictions are related to the long wave features of this radiation. In this respect, we note that the developed approach can be extended to a wider range of parameters of the external field, including into the consideration arbitrary polarizations. It is also important to consider in a similar manner such processes as photoproduction of the eh-plasma under the action of quantum electromagnetic field, cascade processes, the Breit–Wheeler process (e.g., [82]) and so on. We recall that a radiation of quasiclassical plasma waves in the graphene obtained in the frame work of the nonperturbative KE approach had got sufficiently reliable experimental testing (e.g., [65] and having there references).

In the conclusion, it should be noted that the new extended formulation of KE approach contains some model assumptions. Further calculations in the framework of the formulation may check the validity of these assumptions. Such checking would be important for applications of strong field QED in $3+1$ dimensions, where the situation is more complicated.