QED cascades induced by circularly polarized laser fields

The results of Monte-Carlo simulations of electron-positron-photon cascades initiated by slow electrons in circularly polarized fields of ultra-high strength are presented and discussed. Our results confirm previous qualitative estimations [A.M. Fedotov, et al., PRL 105, 080402 (2010)] of the formation of cascades. This sort of cascades has revealed the new property of the restoration of energy and dynamical quantum parameter due to the acceleration of electrons and positrons by the field and may become a dominating feature of laser-matter interactions at ultra-high intensities. Our approach incorporates radiation friction acting on individual electrons and positrons.


I. INTRODUCTION
The dramatic progress in laser technology has enabled a novel area of studies exploring laser-matter interactions at ultra-high intensity [1]. The intensity level of 2 × 10 22 W/cm 2 has recently been achieved [2] and two projects [3,4] aiming at intensity levels up to 10 26 W/cm 2 have been supported and are under way. Furthermore, several original proposals have been suggested e.g. [5][6][7], which reach even higher intensities with almost the present level of technology. One of the key phenomena of laser-matter interactions, that probably dominates at ultra-high intensities of our interest, is the occurrence of QED cascades [4,[8][9][10]. These cascades (also called avalanches, or showers) are caused by successive events of hard photon emissions and electron-positron pair photoproduction by hard photons. As predicted in [10] based on qualitative estimations, the cascades may arise as soon as the laser field strength exceeds the threshold value of E * = αE S , where α = e 2 / c ≈ 1/137 is the fine structure constant and E S = m 2 c 3 /e = 1.32 × 10 16 V/cm is the characteristic QED field. Such a field strength corresponds to an intensity of ∼ 10 25 W/cm 2 .
Previously QED cascades have been observed and studied as a part of Extensive Air Showers (EAS) in the context of the passage of ultra-high energy particles, that originate from Cosmic Rays, through the atmosphere [11][12][13]. However, similar processes can be observed in the external electromagnetic field as well. In this case, Bremsstrahlung is replaced by the non-linear Compton scattering and Bethe-Heitler process is replaced by the non-linear Breit-Wheeler process. The latter processes are well studied both theoretically [14][15][16][17][18][19] and in laser experiments [20] and are probably of great importance for astrophysics (see, e.g., [21]).
An important novel distinctive feature of the cascades in the ultra-strong laser field, compared to situations ever studied previously, is that the laser field is not only able to be a target for ultra-relativistic electrons and hard photons, but can also accelerate the charged particles to ultra-relativistic energies.
As a result, the cascades can be produced even by initially slow electrons or positrons, if they were somehow injected into the strong field region. Moreover, the mean energy of the particles is no longer decreasing in the course of the cascade development due to its redistribution among the permanently growing number of the created particles, but rather is restoring at the expense of the energy extracted from the laser field. This must lead to a vast increase of the cascade yield, as compared to the cascades in media or in strong magnetic fields. In this case the cascade multiplicity would be restricted either by the duration of stay of the particles in the focal region of the laser field, or even, under more extreme conditions, by the total energy stored in the laser field. In the latter case the focused laser pulses would be depleted by cascade production.
As it will be explained in more details below, the restoration mechanism works if the particles can be accelerated transversely to the field. It was conjectured [8,10] on the basis of qualitative analysis for the model of a uniformly rotating electric fields, that this may be indeed the case.
In the EAS theory, the 1D approximation is often used because spreading in transverse direction is inessential for ultrarelativistic particles and has no significance for that problem. Besides, the cascade equations can be solved in this case analytically within the ultra-relativistic approximation by means of the Mellin transform [12,13]. The results of such analytic theory are in good agreement with both experiments [11,12] and direct Monte-Carlo simulations [22]. The attempts to treat the cascades in strong magnetic fields on similar grounds are also known [23]. However, though the 1D approximation remains valid, the cascade equations can not be simplified via the Mellin transform unless some further approximation is made. According to Monte-Carlo simulations [24], the resulting analytical approach works much worse here than in the case of 1D approximation for EAS. In our case of cascades arising in a laser field, the structure of the cascade equations (see Appendix A) is the same as for the mag-netic field, but it is impossible to incorporate restoration mechanism within the 1D approximation in momentum space. This means that our problem is essentially twoor three dimensional.
In this work we report on the first results of the Monte-Carlo simulations of cascades produced by initially slow electrons in a uniformly rotating homogeneous electric field. Such a field can be obtained practically at the antinodes of a standing electromagnetic wave. The choice of the field model is uniquely specified by the existence of reasonable qualitative estimations for scaling of the basic cascade characteristics for this particular case [10]. Our goal was to prove explicitly the existence of the restoration mechanism and to test the estimations [10] by direct numerical simulations.
The paper is organized as follows. In Sec. II, which can be considered as a technical introduction, we review and collect the known information on the elementary quantum processes: single photon emission by electrons and pair creation by hard photons in strong fields of arbitrary configuration. Though this information is not completely new, it is of essential importance for our presentation and is spread among the literature on the subject. After that, in Sec. III we present the reasoning in favor of the energy restoration mechanism for cascades in electromagnetic fields. In Sec. IV we formulate the assumptions of our model, present the details of our Monte-Carlo routine and discuss the results obtained by numerical simulations. These results are compared to the known estimations. Summary and discussion is given in Sec. V. Finally, in Appendix A, we discuss the cascade equations for our problem and explicitly demonstrate that, contrary to the recent doubts [25], the approach we use takes proper account of radiation friction by ultrarelativistic electrons.

II. QUANTUM PROCESSES WITH HIGH-ENERGY PARTICLES IN A STRONG ELECTROMAGNETIC FIELD
General properties of radiation of ultrarelativistic particles are well known [26]. Due to the relativistic aberration effect the momenta of the products of any decay of an ultrarelativistic particle are directed within the narrow angle ∆θ ∼ γ −1 with its momentum, where γ = 1 + (p/mc) 2 . Thus, radiation of a charged ultrarelativistic particle is visible at a point of observation only for a short period of time τ during which its momentum turns through the angle of the order ∆θ. The momentum turning angle can be estimated by ∆θ ∼ eF ⊥ τ /mcγ, where F ⊥ denotes the characteristic value of the transverse component of the field. Thus, τ ∼ mc/eF ⊥ . The characteristic frequency Ω of classical radiation can be most simply estimated in the proper reference frame of the particle, by transition to which the duration τ transforms into τ ′ = τ /γ. Thus, Ω ′ ∼ τ ′ −1 ∼ eF ⊥ γ/mc. In the laboratory frame, due to the Doppler up-shift, we would thus have Ω ∼ (eF ⊥ /mc)γ 2 . Such a scaling with γ is typical for congeneric problems and arises e.g. in the theory of synchrotron radiation.
The parameter χ = Ω/(γmc 2 ) = F ⊥ γ/E S [27], being the ratio of the classically estimated mean energy of an emitted photon to the energy of the radiating particle, determines whether the process of radiation is controlled by classical electrodynamics or QED. Namely, if χ ≪ 1, then quantum recoil is inessential and radiation is classical, whereas if χ 1 then it must be quantum. The parameter χ is Lorentz and gauge invariant and is precisely defined as χ = e /(m 3 c 4 ) −(F µν p ν ) 2 , where F µν is the strength tensor of the field and p µ is the 4-momentum of a particle. In what follows, we assume that F ≪ E S . On the other hand, we assume that the field is of relativistic strength in the sense that the dimensionless field amplitude a 0 = e −A µ A µ /(mc) ≫ 1, where A µ is the 4-vector of the field potential. The latter means, in particular, that it varies on the scale that exceeds τ and thus can be considered constant with respect to the decay processes.
Two different theoretical approaches have been developed in order to study photon emission by ultrarelativistic (γ ≫ 1) charged particles in electromagnetic fields of ultrarelativistic (a 0 ≫ 1, χ 1) but still subcritical (F ≪ E S ) intensities. Nikishov and Ritus (NR) have calculated the appropriate quantum amplitudes in terms of Volkov solutions in a constant crossed field (E 2 − H 2 = 0 and E · H = 0) of arbitrary strength [15,27]. Their results can be applied directly to our problem because, as they have pointed out especially, under the abovementioned conditions any field looks locally as constant and crossed, the latter in the sense that both field invariants (E 2 − H 2 )/E 2 S and E · H/E 2 S are much less than χ 2 . The other approach by Baier and Katkov (BK) [19,27] is based on the observation that the motion of a particle in between two acts of photon emission (which corresponds to free lines of Feynman diagrams) can be considered classically if F ≪ E S , so that all the relevant quantum corrections are reduced to quantum recoil and, possibly, to field-spin interactions. In a sense, the BK approach is equivalent to the replacement of the aforementioned Volkov solutions by the localized wave packets moving along the classically prescribed trajectories. Both approaches are based essentially on the same approximations and provide the same results for the energy spectra of emitted photons.
The energy distribution of the probability rate for photon emission by ultrarelativistic electrons in an electromagnetic field is given by [15,19,27] dW rad (ε γ ) where x = (χ γ /χ e χ ′ e ) 2/3 , Ai(x) = (1/π) ∞ 0 cos(ξ 3 /3 + ξx) dξ is the Airy function, ε γ and ε e are the energies of the emitted photon and the initial electron, respectively. χ e , χ ′ e = χ e − χ γ and χ γ (0 < χ γ < χ e ) are the dimensionless quantum parameters for the electron before and after emission, and for the emitted photon, respectively. In terms of the field strengths, this parameter is represented as where ε and p are the energy and the momentum of the corresponding particle. Note that expression (1) suffers from the infrared singularity at ε γ → 0. However, in our constant field approximation this singularity dW is weaker than the usual O(ε −1 γ ) scaling of the infrared behavior of perturbative QED [27,29]. In particular, the total radiation probability rate is infrared convergent in our approximation. The infrared sector, however, is not important for the parameters considered in our paper, because most of the emitted radiation are found to have much larger frequencies than the frequency of the driving field.
The energy distribution of the probability rate for direct pair creation by hard photons (ε γ ≫ mc 2 ) is given by [15,19,27] where the indices "γ" and "e" this time refer to the initial photon and to the created electron, respectively. For the created positron, we have χ ′ e = χ γ − χ e (0 < χ γ < χ e ). Formula (3) is completely symmetric with respect to electron and positron remaining unchanged under the replacement χ e ↔ χ ′ e . Similarity between formulas (1) and (3) is explained by the fact that these two processes are related by the cross-symmetry [27].
The total probability rates for both processes cannot be written in terms of known special functions and should be obtained by numerical integrations. However, they allow simple asymptotic expressions in the limits of small and large χ e , χ γ , respectively. Namely, we have and Eqs. (5a), (6a) in these formulas describe the quasiclassical regime. For small values of the quantum parameter χ γ , the probability rate for pair photoproduction W cr is suppressed exponentially, in accordance with the essentially quantum nature of this process. At the same time, W rad remains O( −1 ), thus providing a finite classical limit for the mean radiated intensity I rad . As for the limit of large χ, both rates (5b) and (6b) differ only by a numerical factor of the order of unity. Given the energy and the momentum of the electron before emission, Eq. (1) determines the probability distribution for the energy ε γ of the emitted photon. Under our assumption γ e ≫ 1, the momentum of this photon is given by p γ = (ε γ /p e )p e . The energy and the momentum of the electron after emission should be determined from the conservation laws. In the electromagnetic background, they are of the form p µ e + q µ = p ′ e µ + p µ γ , where q µ is the four-momentum extracted from the field. The exact value of this q µ essentially depends on the global structure of the field. This is because the whole spacetime contributes to the integrals in the QED matrix element that yield delta-functions expressing the conservation laws. For example, in a crossed constant field with the Poynting vector directed along the z-axis the conserved quantities are ε/c − p z , p x and p y . In a constant electric field, the canonical momentum is conserved. In a constant magnetic field, directed along the z axis and with symmetric gauge, the conserved quantities are the number of the Landau level, the angular momentum, and p z . Nevertheless, for ultrarelativistic particles there is actually no difference between these possibilities, and either of them can be adopted with the accuracy of our approximation. The reason is that q eF τ mc ≪ p e , p ′ e , p γ . In particular, we can assume p ′ e = p e − p γ . The same argument can be applied to pair creation by hard photons as well.
In addition to one-photon emission and direct pair photoproduction reviewed above, there exist more complicated higher-order processes, such as e.g. the twophoton emission e − → e − γγ or the trident process e − → e − e − e + . Their specific feature is that the intermediate particle is off the mass shell, i.e. is virtual. However, in strong field situations of our interest the two-step processes dominate [8,9,20]. For this reason, we do not consider higher-order processes in the sequel.

III. BASIC ESTIMATIONS FOR CASCADE PRODUCTION IN A ROTATING ELECTRIC FIELD
Since there is no difference whether an electron or a positron initiates a cascade, we assume in this section that our cascade is initiated by a positron (e > 0).
Consider a positron in a homogeneous, uniformly rotating electric field The equation of motioṅ with the initial condition p(t 0 ) = p 0 , can be easily solved: Here, a 0 = eE 0 /mωc is the dimensionless field amplitude. Let us assume first that the positron is at rest (p 0 = 0) initially (t 0 = 0). Equations (7) and (9) show that the energy and the quantum parameter χ of the positron for this case depend on time as Both quantities are increasing initially. They are oscillating with the period 2π/ω of the rotation of the field. The amplitudes of these oscillations, ε m ≈ 2mc 2 a 0 , are proportional to a 0 and a 2 0 , respectively, and are quite large under our basic assumptions. For example, χ m approaches unity already at a 0 ∼ a 0c = 500 for an optical rotation frequency of ω = 1eV. This corresponds to the field strength E 0 ∼ 10 −3 E S ∼ 10 13 V/cm and the intensity 10 24 W/cm 2 . Since χ m ∼ 1 under such conditions, the positron, according to the preceding section, is able to emit a hard photon with χ γ ∼ χ e ∼ 1, which, in turn, can create an electron-positron pair. However, at such intensities a new generation of pairs is typically produced on the time scale π/ω, and the whole pair generation process may be rather sensitive to peculiarities of the field model. As we discuss below, stable cascade formation is expected at higher intensity levels.
The formulas (10a) and (10b) become especially simple for stronger fields a 0 ≫ a 0c , because in this case the value χ e ∼ 1 is being reached within just a small fraction t acc of the rotation period. Namely, we have Eq. (11a) is easy to understand, because initially the positron is accelerating almost along the field. In order to understand Eq. (11b) better, let us note that in the case p(0) = 0, according to Eqs. (7) and (9), the momentum of the positron constitutes exactly the angle ωt/2 with both the instant direction of the field and the x-axis. Intuitively, this is because due to its inertia the particle does not follow the rotation of the field precisely. As a consequence, the transverse component of the field with respect to momentum of the particle increases as we arrive immediately at Eq. (11b). Qualitatively, the same growth of the energy and the parameter χ with time has been observed for generic field configurations [10].
As it follows from Eq. (11b), the quantum parameter χ e becomes of the order of unity over the period of time t acc , Here we have introduced a new dimensionless field inten- , which is appropriate for the cascade problem [10]. The parameter µ is related to the commonly accepted parameter a 0 by µ = ( ω/αmc 2 )a 0 . According to Ref. [10], the cascades can be caused by initially slow particles if µ 1.
In the course of hard photon emission, the value of the quantum parameter χ e is shared between the positron and the emitted photon [35], χ e ≈ χ γ + χ ′ e . If χ e 1 then both χ γ and χ ′ e are less than χ e but are of comparable value χ ′ e ∼ χ γ χ e . Although propagation of the resulting hard photon is not affected by the field, nevertheless its χ γ continues to increase after emission just due to rotation of the field.
In order to understand better what must happen after the first hard photon emission, let us come back to Eq. (9) and consider the general initial condition. It is easy to see that in the case H = 0 the sign of the derivative of the quantity (2) is determined completely by the expression −e(p · E)(p ·Ė). The zones in the p 0 -plane in Fig. 1, whereχ e > 0 is valid for positrons at time t = t 0 are shaded. In the shaded areas χ e can be expected to increase for some time. The time t 0 can be identified with the creation time of a new pair.
Since the momentum of a primary positron is confined to the shaded region to the right, and the new secondary particles are created with momenta parallel to the momentum of parental particle in our approximation, we see that momenta of the secondary particles also lie in the shaded sector. Thus, the newly created positrons are accelerating withχ(t) > 0, while the newly created electrons are initially decelerating (χ < 0). However, they are quickly turned back by the field and also get accelerated. The same will be true for successive pair creation processes as well.
In spite of the complexity of the picture of the cascade development (see Fig. 2), some general estimations for it can nevertheless be obtained [10] in the high-field limit µ ≫ 1. The idea is that, due to similarity of the rates (5b) and (6b), and also because the variation of the angles between the momenta of all particles (positrons, electrons, and photons) and the field is determined by the same temporal scale ω −1 , there is actually no need to distinguish between all three sorts of particles. So, an order of magnitude estimation can be provided within the model of a simple doubling chain process.
Let us denote by t e the typical lifetime for electrons and positrons with respect to hard photon emission. The same quantity up to an order of magnitude defines the lifetime of photons with respect to pair creation. The lifetime t e , together with the typical energy and the value of quantum parameters of the particles, as well as with the angle between their momenta and the field, can be estimated as [10] t e ∼ αmc 2 µ 1/4 Under the condition µ ≫ 1, as is assumed here, we have θ ≪ 1 and χ ≫ 1. The latter inequality approves the choice of the asymptotic expressions (5b) and (6b). In addition, we have the following hierarchy of the time scales t acc ≪ t e , which assures that exactly hard photons with χ γ 1 are typically emitted. Within the framework of the doubling chain process model, the number of pairs (multiplicity of the cascade) must grow exponentially, In the next section, we are checking the estimations (13) and (14) by direct Monte-Carlo simulations.

IV. DESCRIPTION OF MONTE-CARLO APPROACH AND NUMERICAL RESULTS
In our simulations we are using a Monte-Carlo approach for the integration of the cascade equations [see the Eqs. (A1), (A2)]. We trace the motion of the electrons and positrons in between the photon emissions classically, whereas for hard photons we exploit the ray tracing approximation in between their emission and conversion into pairs. Even though there exists the exact analytical solution (9) for equations of motion (8) for positrons and electrons, we are integrating Eq. (8) numerically for each of the particles. This is done in order to incorporate the probabilistic events of photon emission and pair creation in the routines as described below, as well as for the purpose of future generalization to more realistic field configurations.
Our numerical algorithm works as follows. At each time step t i < t < t i +∆t we are calculating the momenta of all the particles created at preceding time steps by p i+1 = p i + q i E i+1/2 ∆t, where E i+1/2 = E(t i + ∆t/2) and q i = +e, −e, 0 for positrons, electrons and photons, respectively. The event generator determines which of the electrons or positrons is going to emit a photon at this time step and whether any of the present photons is going to produce a pair.
Let us explain our event generator for photon emission in more details (see also Refs. [22,28]). Starting from p i and E i = E(t i ), we attach the value χ i at time t i using Eq. (2) to each electron and positron and compute the total probability rate W rad (see Eq. (4)). In order to isolate the infrared singularity, we set the lower limit of integration to ε min . For each electron and positron, we assume that it emits a photon between t i and t i+1 if r < W rad ∆t, where r (0 < r < 1) is a uniformly distributed random number. If the above inequality is fulfilled, then the energy ε γ of the emitted photon is obtained as the root of the sampling equation where r ′ (0 < r ′ < 1) is an independent random number. The time step ∆t, which remains fixed in the course of computation, must be chosen such that ∆t ≪ W −1 rad , W −1 cr holds. The direction of propagation of the newly emitted photon is parallel to the momentum p i of the parental electron or positron, whose momentum after emission we find from the conservation law as discussed in Sec. II. For pair creation, the event generator works similarly, apart from the fact that there is no need for the regularization parameter ε min . Within the constant crossed field approximation applied here we assume that ε γ ≫ mc 2 . However, the photons with energies ε γ mc 2 are not able to create pairs in a subcritical field, for which χ γ ≪ 1 holds. Therefore, we can completely neglect emission of soft photons in our problem. Based on this reasoning, we currently set the lower integration limit ε min to mc 2 for both W rad and the sampling Eq. (15).
As a benchmark for our code we have simulated the development of a cascade initiated by a high-energy (ε 0 = 2 × 10 5 mc 2 ) initial electron in a constant homogeneous transverse field with E 0 = 0.2E S . Our results are averaged over 10 3 simulation runs. In this particular simulation, the curvature of trajectories of electrons and positrons has been neglected, so that the results of our simulations can be directly compared with previous simulations of cascades produced by high-energy electrons in a magnetic field [24]. Comparison of cascade profiles obtained in both simulations is given in Fig. 3 by the solid red and dashed green lines, respectively. The figure represents the number of pairs with an energy exceeding 0.1% of the energy of the primary electron versus the elapsed time. In our notation the reference characteristic radiation time t rad as adopted in Ref. [24] in our notation is t rad = 3.85 × (γ in /αχ 2/3 in ) × ( /mc 2 ) = 5.64 × W −1 rad,in , where the subscript "in" refers to the initial data for primary electron. We see that our results are in reasonable agreement with the paper [24].
We have also implemented and tested a different event generator, which provides significant speed up due to the absence of numerical integrations. The idea is to exploit some explicit algebraic fits for the energy spectrum (1), and to exchange the order of testing the occurrence of photon emission and of sampling its energy. In this alternative version of the algorithm, within each time step one first samples the possible energy of an emitted photon just as an uniformly distributed random quantity, ε γ = ε e r ′ in the above notation. After that, photon emission is assumed to take place if r < [dW rad (ε γ )/dε γ ]ε e ∆t. In this case, the time step must satisfy the condition ∆t < [ε e dW rad (ε * γ )/dε γ ] −1 for all appearing electrons and positrons, where ε * γ is the photon energy that corresponds to the maximum of the emission spectrum (1). The same scheme can be applied to the simulation of pair photoproduction as well. Note that in this case there is no need to introduce the energy cutoff ε min , although this may serve as a useful trick if one wants to restrict the number of soft photons that are traced by the code. The test of the modified event generator is included by a dashed blue line in Fig. 3. This test demonstrates that both versions of the event generator are in fact equivalent.
The results of our simulations are collected in the figures 2, 4-7. Fig. 2 represents a typical spatial picture of the formation and development of a cascade initiated by a positron. The electrons and positrons are deflected by the field in opposite directions, whereas the directions of propagation of photons are distributed randomly, as could be expected. For the rest of the paper we assume for all simulations in an uniformly rotating field that at t = 0 we have a single electron at rest (p e = 0) and no photons and positrons. The typical evolution of the quantum dynamical parameter χ e of the primary electron is illustrated with the left panel of Fig. 4. Before the emission of a first photon, the electron is gaining energy and its parameter χ e is growing as the square of time in accordance with Eq. (11b). After the first photon emission, which for our parameters happens typically on the time scale t e smaller than ω −1 , the curves become stochastic and consist of smooth sections with typical growth of χ e due to acceleration by the field. These sections are separated by sudden breakdowns resulting from recoils due to successive photon emissions. Since these recoils are random, the three curves in the figure corresponding to independent simulation runs deviate at later times. After the transient period which typically lasts for several lifetimes t e , the momentum losses due to quantum recoils are coming into equilibrium on the average with the trend of acceleration by the field. After that, function χ e (t) for an individual electron describes a stationary stochastic process.
As it was predicted in Ref. [10], the development of a cascade results in exponential growth with time of the total numbers of secondary hard photons and electronpositron pairs. This is illustrated with the right panel in Fig. 4. The plot N e − e + (t) is given by a random stairway, with each stair corresponding to creation of a single pair.
The successive stairs are well separated initially, when the total number of pairs remains small. At later time with the number of pairs growing rapidly the stair-like structure of the lines in the plot becomes invisible and straight lines are obtained. Although these straight lines for independent simulation runs are typically different, mostly because emission of the first photon starts randomly from one simulation run to another, nevertheless their gradients are varying weakly in different runs and can be used to determine the growth rate Γ in Eq. (14). For example, the growth rates extracted from the curves 1-3 at Fig. 4 are Γ = 4.62, 4.84 and 4.90, respectively.
We have studied the averages of the quantities χ e , ε e and θ over the cascade. For example, temporal evolution of the mean value where N e − (t) is the instant number of present electrons and χ e i (t) is the instant value of the quantum dynamical parameter for the i-th electron, as depicted in the left plot of Fig. 5. One can see that at later times the random fluctuations are smoothed out and the quantity (16) stabilizes acquiring some definite constant value which is independent of the simulation run. The same behavior was observed for ε e and θ , which are defined in a manner similar to Eq. (16). The typical evolution of the averaged energy of all the components of the cascade is represented in the right plot of Fig. 5. At later times, the mean energies of electrons and positrons coincide as is expected from symmetry consideration, whereas the mean photon energy typically remains smaller. At the same time, the energy spectrum of created electrons and positrons is wider than the photon spectrum (see Fig. 6). Both features are explained naturally by the fact that in our setup the hard photons (χ γ 1) are quickly converted into pairs which survive, whereas soft photons (χ γ 1) are stable with respect to pair photoproduction and hence are accumulated. In the high energy region, all the spectra are likely to show exponential decay.
One of our main tasks was the investigation of the validity of estimations (13) and (14) which were suggested previously in Ref. [10] and are of crucial importance. In particular, Eq. (14) was serving as an argument for the estimation of the maximal value of the intensity attainable with focused laser fields. In order to test these estimations we have performed parametric studies of the stabilized values of the quantities χ e , ε e , θ and of the increment Γ. These results are presented in Fig. 7. The left plot demonstrates the dependence of the ratios of the quantities obtained by simulations to their estimations (13) on µ for fixed rotation frequency ω = 1eV. It is clear from the figure that for large values of µ each ratio acquires some definite limit of the order of unity. According to our results the formulas (13) are valid up to some numerical coefficients of the order of unity, which vary no more than twice in the whole range µ > 1. The results of simulations for Γ are compared with Eq. (14) for two different rotation frequencies ω = 0.66eV and ω = 1eV on the right panel of Fig. 7. One can see that for large values of µ the estimation (14) is justified with good accuracy even without any correction factor. For µ 30 formula (14) overestimates Γ but not more than by half of an order of magnitude. This may be nevertheless crucial for the estimation of the total cascade yield due to its exponential dependence on Γ. For the particular value µ ≈ 10, which was exploited in Ref. [10], formula (14) overestimates Γ by approximately a factor of 1.5. This, however, can be compensated in principle by simultaneous underestimation of the escape time in Ref. [10].
In order to apply the results of our simulations to estimate the cascade yield by a realistic focused laser field, we assume that the appearing electrons and positrons are pushed away as a whole from the focus by the ponderomotive potential in radial direction with almost the speed of light. Assuming the Gaussian profile of the focused beam we can write µ(t) = µ 0 e −c 2 t 2 /w 2 0 , where µ 0 is the value of the parameter µ at the center of the focus and w 0 is the focal radius. Then the total number of pairs produced by the cascade can be estimated to the The remaining integral defines the effective time of escape from the focus and equals √ π(w 0 /c), i.e. is √ π ≈ 1.77 times larger than it was assumed in Ref. [10]. This correction almost totally cancels the overestimation of Γ by formula (14) at µ ≈ 10 that we have observed in our simulations. Thus, we hope that the quantitative predictions in Ref. [10] must remain unaffected by our corrections.
One can see that we are currently neglecting the complicating details in our code such as, e.g. the elastic collisions, the Compton scattering, and the annihilation processes. Such phenomena must become important only at later times, when the plasma is dense enough. Though successive collisions and annihilations of the electron and positron from the same created pair may be important [30], we are currently ignoring these effects for simplicity [36]. We ignore the higher-order processes, such as twophoton creation and the trident processes as well (see the remark at the end of Sec. II). All these assumptions are natural and commonly accepted in the present cascade theory [23], even though they may be revised in future studies.
Due to limitations of computer power, we currently stop our simulations after the creation of around N e − e + 10 4 pairs. This was shown to be enough to estimate the growth rate Γ, as well as to average the characteristics of a cascade over the ensemble of pairs with reasonable accuracy. In the time interval of simulation these pairs occupy a volume of the order d 3 , where d ∼ c/ω ∼ 1µm. This corresponds to the pair density n e + e − ∼ 10 16 cm −3 . Typical values of the γ-factor for electrons and positrons are γ ∼ 10 3 ÷10 4 , see Fig. 6. This corresponds to their energies ε e = γmc 2 ∼ 0.5 ÷ 5GeV. Assuming the temperature T ∼ ε e /k ∼ 10 13 K, we can estimate the Debye screening radius r D ∼ kT /e 2 n e + e − ∼ 1cm ≫ d. The relativistic plasma frequency Ω pe ∼ c/r D ∼ 10 10 sec −1 remains about five orders of magnitude smaller than the optical frequency. For these reasons we have completely neglected Coulomb interaction between electrons and positrons and all the accompanying collective plasma effects in the present simulations. However, the density of pairs is growing exponentially, n e + e − (t) ∝ e Γt , and hence r D ∝ e −Γt/2 and Ω pe ∝ e Γt/2 . After a relatively short period of time 2π/ω, when the number of pairs becomes macroscopic (∼ 10 11 ), the quantities r D and Ω pe attain the values d and ω, respectively, and the collective plasma effects may come into play. Within our approximation of a homogeneous field, the total number of created pairs would be restricted by the screening of the external field by the self-field of arising plasma.
Let us note, that despite some doubts in the literature [25], the radiation friction is taken into account properly in our version of the algorithm by the recoils happening at the times of photon emission (see, e.g., Ref. [31] and the Appendix A in our paper). Hence, there is no need to include an additional radiation friction force in the equations of motion (8) for electrons and positrons, otherwise this would cause double counting. Moreover, our approach transfers the concept of classical radiation friction into the quantum domain in a correct fashion. It can be asked how the classical continuously acting radiation friction can be recovered from the sudden jumps of momentum similar to those in Fig. 4. In fact this happens on the average with respect to the ensemble of Monte-Carlo realizations, since the moments of successive photon emissions are distributed randomly. At later times, when the number of created pairs becomes large, the cascade forms a representative ensemble itself and there is in principle no need for taking the average over independent realizations.

V. SUMMARY AND DISCUSSION
In this paper we have presented the first results of numerical simulations of the formation and development of electron-positron-photon cascades by initially slow electrons in a uniformly rotating homogeneous electric field. In such a situation the cascades reveal a principally new feature, i.e. the restoration of the energy and the dynamical quantum parameter due to the acceleration of electrons and positrons by the field. This feature may be of crucial importance for the whole physics of laser-matter interactions in the strong field domain, as it was demonstrated in Ref. [10]. We have explicitly identified this restoration mechanism in the course of our simulations. Also, our simulations clearly confirm the qualitative analysis of Ref. [10], including the basic scalings (13) and the estimation (14) for the cascade yield. So, they can be used to fix the remaining numerical correction factors in Eqs. (13) and (14), which turn out to be of the order of unity.
The numerical approach that we adopt is based on Monte-Carlo simulations of the cascade equations. We have shown explicitly that contrary to some recent doubts in literature [25] such an approach incorporates radiation friction acting on individual electrons and positrons and, moreover, is doing this in a manner which is consistent with intense field QED.
The codes designed for our task can be readily adopted for simulating cascades in the laser fields with more realistic configurations, such as tightly focused Gaussian beams and pulses. This is required in order to make more definite predictions on the impact of cascade production for possible future experiments, as well as for further corrections of the maximal value of intensity that can be attained with optical lasers [10]. Simulation of cascades in a focused laser field will be presented in a separate publication. However, let us make several brief comments about cascades in focused laser fields, possible experimental scenarios, and some yet unresolved technical problems that may require further studies.
The restoration mechanism arises due to the curvature of the trajectories of the charged particles across the field and may be sensitive to its polarization. Although we expect that this mechanism must work for generic field configurations (e.g., for generic tightly focused laser fields), there may exist several particular configurations for which the restoration mechanism does not work. For example, in an arbitrary constant electromagnetic field or a circularly polarized propagating plane electromagnetic wave the dynamical quantum parameter χ e is conserved exactly in the course of motion. In the case of a generic propagating plane wave the amplitude of oscillations of the parameter χ e for an initially slow electron does not exceed E 0 /E S , i.e., always remains smaller than unity. Another example is a linearly polarized oscillating electric field [25], since in this case the initially slow particles are accelerated strictly along the field and hence the growth of the transverse component of the field is absent. In some intermediate cases, e.g. for elliptical polarization or a weakly focused Gaussian beam, restoration of χ e must exist but may be less effective than in the case of circular polarization. However, in all the cases at least the usual cascades would be caused by external high-energy electrons or hard photons passing through the high-field region transverse to the field. In this case, the cascade yield remains microscopic and would be determined by both the initial energy of an external energetic particle and the laser field strength.
In order to initiate a cascade in a tightly focused laser field, it is required to inject a primary particle into the center of a focal region. This task may be not trivial because the focal region is surrounded by a ponderomotive potential wall of the characteristic height U 0 ∼ mc 2 1 + a 2 0 ≈ mc 2 a 0 . The external high-energy electrons will be most likely deflected rather than penetrate inside. The most direct and elegant scenario is based on the exploration of pairs that are created spontaneously from vacuum by the laser field itself [10], since they are appearing exactly at the center of the focus as required. However, this possibility implies high intensities 10 26 W/cm 2 . Another possible resolution would be the initiation of cascades by energetic γ-quanta. In our opinion, the final conclusion whether or not cascades with macroscopic yield can arise in generic real experiments exploring laser-matter interaction at intensities lower than ∼ 10 26 W/cm 2 requires further studies. We note that for cascades that arise in the course of the interaction of high-intensity laser radiation with material targets it may be necessary to take the impact of ordinary cascades in matter into account as well [33].
If the cascade yield attains macroscopic values (N e − e+ ∼ 10 11 ), the self-field of the electron-positron plasma becomes comparable to the guiding field. In this regime screening of the external field and its absorption by the electron-positron plasma self-field will restrict further pair production. Such a regime can be simulated by combining our codes with the Particle-in-Cell (PIC) method [34]. We hope to address this problem in one of our next publications.
differ from the standard equations of EAS [12,13] only by the addition of the second term to the LHS of Eq. (A1), which takes electron and positron acceleration into account. Here, f ± and f γ are the distribution functions for positrons, electrons and photons, respectively. In our approximation of photon and pair emission in strictly forward direction the differential probability rates are of the form so that the integrals standing on the RHS of Eqs. (A1), (A2) are essentially one-fold. Nevertheless, in the main case of interest, the direction of the field E(t) varies in time, so that the problem does not reduce to 1D. Note that the scalings (13) can be obtained from (A1) and (A2) in the limit of large χ via dimensional analysis. It is interesting to note that in the view of the approximation (A3) it follows from the Eqs. (A1), (A2) that These identities ensure that in our approximation the momentum and the energy are extracted from the field only during the acceleration of electrons and positrons, i.e., both energy and momentum of the electron-positronphoton plasma are conserved during photon emission and photons to pairs conversion. The first two terms on the RHS of the Eq. (A1) describe the influence of photon emission on electrons and positrons. Let us now demonstrate that classical radiation reaction is taken into account properly by these terms.
For this aim let us assume in what follows that the electrons are slow in the sense that distributions f ± are restricted to such momenta for which χ e± ≪ 1 (though we continue to assume that they are ultrarelativistic). In this case their motion can be described in completely classical terms. Accordingly, let us skip the third term on the RHS of the Eq. (A1), which is responsible for pair production. Once this is done the total numbers of positrons and electrons N ± = f ± d 3 p e are conserved.
The relation between the variables χ γ and x in Eq. (1) can be expressed in the alternative form It is clear from this formula by taking into account that the spectrum (1) of emitted photons is effectively concentrated in the range x 1, that As a consequence, p γ ≪ p e in the remaining integrals on the RHS of Eq. (A1). Thus, the expansion w rad (p e + p γ → p γ )f ± (p e + p γ ) − w rad (p e → p γ )f ± (p e ) ≈ p γ · ∂ ∂p e [w rad (p e → p γ )f ± (p e )], is valid. Consider the average momenta P ± (t) = (1/N ± ) p e f ± (p e , t) d 3 p e . By multiplying Eq. (A1), truncated as described above, by p e and integrating it over p e by parts we havė where R(p e ) = − p γ w rad (p e → p γ ) d 3 p γ (A10) is the mean rate of momentum losses of electrons and positrons due to photon emission, and R ± (t) = 1 N ± R(p e )f ± (p e , t) are its averages over the momentum distributions of positrons and electrons, respectively.
From this point, let us pass to the integration over the variable x. In the approximation χ e ≪ 1 we have χ γ = x 3/2 χ 2 e . Contribution to the integral in Eq. (A11) comes from the range x ∼ 1. We can neglect the term χ γ √ x in the brackets in the expression (1) and, in addition, replace the upper limit of integration over x by infinity. After these manipulations, we have This is exactly the main contribution to the Landau-Lifshitz (LL) force [26] for ultrarelativistic electrons.
Other contributions to LL equation do not appear in our derivation only because we are essentially using the constant crossed field approximation.
In the quantum case χ e ≫ 1 the expansion (A8) is not valid. Thus, radiation friction in the quantum regime cannot be described by the concept of classical force in principle, as it was attempted to do e.g. in [9]. In addition to the advection term in the transport equation, which could be ascribed to the radiation reaction force as above, spreading of the distribution functions in momentum space becomes important as well. This spreading is associated with quantum fluctuations and is observable, e.g., as quantum excitation of synchrotron and betatron oscillations [32].