High-density electron–ion bunch formation and multi-GeV positron production via radiative trapping in extreme-intensity laser–plasma interactions

Multi-petawatt laser systems will open up a novel interaction regime mixing collective plasma and quantum electrodynamic processes, giving rise to prolific generation of gamma-ray photons and electron–positron pairs. Here, using particle-in-cell simulations, we investigate the physics of the interaction of a 1024 W cm−2 intensity, 30 fs duration, circularly polarized laser pulse with a long deuterium plasma at classically overcritical electron density (1022 cm−3). We show that radiative trapping of the plasma electrons causes a high-density (∼5 × 1023 cm−3), quasineutral electron–ion bunch to form inside the laser pulse. This phenomenon is accompanied by up to ∼40% energy conversion efficiency of the laser into gamma rays. Moreover, we find that both the radiation-modified Laplace force and the longitudinal electric field exerted on the positrons created by the multiphoton Breit–Wheeler process can accelerate them to GeV-range energies. We develop a theoretical model, the predictions of which provide a good match to the simulation results. Finally, we address the influence of the ion mass, showing that the laser absorption and positron acceleration is enhanced with deuterons compared to protons.

The radiation backreaction on the laser-driven electrons-recently evidenced experimentally in the collision of a laser-wakefield accelerated electron beam (of energy E e 500 MeV) with an intense laser pulse (I L 10 20 W cm −2 ) [21,22]-can dramatically alter their dynamics [23][24][25][26][27], and thus modify the characteristic phenomena of (moderately) relativistic laser-plasma interactions, e.g., hole boring [28][29][30][31] or filamentation instabilities [32]. Moreover, it can trigger novel mechanisms such as radiation-mediated, inverse Faraday B-field generation [33], or 'anomalous' radiation trapping, whereby the laser-driven ultrarelativistic electrons tend to be confined around the extrema of the laser electric field [34][35][36]. Despite extensive investigations in past years, many aspects of QED-modified, ultraintense laser-plasma interactions remain poorly understood, particularly the interplay between the highly radiative energized electrons and the ions, both species then sharing similar relativistic inertia.
In this article, through advanced particle-in-cell (PIC) simulations and modelling, we demonstrate that laser-plasma interactions at extreme intensity (10 24 W cm −2 ) and classically overcritical density (10 22 cm −3 ) can generate, via radiative trapping, highly compressed (∼ 5 × 10 23 cm −3 ) electron-ion bunches inside the laser pulse, enabling copious production of gamma-ray photons and electron-positron pairs. We further reveal that the positrons that result from the Breit-Wheeler process can be efficiently accelerated in the forward direction by the radiation-modified Laplace force, together with the induced charge-separation field. Notably, we find that ∼ 10% of the positrons are forward-accelerated up to multi-GeV energies. We propose a detailed analytical description of those effects, which is in good agreement with the simulation results. Finally, we show that both the radiative trapping efficiency and positron yield increase with decreasing charge-to-mass ratio of the plasma ions.
The remainder of the article is organized as follows. In section 2, we detail the parameters of the PIC simulations supporting our study. Section 3.1 describes the typical sequence of events revealed by our simulations, leading to the formation of a very dense electron-ion bunch within the laser pulse, along with the depletion of the latter in the course of its propagation. In section 3.2, we construct an analytical model that captures well the simulation results. In section 4, the properties of the beam of positrons created during the interaction are examined, and we show that a substantial number of them are forward-accelerated up to multi-GeV energies with low divergence, an observation explained in the light of our analytical model. Section 5 addresses the influence of the ion mass on the previously discussed phenomena. Our conclusions are summarized in section 6.

Simulation setup
We have performed 2D3V (two-dimensional in configuration space, three-dimensional in momentum space) PIC simulations using the epoch code [37] describing the interaction of a circularly polarized pulse with a semi-infinite deuterium plasma. The laser intensity profile is Gaussian in time and space with 30 fs full width at half maximum (FWHM) duration and 3.75 μm FWHM width. The maximum laser intensity is 10 24 W cm −2 , corresponding to a normalized potential vector a L ≡ eE L /m e cω L = 850 (E L and ω L are the laser electric field strength and angular frequency, respectively). The initial electron plasma density is set to 10n c , where n c = 0 m e ω 2 L /e 2 is the critical density at the laser wavelength λ L ≡ 2πc/ω L = 1 μm and 0 is the vacuum permittivity. While such a dense plasma would be opaque to mildly relativistic laser pulses, it turns relativistically transparent at the extreme laser intensity considered here. The laser pulse starts irradiating the plasma at t = 0 and the rise time of the pulse is 60 fs. In order to assess the influence of the charge-to-mass ratio of the plasma ions on the laser-plasma interaction, we have also run a simulation with a hydrogen plasma (see section 5).
As we only consider singly charged ions, it is reasonable to assume that the Bethe-Heitler and Trident pair production processes that are mediated by Coulomb fields can be safely neglected. The Trident process, however, can also arise in strong electromagnetic fields. In fact, it is now increasingly viewed as comprising two processes [38]: the one-step process, in which the intermediate photon is virtual (traditionally referred to as the electromagnetic Trident [11]), and the two-step process, in which a real gamma-ray photon is emitted via nonlinear Compton scattering and creates a pair via the multiphoton Breit-Wheeler process. Only the latter two-step process is described in our PIC-QED simulations. Although the one-step process can be significant for ultrashort laser pulse durations (i.e., less than the laser period, T L = 2π/ω L 3.3 fs), it has been demonstrated in references [39,40] that it becomes negligible compared to the two-step process under conditions relevant to the forthcoming multi-petawatt lasers (i.e., I L 10 23 W cm −3 and 5-10T L pulse durations). Therefore, the succession of nonlinear Compton scattering and multiphoton Breit-Wheeler processes (i.e., the two-step Trident process) should dominate pair production for the laser parameters considered in this work.
The simulation box is defined by a spatial grid of dimensions 200 × 57.6 μm 2 using 20 000 × 5760 mesh cells, and the target is represented by 9.2 × 10 8 macroelectrons and the same number of macroions. Radiation reaction (RR) is modeled in both the classical (deterministic) and quantum (stochastic) regimes [37]. Throughout the text, for any vector quantity, the subscripts ' ' and '⊥' shall refer to the components in the laser polarization plane that are, respectively, parallel and perpendicular to the (rotating) laser electric field, while the subscript 'x' refers to the component parallel to the laser propagation axis. The momentum components of a particle can thus be written as ⎛ where φ is the phase of the electromagnetic pulse, and primes denote quantities computed in the frame moving at the laser-energy propagation velocity, hereafter called the comoving frame.

Overall interaction scenario
Early in the interaction, the plasma electrons are pushed ahead in the laser's rising edge by the ponderomotive force. The resulting charge separation between the electrons and ions induces a longitudinal plasma wave. This wave has an electric field E x > 0 within its first half-period from the electron density peak, and so accelerates the ions in the forward direction. When the associated electrostatic potential exceeds the ion kinetic energy in the frame comoving with the plasma wave (or, equivalently, the laser front), most of the ions get reflected. A quasi-steady double-layer structure then arises at the front edge of the laser pulse, which acts effectively as a piston upon the traversed plasma [35], as shown in figure 1(a), which plots at time t = 45T L the longitudinal lineouts of the electron (n e ) and ion (n i ) number densities, both normalized to n c , along with those of the longitudinal (E x ) and transverse (E y ) electric fields measured in units of the laser field strength E L (and multiplied by 300 for visibility). The laser front is then located at x 28 λ L . When RR is turned off, this structure can be maintained over relatively long times (depending on the laser pulse duration), as shown in figure 1(b). The peaks in n e and n i at x 28.5λ L and x 27.5λ L , as well as the coherent E x > 0 field induced in-between, are clear signatures of the double layer at the laser edge. When RR is enabled, and at laser intensities I L 10 23 W cm −2 , i.e., high enough to initiate the radiation-dominated interaction regime [41], the plasma wave excited at the laser edge is strongly altered: the electron backward motion (towards the laser field) is hindered by the RR force [42], which tends to enhance the potential barrier, and thus facilitate ion reflection off the laser piston [35]. However, when the laser intensity exceeds I L 5 × 10 23 W cm −2 , as is the case here, rapid reflection of the returning electrons may not necessarily happen due to the straggling effect, which, in the present context, is the possibility for electrons to penetrate deeper into the laser wave, and emit photons more energetic, than would be allowed classically [43,44]. The potential barrier at the laser edge may then drop below the ion reflection threshold, hence quenching the local piston structure. Eventually, the electrons leaking into the laser wave will find themselves trapped at the electric-field nodes around the pulse maximum [36]. As they are carried along with the laser pulse, they will accumulate, and so will the ions subject to the charge-separation field, causing a highly compressed, quasineutral plasma bunch to form within the laser pulse, as is clearly seen in the density profiles of figure 1(a). The laser edge is then located at x 28λ L , where the electron density shows a first peak. Figure 1(c), which displays the longitudinal (x − p x ) particle phase spaces, reveals that the E x > 0 field (green line) induced just behind the laser edge drives forward the ions (black dots) to a normalized momentum γβ x 0.7, yet without causing ion reflection. By contrast, the highly compressed (n e n i 100n c ) plasma bunch formed within the laser pulse (24 x 26λ L [figure 1(a)] carries a potential barrier capable of reflecting, at its backside (x 24λ L , where the laser field peaks, see blue dashed line) a significant fraction of the ions, thus creating a shock-like structure. A similar, albeit denser, structure is observed to arise when running the same simulation in 1D3V geometry [see figure 1(d)].
The overall interaction dynamics is illustrated more quantitatively in figures 2(a) and (b), which plot the temporal evolution of the particle densities around the laser pulse maximum [figure 2(a)] and of the maximum strengths of the E x and E y electric fields in the simulation domain [figure 2(b)]. The longitudinal profiles of those quantities and particle phase spaces are displayed at time t = 45T L in figures 1(a) and (c). The interaction dynamics is found to proceed through the following three main stages (indicated by numbers in figure 2).
(1) For t 13T L , the electron density near the laser peak starts increasing (up to n e 100n c ) due to radiative cooling [45], while the ion density remains unperturbed.  . The compressed bunch is clearly visible around (x, y) (40, 0)λ L , leaving behind a tenuous plasma channel. This indicates that the bunch eventually turns into a highly reflective piston, as corroborated by the phase space of figure 3(b), recorded at the same time. A reflection coefficient of 80% is measured in the 45 t 65T L time interval. The reflected deuterons then attain longitudinal momenta as high as γβ x 3. During this stage of the interaction, the laser pulse undergoes strong depletion [46], as shown by the steady decrease in the E y field maximum [ figure 2(b)]. We also note that from figure 2(b) the positron density has then reached a saturated value, n p n c .

Analytical model
Here, we propose a model for the dense plasma bunch formed within the laser pulse, taking explicit account of the effect of RR on the laser-plasma interaction.
In the following, we shall assume that the bunch rest frame coincides with the frame moving at the laser-energy propagation velocity. Since the bunch behaves like a piston, its propagation speed, β b , can be estimated in the standard way, that is, by equating the radiation pressure with the ion momentum flux [47,48]: In the above, n i0 is the (unperturbed) ion density, m i the ion mass, m e the electron mass, and R the laser reflection coefficient in the comoving frame. In the interaction region of interest here, the strong radiative losses experienced by the electrons lead to R 1, so that β b 0.69, in very good agreement with figure 2(c). The kinetic energy of the reflected ions is predicted to be E i,r = 2m i c 2 B 2 /(1 + 2B) 6900m e c 2 , corresponding to a longitudinal momentum p i,x 5900m e c, consistent with figure 3(b).
An alternative estimate of the effective laser propagation speed can be obtained from the dispersion relation of a circularly polarized electromagnetic wave through an absorbing (due to radiative losses) unmagnetized electron-ion plasma. Neglecting thermal effects, this relation can be written as [42,49] where ω pν = n ν Z 2 ν e 2 / (m ν 0 ) and γ ν are the (nonrelativistic) plasma frequency and mean Lorentz factor of particle species ν, in the comoving frame. See the appendix for a derivation of equation (3b) from the simplified Landau-Lifshitz equation. The G term quantifies the dissipation induced by RR. It is defined in terms of the (Lorentz invariant) radiated power P γ as G = P γ γ e m e c 2 (Ω) , ( 4 a ) where χ e is the (mean) electron quantum parameter, α f = e 2 /(4π 0 c) 1/137 is the fine-structure constant, and λ c = h/m e c the Compton length. First-order quantum corrections to the radiated power are contained in the Gaunt factor g(χ e ) 1 + 4.8 (1 + χ e ) ln (1 + 1.7χ e ) + 2.44χ 2 e −2/3 [50]. It is worth noting that Bashinov and Kim [42] used the Lorentz-Abraham-Dirac [51] equation to compute G (written down as δγ 3 e in their paper). Thus, they obtained a Ω 2 term close to ours, except for the substitution G → δγ 3 e where γ e is defined through a L = γ e 1 − iδγ 3 e v ⊥ . A direct consequence of this relatively simple formulation is that Ω 2 ∼ Ω 2 when δγ 3 e ∼ 450. This contrasts with our assertion that for a L = 850, Ω 2 < Ω 2 . Such a discrepancy is due to the fact, that in our case, the longitudinal motion of ions (Γ b > 1) as well as the quantum correction g (χ e ) to the radiated power (4b) have been accounted for, leading to Ω 2 < Ω 2 for a L = 850. To make analytical progress, we rewrite the real part of Ω involved in equation (4a) as where (Ω 2 ) = ω pe 2 /[γ e (1 + G 2 )] + ω pi 2 /γ i ω 2 pe [1/γ e (1 + G 2 ) + Z i m e /(m i γ i )], assuming quasineutrality within the bunch (n e Z i n i ). We have also introduced X = (Ω 2 )/ ( . We now assume that the laser comoving frame corresponds to the group velocity deduced from the above dispersion relation, where the frequency (ω) and wavenumber (k) are taken to be real and complex, respectively. While in the general case of an absorbing medium the group velocity differs from the effective energy propagation velocity [52][53][54], we expect this difference to remain weak provided (Ω 2 ) (Ω 2 ), as is relevant to our study. The wavenumber can be expressed as k = (n + iκ)ω/c, where n + iκ is the complex refraction index. Introducing the transverse plasma permittivity, (ω) ≡ (n + iκ) 2 = 1 − Ω 2 /ω 2 , one has , the plasma bunch Lorentz factor. There follows Under our conditions, H remains 1, so that β g √ ( ). In the frame moving at the laser group velocity, which is assumed to coincide with the bunch rest frame, the laser wavenumber vanishes (k L = 0), and therefore so does the laser magnetic field. As a consequence, the laser electrons experience only a spatially uniform laser electric field, rotating in time [55], at the frequency ω L = (Ω), given by equation (5). Thus, the mean quantum parameter χ e can be written as To estimate p e,⊥ , we further suppose that, in the comoving frame, the gyrating electron motion is purely transverse (p e,x = 0). This configuration was addressed in reference [56]; using their equations (31) and (32), one can derive the approximate formula: The dimensionless parameter rad = e 2 (Ω) /(6π 0 m e c 3 ) 6.2 × 10 −24 (Ω) characterizes the RR effect [56]. The mean quantum parameter can then be expressed as where we have approximated √ (ω 2 ) ω L /Γ b . Performing a 3rd-order Taylor expansion of g(χ e ), the above equation simplifies to 2.3χ 4 e + 9.9χ 3 e − 2.7χ 2 e + χ e where χ 0 = a L / rad ( ω L /m e c 2 ) is the value of the electron quantum parameter without quantum (i.e., g (χ e ) = 1) and radiative (i.e., G = X = U = 0) corrections, and where longitudinal motion is absent (i.e., Γ b = 1). For our parameters (a L = 850, Γ b 1.4), the solution to equation (13) is χ e 0.4, consistent with the χ e distribution measured within the bunch [see inset of figure 1(c)].
To conclude this part, knowing χ e , one can evaluate A, and then, from β g √ ( ), deduce the bunch density in the lab frame (noting that n e = Γ b n e ): Setting γ i 1, the above formula predicts n e 500n c , consistent with our 2D and 1D simulation results [see figures 1(a) and (d)].

Positron acceleration
From t 20T L , the multiphoton Breit-Wheeler process causes copious production of electron-positron pairs at the center of the laser pulse, saturating at a maximum positron density n p n c by t 45T L [see figure 2(a)]. Their typical energy depends on that of the γ-ray photons produced by nonlinear inverse Compton scattering. The latter is characterized by the quantum emissivity F(χ e , χ γ ) [with χ γ ( ω γ /m e c 2 γ e )χ e the photon quantum parameter], which itself depends on the reduced variable y ≡ 1/γ 2 e γ e m e c 2 / γ e m e c 2 − ω γ 2/3 ω γ /ω e 2/3 where ω e = |p e × F L |/p 2 e = m e cω L a L p e,⊥ /p 2 e is the instantaneous electron rotation frequency [50,57]. Assuming that most of the high-energy photons fulfill y 1 [58], and that the average positron Lorentz factor is γ p ω γ /(2m e c 2 ), the mean energy of the positrons (at creation time) is expected to be Under our conditions, this formula predicts an approximate positron energy of γ p 400.
Assuming that the positrons are created with vanishing longitudinal momentum (p x = 0), it has been shown [56] that the positron momenta in the laser polarization plane satisfy Using There follow the parallel and perpendicular positron (normalized) velocities: Once created, the positrons get accelerated in the forward direction, as shown in the phase space of figure 3(b) (green dots). This is due to the combined effect of the longitudinal electrostatic field (E x ), associated with the plasma density gradients, and the radiation-modified Laplace force, originating from the particles' developing a momentum component parallel to the laser electric field [59,60]. The energy-angle To explain this difference, we estimate the longitudinal force experienced by the positrons in the laboratory frame as where The second term in the right-hand side of equation (21) corresponds to the radiation-modified Laplace force [59,60]. The energy gain of the positrons resulting from the longitudinal electrostatic and radiation forces is where l p denotes the acceleration length. To take into account the spatial profile of the laser electric field, we replace a L byā L = 1 ln 2 and τ 11.4T L the laser pulse FWHM duration). Finally, the maximum energy gain of the forward-accelerated positrons can be expressed as For a numerical application, we take l p 4λ L ,ā L 765 and a x ≡ eE x /m e cω L 0.2a L as estimated at t = 65T L . The above formula then predicts ΔE p 4 GeV, consistent with the maximum positron energies found in figure 4(a). Although the energy gain is due to both F rad,x and eE x , the latter charge-separation field tends to reduce the angular emission of the forward-accelerated positrons, see figure 4(a). Note that for a vanishing longitudinal motion in the simulation frame (β b = 0), the radiation-modified Laplace force becomes zero, thus lowering the maximum energy gain by a factor ∼ 2-3.

Influence of the plasma ion mass
Finally, it is of interest to examine the influence of the ion mass on the above-discussed phenomena. To this purpose, we compare in figures 5(a)-(d) the electron and ion density profiles extracted from simulations run with either a deuterium (a) and (b) or a hydrogen (c) and (d) plasma of same electron density. While a dense plasma bunch is also observed to form around the laser axis in the H case, it is less compressed and of narrower extent than in the D case, and also located closer to the foot of the laser pulse. We ascribe these differences to the reduction in the RR effects with decreasing ion mass. Indeed, those are quantified by  albeit modestly, the laser energy conversion efficiency into γ-ray photons (from ∼ 40% to ∼ 30%), and more strongly, into the positrons (from ∼ 0.2% to ∼ 0.1%, along with a sixfold decrease in the positron number), see figure 6(a). These variations are consistent with previous works [61], yet they are not accompanied by markedly different laser energy conversion efficiencies into electrons and ions, the temporal evolution of which is plotted in figure 6(b). While the photon spectrum is significantly colder in the hydrogen plasma [figure 6(c)], the high-energy parts of the positron spectra prove surprisingly similar in the two cases [figure 6(c)]. Note also that the mean energy of the positrons ( 200 MeV) and its weak dependence on the ion mass agree with the approximate formula (15).

Conclusions
By means of 2D PIC simulations and analytical modelling, we have examined in detail how RR alters the interaction of a 10 24 W cm −2 , 30 fs laser pulse with a 10-times critical density plasma. Our analysis has revealed that RR induces the formation of a piston-like, dense plasma bunch within the laser pulse. Through their interaction with the laser field, part of the high-energy photons radiated by the bunch electrons subsequently convert into electron-positron pairs. Drawing upon previous works on radiation-modified laser-plasma interactions [15,42,56], we have developed an analytical model which predicts to satisfactory accuracy the mean density, velocity and radiative efficiency of the plasma bunch. Moreover, we have shown that, as a result of the charge-separation field and the radiation-modified Laplace force, the positrons undergo efficient forward acceleration so that their final GeV-range energy is about twice that of the laser-driven electrons. Those results, evidenced in a deuterium plasma, turn out to be sensitive to the ion mass: the use of a hydrogen plasma leads to a less stable and dense plasma bunch, along with less efficient pair generation. A brief remark is in order regarding the impact of the reduced 2D geometry of our simulations. Since circular laser polarization has been considered, one expects the electron-ion bunch to be approximately cylindrically symmetric. Although adding a third spatial dimension may somewhat alter the level of electron energization, and therefore the strength of RR, it should yield little qualitative change to the sequence of processes leading to the bunch formation, similarly to what we observed when moving from 1D to 2D simulations.
Finally, the results of this work will be useful for the design of high-field plasma experiments at the soon-to-be-operational multi-petawatt laser facilities [3][4][5], especially in view of producing dense pair plasma beams of relevance to laboratory astrophysics [62].

Appendix. Electromagnetic dispersion relation in the radiation-dominated interaction regime
Let us consider a circularly polarized, electromagnetic (laser) plane wave, propagating through a plasma along x > 0, and characterized by the potential vector where φ = ω L t − k L x is the phase. Using Maxwell's equations and the Lorenz gauge, one readily obtains the wave equation satisfied by A: where j = j e + j i = −en e v e + Zen i v i is the total current density, with contributions from the plasma electrons (j e ) and ions (j i ). We have introduced e the elementary charge, n e and n i the electron and ion numeber densities, v e and v i the electron and ion mean velocities, and Z the ion charge state. We now consider the comoving frame of the laser wave, of velocity β g in the laboratory frame. In the following, primes will denote quantities measured in the comoving frame.
In the comoving frame, the laser frequency and wavenumber fulfill ω L = (Ω) and k L = 0, resulting in B L = 0. In this frame, the laser wave thus turns into a purely time-dependent, spatially homogeneous, rotating electric field [55]. The electron and ion equations of motion can then be written as dp e dt = e dA dt − P γ γ e m 2 e c 3 p e (A.3a) dp i dt = −Ze dA dt (A.3b) where P γ is (the Lorentz invariant) radiated power defined in equation (4b) and A = A has been used. where Ω is the (Lorentz invariant) effective plasma frequency defined by: