Hydrodynamical model of QED cascade expansion in an extremely strong laser pulse

Development of the self-sustained quantum-electrodynamical (QED) cascade in a single strong laser pulse is studied analytically and numerically. The hydrodynamical approach is used to construct the analytical model of the cascade evolution, which includes the key features of the cascade observed in 3D QED particle-in-cell (QED-PIC) simulations such as the magnetic field predominance in the cascade plasma and laser energy absorption. The equations of the model are derived in the closed form and are solved numerically. Direct comparison between the solutions of the model equations and 3D QED-PIC simulations shows that our model is able to describe the complex nonlinear process of the cascade development qualitatively well. The various regimes of the interaction based on the intensity of the laser pulse are revealed in both the solutions of the model equations and the results of the QED-PIC simulations.


I. INTRODUCTION
There is strong evidence that the development of QED cascades is an immanent feature of interaction of an extremely strong electromagnetic fields with matter in a majority of configurations [1][2][3][4][5][6][7][8][9][10][11] . The essence of the cascade is emission of the high-energy photons by the ultrarelativistic particles (the nonlinear Compton scattering) and the subsequent decay of these photons into the electron-positron pairs (the Breit-Wheeler process) which leads to the multiplication of particles. These processes are believed to play an important role in many astrophysical phenomena like cosmic ray showers 12 , processes in pulser magnetosphere 1,2,13 , γ-ray bursts 14 etc. The variety and complexity of the e + e − plasma structures produced in QED cascading makes it clear that there is no simple way to tackle the problem. One of the reasons behind it is that the emission of the photons by the electrons and the positrons greatly alters the dynamics of the latter -the effect known as the radiation reaction. Accurate description of the radiation reaction is a long standing problem of both classical and quantum electrodynamics [15][16][17][18] .
With the upcoming multi-PW level laser facilities such as ELI 19 and Apollon 20 these processes are expected to be observed in the laboratory light-matter interaction experiments. An extensive search for the optimal configuration of such experiments is being conducted nowadays [21][22][23][24][25][26][27] . As already mentioned the QED cascade is a significantly nonlinear phenomenon that is why its analytical study is complicated. And while PIC simulations serve as a starting point for most of the theoretical research and can give a valuable insight into the nature of the discussed processes, even deriving phenomenological laws or scaling can be extremely time consuming as it usually requires scanning over the multi-dimensional map a) Electronic mail: asams@ipfran.ru of the parameters. These dependencies though are crucial to design the experiments with the next generation of laser facilities. Various schemes are proposed in order to lower the threshold of QED cascading: multiple laser pulses with small number of seed particles, laser-beam interaction, etc. 8,24,25,[28][29][30][31][32][33][34] The key reason for such configurations to be optimal is that the governing parameter of the QED processes χ which is the ratio of the transverse component of the effective electric field experienced by the relativistic particle in the rest frame to the critical Sauter-Schwinger field 35 is maximaized in these scenarios. There is a crucial difference though between these two approaches. In the first case the cascade gains its energy from the electromangetic field and this is so-called A-type cascade, while in the second one the cascade energy is mostly limited by the initial energy of the seed and this is called the S-type cascade 36 .
A specific configuration of QED cascade development which occurs in the interaction of the extremely intense laser pulse with the solid target is recently explored by 3D QED-PIC simulations 37 . The peculiar mechanism of that cascade development makes it difficult to attribute it to either of A-type or S-type cascade. For the sake of convenience we will briefly describe the core mechanism of the discussed cascade (see Fig. 1 for the visual schematics of the process). The main feature that allows the cascade to sustain itself is the fact that the collective motion of the electrons and the positrons alters the laser field so that its magnetic component becomes larger than the electric one while staying mutually nearly perpendicular. In such fields the particles drift along the direction of the laser propagation with the drift velocity and rotate in the plane perpendicular to the direction of the magnetic field, so their trajectories are helical (see Supplemental material 38 ). The particles radiate gamma-quanta along the direction of the instant tangent to the trajectory. Occasionally that direction may be opposite to the direction of the laser propagation. Such gamma-quanta eventually leave the electron-positron plasma and get to the vacuum where strong laser field is present and where they are highly probable to photoproduce new electronpositron pairs. The newly created pair is then accelerated by the laser pulse and pushed towards the plasma region and the process repeats. As a results the QED cascade continuously expands towards the laser pulse building a pairs plasma 'cushion'. This proccess is similar to the propagation of the ionization front in the microwave gas discharge 39 The models describing electrodynamics of the cushion plasma 37,40 and QED cascade evolution 37 are proposed. However these models are not self-consistent. The model proposed in Ref. 40 does not take into account particle multiplication because of cascading while the model in Ref. 37 neglects the laser field depletion because of absorption in the cascade plasma. In this paper we construct the self-consistent model which describes spatio-temporal evolution of both the laser field and the cascade plasma. The goal of this paper is twofold. First, we develop a reduced model which requires much less computational resources than that for 3D QED-PIC simulations. Second the model development and its verification allows us to understand and evaluate the role of different physical processes behind QED cascading in a single laser pulse.
The paper is organized as follows. In Sec. II we formulate the problem and start from the kinetic equations including the QED processes. Next we propose a set of general simplifications which allow to greatly reduce the complexity of the equations. According to these simplifications self-consistent hydrodynamics equations are derived. In Sec. III the method of the numeric solution of the derived equations is discussed and the results of the solution are compared to the results of the QED-PIC simulations. The obtained results are summarized and discussed in Sec. IV. In Appendix A we examine the problem of the electron acceleration in the plane wave and in Appendix B we derive the electrodynamical properties of the dense electron-positron plasma. Results of both these problems are used in the model equations derivation.

II. QED CASCADE DEVELOPMENT MODEL DERIVATION
Similar to Refs. 28,41 we start our analysis from the kinetic equations for electrons, positrons and gammaquanta, assuming that the QED cascade is in the selfsustaining stage so the seed particles (e.g. the electrons and the ions of the target) do not contribute to it. The kinetic equations and the Maxwell's equations can be written as follows where f e ± ,γ (t, r, p) are the distribution functions of the electrons, the positrons and the gamma-quanta respectively, v is the particle velocity which is equal to p/ 1 + p 2 for the electrons and the positrons and to p/p for the gamma-quanta, w rad (p , p)dp is the probability in the time unit for the electron or the positron with the momentum p to emit the gamma-quant with the momentum p − p, w pair (p , p)dp is the probability in the time unit for a gamma-quant with the momentum p to photoproduce the electron with the momentum p and the positron with the momentum p − p. We use a common relativistic normalization where the electric and the magnetic fields are normalized to the value of m e cω L /e, where m e and e > 0 are the mass and the charge of the electron, c is the speed of light and ω L is the characteristic frequency of the external field, the particle number densities are normalized to the critical density n c = m e ω 2 L /4πe 2 , the energies and the momenta are normalized to m e c 2 and m e c respectively, the coordinates and the time are normalized to c/ω L and 1/ω L respectively, thus the velocities are normalized to the c.

A. Model assumptions
Let us apply several assumptions to simplify the model. First of all, as we investigate interaction with the plane EM wave, the problem can be considered spatially one-dimensional. If we also restrict ourselves to the circularly-polarized laser pulses then the symmetry relative to the rotation along the axis of the pulse propagation can be also utilized. This simplifications lead to the distribution functions becoming dependant on three variables (excluding time) rather than six: f (t; r, p) = f (t; x, p, θ)/2π, where p is the momentum of the particle and θ is the angle between the momentum and the x-axis.
Secondly, we will assume all distribution functions to be locally monoenergetic, i.e. f ∝ δ(p − p(x))/p 2 , where p(x) is the mean value of the momentum of the particles located in the small vicinity of x. We denote the mean energy of gamma-quanta as ε γ and the mean energy of pairs as ε p , assuming that particles are ultrarelativistic and thus ε 2 p = 1 + p 2 p ≈ p 2 p . While the monoenergetic approximation is quite strong, we suppose that the mechanism of the cascade development explained in Sec. I does not rely on any particular feature of the particles energetic spectra. So we argue that accounting energetic spectra evolution in our model will cause only quantitative changes rather than qualitative ones while greatly complicating the equations. It will be discussed later that this assumption is valid for the pairs which enter the plasma region with approximately equal energies. As for gammaquanta we actually use a two-stream approximation, i.e. we separate gamma-quanta into ones that are emitted in the vacuum region and propagate mostly along the direction of the laser pulse propagation (x-axis) and thus do not contribute to the cascade (we designate them as decoupled gamma-quanta) and ones that are emitted in the plasma region in many different directions and provide a positive feedback needed for the cascade development (we designate them as either involved gamma-quanta or simply gamma-quanta). As Fig. 2 shows the energy spectrum of the gamma-quanta is broad, while if we exclude the decoupled gamma-quanta, then the width of the spectrum diminishes significantly which justifies our assumption. As the decoupled gamma-quanta affect the cascade development only by taking away some portion of the total energy, their spatial distribution is irrelevant for the cascade but it will be calculated for more explicit comparison with the results of the QED-PIC simulations.
In order to omit integration over energies and azimuth angle ϕ we redefine f as follows where n(x) is the particles density distribution, and Φ(θ) Figure 2. The mean energy εγ of the gamma-quanta located in the vicinity of the x-coordinate calculated from all the particles (red line) and from the particles with the velocity along the x-axis less than 0.5 which supposedly include only the involved gamma-quanta (green line). The error bars depict the standard deviation. The data is taken from the results of the PIC simulation for the time instance ct/λ = 18. The simulation parameters are discussed in Sec III. The initial conditions are the same as in Fig 6. is the particles momenta angular distribution such that π 0 Φ(θ) sin θdθ = 1.
where N is the total number of particles. Assuming monoenergetic distribution hints that we can use a hydrodynamical approach by calculating moments of the distribution functions from Eqs.(1)-(2). The following quantities can be introduced W pair (χ γ , ε γ ) = w pair (p, p )dp , W rad (χ p , ε p ) = w rad (p, p )dp , where W pair , W rad , I rad are the total probabilities of the pair photoproduction (or gamma-quant decay), gammaquanta emission and the intensity of the gamma radiation, respectively 18 . Note that these quantities depend on the Lorentz-invariant QED parameter χ where ε is the particle energy, E S = eE S /m e cω L = m e c 2 / ω L and E S = m e 2 c 3 / e is the critical Sauter-Schwinger field 35 .
Generally the hydrodynamics equations have form of a continuity equation, i.e.
where D α and F α are the density and the flux of some value α and S[α, β] is the source responsible for change of the value α in the process β. Note that despite the fact that we specified the energy distribution of the particles, to calculate the sources S[α, β] one also needs to know the angular distribution of the particle which will be discussed below. We suppose that the following set of the equations can quantitatively describe the cascade development ∂ ∂t where n p = n e + = n e − is half the density of the electron-positron plasma assuming that cascade plasma is quasineutral, v x is the mean velocity of the pairs which is calculated below and v γ and v 2 γ are the mean longitudinal velocity and mean square longitudinal velocity of the gamma-quanta which are calculated from their angular distribution as follows Eq. (18) Fig. 1]. Σ γ is the total energy of the decoupled gammaquanta. The factor ψ vac (ψ pl ) is equal to 1 in the vacuum (plasma) region and is equal to 0 in the plasma (vacuum) region. These factors will be specified below. Note that ψ vac + ψ pl = 1. For the sake of convenience these factors will be omitted where it is clear which region is being discussed.

B. Electromagnetic field configuration
In the vacuum region the fields are close to the plane wave, i.e. the magnitudes of the electric and the magnetic fields are equal E ≈ B and they are mutually perpendicular E · B ≈ 0. According to the results of the 3D QED-PIC simulations the electric and the magnetic fields inside the plasma region stay almost mutually perpendicular and magnetic field magnitude is everywhere larger than the electric one B > E. The spatial distribution of the electromagnetic field has a characteristic scale of λ in both regions. In such field configuration the charged particle drifts perpendicular to the both the electric and the magnetic field with the velocity Assuming that the laser pulse propagates along the xaxis and thus the fields lay in the yz-plane this velocity is directed along the x-axis. Assuming that the fields are mutually perpendicular we get In the vacuum region the particles move in the intense EM plane wave with equal electric and magnetic fields E = B. If the initial energy of the particle ε is smaller than the field amplitude E then on the timescale much smaller than the laser period the particle longitudinal velocity tends to the speed of light, i.e. v x ≈ 1 = E/B (see Appendix A). So we assume that Eq. (22) is valid in both the vacuum and the plasma region.

C. Coupled gamma-quanta distribution function
As discussed in Sec. I the involved gamma-quanta are emitted by the pairs during their motion along the helical trajectories in the plasma region (see also Supplemental material 38 ). Because of that the angular distribution of the gamma-quanta is quite wide. We will assume that it is smooth and can be described by the single parameter. This parameter is the velocity of the instantaneous reference frame in which the angular distribution of the gamma-quanta located in the vicinity of the x-coordinate is uniform. Thus in the laboratory reference frame the total distribution function of the involved gamma-quanta has the following form The mean velocity v γ and the mean square velocity v 2 where ath(x) is the inverse hyperbolic tangent function. The results of the QED-PIC simulations show that Eq. (24) is a good approximation for the angular distribution of the involved gamma-quanta [see Fig. 3 (a), The χ parameter of the gamma-quanta propagating in the crossed electric and magnetic fields which is the case for both the plasma and the vacuum region can be calculated as follows where we used (22).
Having specified the distribution function of the gamma-quanta completely it is possible to calculate sources S[α, pp] which correspond to the process of the pairs photoproduction

D. Pairs dynamics in the vacuum region
Let us first examine the electrons and the positrons located in the vacuum region, where the number of particles is small so the collective plasma effects can be neglected. Thus the electromagnetic field in that region coinsides with the field of the incident radiation. In our case these electrons and positrons that are born in the vacuum region move in the field of a plane EM wave. Dynamics of a single particle in the plane wave is discussed in Appendix A. Computing the sources S[α, β] in the rhs of the (13) -(18) requires knowledge of the particles distribution function. Although approximate particles trajectories can be found analytically, deriving explicit expression for the distribution function is practically unfeasible due to the fact that particles are being born at different time instances with varying initial conditions. But few observations will allow us to estimate S[α, β] using the different approach. The first observation is that a strong EM plane wave pushes particles in the direction of its propagation, i.e. the x-axis. So after some time interval independently on its initial direction the particles momenta orient parallel to the x-axis, so that v x ≈ 1 [see Fig. 1 (c)]. We neglect the duration of this momenta orientation, which allows us to approximate the fluxes of the particles density and particles energy density by simply multiplying these densities by the velocity v x ≈ 1.
The purpose of the continuity equations in the vacuum region essentially is to provide the values of the particles densities, energies and magnitude of the electric field at the plasma-vacuum interface (which we call occasionally the cascade front). So we are interested in the total contribution to these values by each particle during its motion from the moment of its birth in the vacuum region up to the moment of reaching the plasma boundary. We do this by defining the sources S[ε, β] in the following way where ∆ε β is the total change of the energy due to the process β of a single particle born at the point in space with the coordinate x at the time instance t while the particle stays in the vacuum region.
The energy gained by the single particle in the plane wave can be evaluated as where ε 0 is the initial particle energy and µ is the parameter determining the time the particle spends in the vacuum region (see Appendix A). It depends on the time of the particle birth and the time when the particle crosses the plasma boundary. Finding the latter requires either building an independent model of the cascade front propagation or finding some heuristics both of which lays out of the scope of the paper. Instead we suppose that µ is constant during the whole cascade development and its value can be approximated by comparing the results of the model with the results of QED-PIC simulations, so µ is the first fitting parameter of our model. Noting that when pairs are born from gamma-quanta their average initial energy is approximately equal to half the energy of the parent gamma quanta ε γ we derive the expression The total energy loss due to photon emission of a single particle can be calculated as follows (see Appendix A) So the source term S[ε p , rad] takes the form Relation between the velocity vector and the magnetic field vector in the reference frame K moving along the x-axis with the velocity vx = E/B.

E. Pairs dynamics in the plasma region
In the plasma region at each point in space there exists an instantaneous reference frame K moving with the velocity v x (x, t) ≈ E/B in which only magnetic field presents. It is convenient to obtain some results in that reference frame. In K frame electrons and positrons move along the magnetic field with the velocity v B and revolve in the plane perpendicular to it with the velocity v ⊥ (see Fig. 4). We will assume that the particles remain ultrarelativistic in that reference frame, so v ⊥ 2 +v B 2 ≈ 1. Moving along the magnetic field results in non-zero average current which has to be accounted in the Maxwell's equations, while revolving in the magnetic field gives zero average current but it is responsible for producing gamma-quanta. In particular the Lorentz-invariant QED parameter χ for pairs can be written as We express the primed values through the values in the laboratory reference frame as follows: where we used the fact that the average particle momentum along the x-axis is equal to γv x and v x = E/B. The final form of the expression for χ is Due to particles revolving around the direction of the magnetic field we can assume that the angular distribution of the particles in the K reference frame is close to uniform. Applying Lorentz transformation we derive the following distribution function of pairs in the laboratory reference frame where Φ is defined the same as in Eq. (24) Same as for the gamma-quanta the results of the QED-PIC simulations demonstrate that Eq. (40) is a good approximation for the angular distribution of the pairs [see Fig. 3 (c), (d)]. Below we will show that the velocity v x can be calculated from the local values of the electric field and plasma density. That is why we do not include the continuity equation for the density of the longitudinal velocity of the pairs similar to Eq. (16). Also in the case of the pairs we choose to neglect the difference between the velocity v of the reference frame in which the particles angular distribution is uniform and the actual mean velocity v calculated from that distribution, the maximum difference between which is less then 0.2, according to Eq. (25). Because χ p does not depend on θ we calculate the sources S[α, rad c ] as follows The total current density of the particles averaged over the characteristic period of the Larmor oscillations τ B = ε p /B is The factor 2 is from the fact that currents of the electrons and the positrons are co-directional. This comes from the observation that in the laboratory reference frame the electric and the magnetic fields are actually not exactly perpendicular. It means that in the reference frame K there exist a small electric field directed along or opposite to the magnetic one (depending on the sign of the E · B product). Presence of that field leads to the average electrons velocity being counter-directed to it and the average positrons velocity being co-directed to it. This is not the case for the longitudinal motion of the pairs which does not depend on the charge sign, thus currents of the electrons and the positrons eliminate each other along the x-axis but sum up in the yz-plane. This fact also hints that the electron-positron plasma is actually a conducting media thus some electromagnetic energy absorption also occurs in that region, though it is significantly smaller than the absorption in the vacuum region which is seen in the QED-PIC simulations [see Fig. 5 (a)] and thus we do not account it in our model. The value of v B averaged over the particles which we denote as ν is the second fitting parameter of our model. We can roughly estimate it by noting that for a single particle the value of v B cannot exceed its initial value during the particle motion. The particles enter the plasma region after being accelerated by the laser pulse with predominantly longitudinal velocity, i.e. velocity along the x-axis, so the initial projection of the particles velocity onto the magnetic field which lies in the yz-plane is small. So we expect our model to give valid results with the values of ν closer to zero rather than unity.
The electrodynamical properties of the medium which response to the plane EM wave consists of inducing the current along the magnetic field is examined in Appendix B. The main conclusion is that the ratio between the electric and the magnetic field in such media can be expressed via its density and the electric field amplitude Validity of Eq. (46) is verified by direct comaprison against the mean values of the longitudinal velocity of the particles computed in the PIC simulation as show in Fig. 5 (b).

III. MODEL FORMULATION AND COMPARISON WITH QED-PIC SIMULATIONS
The last undefined terms are the ψ vac and ψ pl which determine the vacuum and the plasma regions, respectively, in space. We note that the longitudinal velocity v x defined in Eq. (46) actually demarcates these regions: in the vacuum region v x ≈ 1 and in the plasma region v x < 1. So we choose ψ vac and ψ pl in the following way where M ∼ 10 is a constant. We choose the exact value of this parameter by defining the upper threshold for v x above which we assume plasma to be rarefied enough not to cause any collective effects. So we choose this threshold value to be 0.7 and M = 8.
The final set of the cascade model equations is the following ∂ ∂t (ε p n p ) + ∂ ∂x (v p ε p n p ) = W pair n γ ε γ 2 + µE 2/3 ε 1/3 γ G rad − I vac n γ ψ vac − I pl n p ψ pl , ∂ ∂t It is important to note that the energy is conserved in the model, i.e.
The main approximations -Eqs. (33) and (46) -are directly validated based on the results of the QED-PIC simulations (see Fig. 5). PIC simulations were performed using the QUILL code 42 , which enables modelling of the QED effects via the Monte-Carlo method. The initial distrubution of the EM fields has a form of the plane wave with the wavelength λ = 2πc/ω L = 1 µm and the amplitude a 0 propagating along the x-axis with the spatiotemporal envelope given by the following expression a(x, y, z) = cos 2 π 2 x 4 σ 4 x cos 2 π 2 The transverse spatial size of the laser pulse is 2σ r = 18 µm and the pulse duration is 60.5 fs (2σ x = 18.15 µm).
The simulation box size is 30λ × 30λ × 30λ, the grid size is 3000 × 300 × 300. As discussed in Ref. 37  gle laser pulse almost does not depend on the seed, so we choose the seed in the form of a short gammabunch counter-propagating to the laser pulse in order not to introduce the electron-ion plasma to the interaction which is significantly different from the forming electron-positron plasma. The initial seed of that form can be used in our model by initializing v γ (t = 0) ≈ −1.
The density distributions in our model and PIC simulation coinside and are expressed by the following formula: where w γ is the half-width of the bunch and x 0 is the position of its center.
The Eqs.(49)-(54) are solved numerically using the method of lines (MOL): partial derivatives ∂/∂x are approximated with the finite differences to derive the system of the ODEs which is solved using the explicit Runge-Kutta method. This scheme is not internally conservative so the energy conservation is done manually at each integration step by clipping the derivative ∂n γ /∂t so that the total energy does not grow. The relative error gained from that procedure turns out to be acceptably small.
Direct comparisons between the solutions of Eqs. (49) -(54) and the results of the QED-PIC simulations are shown in Figs. 6 -8. Our model coinsides with the results of the full QED-PIC simulations qualitatively good in terms of the distributions of the particles and the electromagnetic field as well as the energy balance. Also we can clearly see the different regimes of the cascade development in both cases.
The first regime is observed when a 0 of the laser pulse is not big enough or the gamma-bunch is not dense enough. In that case the density of the produced electron-positron plasma does not reach the relativistic critical density so that v x ≈ 1, i.e. the collective plasma effects do not occur. In this case the plasma region is not present at all and the newly born particles move in the unaltered field of the laser pulse which is close to the plane wave. As discussed in 7,[43][44][45] in that case χ of the pairs does not grow during its motion in the plane wave. But after each act of the gamma-quant emission it splits between the parent and the child particles so after few generations χ of all the particles becomes negligibly small so the cascading ceases. Thus for small enough a 0 the gamma-quanta of the gamma-bunch decay into pairs levaing the 'trail' of the electrons and positrons which are accelerated forwards and co-propagate with the laser pulse. Although a 0 =2500, ε γ, 0 = 200, ν=0.32, µ=10  In the second regime the cascade develops as discussed in Sec. I. The peak of the pairs density propagates towards the laser with a much slower velocity (relative to the leading edge of the laser pulse) than in the first regime. Moreover the density of the plasma grows in time in contrast to the first regime where the plasma density at each point stays almost the same after the initial gammabunch passes that point. As mentioned in Sec. II E the dense electron-positron plasma actually almost does not absorb the laser field that is why despite the fact that in this regime the total number of the pairs is much larger than in the first regime, the rates of the energy transfer from the EM field to the pairs are close to each other in both regimes.
If the a 0 lays in-between the values of a 0 at which either the first or the second regime is observed, at the initial a 0 =1500, ε γ, 0 = 200, ν=0.35, µ=10  . At some point the density of the pairs becomes large enough to alter the laser propagation and to shift the cascade dynamics to the self-sustained regime. The change between these two regimes is indicated by abrupt change in the velocity of the cascade front. The initial stage (stage of the S-type cascade) can also be seen for larger values of a 0 (see Fig. 6), though it is much shorter and is hardly pronounced in the results of the QED-PIC simulations.
There are some features that are not captured by our model which are worth noting. Firstly, in the PIC simulation there is a distinct tail of the gamma-quanta spatial distribution counter-propagating to the laser pulse. These gamma-quanta have relatively low energy and thus are unable to photoproduce pairs. Our model predicts that the edge of the plasma and gamma-quanta distributions almost competely coincide. The total energy car-ried away by this sort of gamma-quanta is insignificant so this feature is not crucial for the cascade development.
The reason that our model cannot capture this feature is the fact that we assume the distribution functions to be monoenergetic. Higher accuracy can be obtained if we would split the gamma-quanta into several groups with different energies and describe them separately; then this feature would be present in our solutions. But as already mentioned in Sec. II it will greatly complicate the model but will not lead to significant qualitative changes in the solutions. Secondly, both the total number of the pairs and peak plasma density are larger in the QED-PIC simulation than in our model solution for smaller values of a 0 . One of the reasons behind these discrepancies is that our model is 1-dimensional and thus does not describe the laser pulse diffraction. In the 3D-PIC simulations the simulation box is always limited thus the pure plane wave cannot be achieved in the simulations. It leads to the fact that the envelope of the laser pulse evolves so that the magnitude of laser pulse in the vacuum region a 0 =1000, ε γ, 0 = 200, ν=0.35, µ=10

IV. CONCLUSION
We have developed a self-consistent model of the QED cascade development in an extremely intense laser pulse. The complete description of that interaction requires solving the Maxwell's equations along with the kinetic equations for electrons, positrons and gamma-photons. This system of equations is too complex for analytical methods and is usually solved numerically with QED-PIC codes consuming a lot of computational resources.
In order to derive the reduced equations for the computationally light model we adopt some assumptions. The main of them are: quasi-one-dimensional hydrodynamics; locally quasi-monoenergetic distribution function for the particles; the plane wave approximation for the laser radiation. The derived simplified system of the equations is written in the closed form and is solved numerically. Despite the complexity and non-linearity of the cascade dynamics it turned out that a relatively simple one-dimensional model can predict its development qualitatively, e.g. the macroscopic spatio-temporal particles distribution as well as the energy balance in the system in our model and in the results of QED-PIC simulations coincide approximately. These facts serve as the justification of the analytical reasoning behind the model and hence our understanding of the phenomenon. Although there are several discrepancies between the predictions of our model and results of the QED-PIC simulations, the reasons behind them are identified and for some the method to resolve them is discussed.
We stress that the model is suitable for the complete description of the QED cascade development in the single laser pulse. For example various regimes of the interaction based on the intensity of the laser pulse are revealed in both full QED-PIC simulations and our model.
We believe that the model also can be adapted to explore regimes of QED cascades in different environments: laser interaction with the targets including foils, particle beams or gamma-quanta, beam-beam interaction etc. Therefore further extensive test of the developed model is planned to be carried out in future. Let us consider a single electron motion in the field of a plane circularly polarized electromagnetic wave. The vector potential of the wave is chosen as follows A = a 0 (e y cos ξ + e z sin ξ) , The Hamiltonian for the problem is where P = p − A is the generalized momentum of the electron with the momentum p. Deriving the motion equations we get From these Eqs. we derive that Let us define the initial conditions of the electron We can then rewrite Eqs. (A8) -(A10) The Eq. (A3) can be rewritten in the following form The energy gain ∆ε = ε − ε 0 can be analytically found where ξ is the current phase of the electron. If we assume that 1 p 0 a 0 then the second term in the brackets of Eq. (A23) is insignificant for any ξ and ξ 0 compared to the first term and so we can omit it. Another reasoning that allows us not to account the second term is the fact that it depends on the abolute phase in which the particle is located and thus its value averaged over the wave period is equal to zero.
Next we calculate ∆ξ If we assume again that a 0 p 0 then Eq. (A26) can be rewritten as follows Note that ρ 0 < p 0 a 0 and ∆x is the distance along the x-axis between the initial and final position of the electron. In the case of our cascade model the initial position coinsides with the position of the parent gamma-quant decay and the final position is the position of the cascade front. QED-PIC simulations show that the distance ∆x does not exceed several λ or in dimensionless variables it means that ∆x/2π ∼ 1. So in Eq. (A27) we can assume that ∆ξ 1 and leave only the first non-vanishing term in the lhs expansion Substituting this solution to Eq. (A24) we get where we also assumed that ρ 0 ≈ ε 0 (1 − cos θ) The duration of the particle motion can be also calculated and under assumption that ∆ξ ∆x which follows from the assumption p 0 a 0 we get that ∆x ≈ ∆t. It is an obvious conclusion because under our assumptions the particle can be considered ultrarelativistic and its velocity along the x-axis is close to the speed of light.
Below we estimate the radiative losses during the particle motion in the plane wave. The governing parameter χ in the plane wave is caluclated as follows (A32) Where E = −∂A/∂t = a 0 (−e y sin ξ + e z cos ξ) is the electric field and B = ∇ × A = a 0 (−e y cos ξ + e z sin ξ) is the mangetic field. It is seen from Eqs. (A10) and (A32) that in the classical aprroach χ is constant during the particle motion. If we account QED effects then χ parameter changes after each act of gamma-quanta radiation. Although this process can be considered instantaneous, for the sake of estimating the energy loss it is convenient to introduce a continuous force of the radiation friction F rr If we assume the particle to be ultrarelativistic then Eq. (A35) does not explicitly depend on the radiative losses as F rr /ε 2 1. Although due to radiative losses ε changes differently in time as compared to the case without radiative losses, we will assume that the equity v · E = [v × B] x , which follows from the conservation of the γ − p x , still holds. Then the equation for the χ can be significantly simplified where F rr ≡ I rad (χ) and I rad is defined in (10) and expression for which can be found, for example, in Ref. 18 . We assume again that the term v x can be calculated according to the classical approach without accounting radiative losses, i.e.
So finally The total energy loss due to radiation then is calculated as follows The values of ∆ε acc and ∆ε rad obtained using this approach overestimate the corresponding values calculated numerically from the solution of the particle motion equations with the account for the radiation friction [Eqs. (A33) and (A35)], though the order of magnitude is estimated correctly. In the cascade model the value ∆x is unknown because it is the distance along the x-axis between the position of the particle birth and the position where the particle crosses the moving cascade front. This quantity cannot be calculated in our model. Also we are interested only in the total energy change during the time the particle stays in the vacuum region, so the exact dependency of the particle energy on the time is irrelevant. Thus we calculate the values ∆ε acc and ∆ε rad as follows ∆ε acc = µ2 1/3 a where µ is the fitting parameter which determines the characteristic time the particle spends in the vacuum region. We set this parameter equally for all the particles. It turns out that by tuning this parameter both the values of ∆ε acc and ∆ε rad can be estimated with good accuracy (see Fig 9).
Appendix B: Effective dielectric permittivity of the e + e − plasma Let us consider the propagation of a circularly polarized plane electromagnetic wave along the x-axis inside a medium, inhomogeneous along the x-axis, which induces the current j = 2n p v. The Maxwell's equations take the form It is convenient to introduce the following complex variables where v E and v E⊥ are the plasma velocities along and across the electric field respectively. We may also introduce the vector potential a as follows (v E + iv E⊥ ) .
We seek the solution in the form of a plane monochromatic wave with a variable amplitude where both E(x) and κ(x) are real functions and E(x) is the amplitude of the electric field. The final form of the equations is If the plasma is slightly inhomogeneous then we can apply the WKB approximation to solve the problem. If the plasma density distribution has a form of a inhomogeneous slab (as in the case of the cascade development) then this approximation is valid inside the plasma and may be invalid near the edges. Assuming the scale of the plasma inhomogeneity L is larger than the laser wavelength λ we can neglect the term with the second derivative: ∂ 2 E/∂x 2 ∼ E/L 2 k 2 E = (2π) 2 E/λ 2 . So Solving this equation we get where we v 2 ⊥ = v 2 E + v 2 E⊥ = v 2 y + v 2 z . We specify the expression v ⊥ as follows [see Eq. (44)] Note that in the case ν > 0 according to (B15) B > E and thus 1/κ has a meaning of the drift velocity v x , so The comparison of this approximate expression with the exact numerical solution of Eqs.(B12)-(B13) is shown in Fig. 10 for both cases whether the WKB approximation is valid or not.
Since v 2 = 1 is assumed in the derivation of Eq. (B18) it is valid for any reference frame where particles are ultrarelativistic, e.g. the laboratory reference frame.