Impact of the laser spatio-temporal shape on Breit-Wheeler pair production

The forthcoming generation of multi-petawatt lasers opens the way to abundant pair production by the nonlinear Breit-Wheeler process, i.e., the decay of a photon into an electron-positron pair inside an intense laser field. In this paper we explore the optimal conditions for Breit-Wheeler pair production in the head-on collision of a laser pulse with gamma photons. The role of the laser peak intensity versus the focal spot size and shape is examined keeping a constant laser energy to match experimental constraints. A simple model for the soft-shower case, where most pairs originate from the decay of the initial gamma photons, is derived. This approach provides us with a semi-analytical model for more complex situations involving either Gaussian or Laguerre-Gauss (LG) laser beams. We then explore the influence of the order of the LG beams on pair creation. Finally we obtain the result that, above a given threshold, a larger spot size (or a higher order in the case of LG laser beams) is more favorable than a higher peak intensity. Our results match very well with three-dimensional particle-in-cell simulations and can be used to guide upcoming experimental campaigns.


Introduction
Abundant electron-positron pair production is among the most exotic and striking processes of extremely high-intensity laser interaction with light and/or matter that will be made possible by the emerging class of multi-petawatt (multi-PW) laser systems [1]. Under conditions envisioned for forthcoming facilities such as Apollon [2], CoReLS [3], ELI [4], OMEGA-EP OPAL [5], XCELS [6] and zeus [7], the dominant process leading to pair production is nonlinear Breit-Wheeler, i.e. the decay of a high-energy photon, as it interacts with the strong laser field, in an electron-positron pair [8,9].
In the last decade, various simulation campaigns have been conducted to help design the future experiments to efficiently drive this process, an effort that was largely made possible and supported by particle-in-cell (PIC) simulations accounting for quantum electrodynamics (QED) processes [10][11][12][13][14][15]. Several configurations have been considered. Electron-seeded electromagnetic/QED cascades in the collision of two counter-propagating laser pulses, a setup originally proposed by Bell and Kirk [16], have attracted particular attention [11,14,[17][18][19][20] as such cascades were identified as a possible limitation on the attainable intensity of high power lasers [21]. This setup was further built upon considering either the use of plasma channels [22] or the collision of multiple laser pulses [23][24][25]. All these works considered multiple laser beams, but electromagnetic cascades were also predicted considering the direct irradiation of a solid target by a single, extremely intense laser pulse, a situation that leads to the production of dense pair plasmas [26] (see also [27] for an application to laboratory astrophysics). In these cases, the solid target, an overdense plasma acts as a mirror and pair production is efficiently achieved in the field of the incident and reflected laser light at the target front [28].
All the previously discussed configurations considered cascades seeded by electrons initially at rest in the laboratory frame. The seed electrons, strongly accelerated in the electromagnetic fields of laser pulses with extreme intensity (typically ∼ 10 24 W cm −2 ), gain energies at the GeV-level (as expected, at such intensities, from the ponderomotive scaling) and are the source (via nonlinear Compton scattering) of the high-energy photons triggering pair production. Another possible path toward pair production in the laboratory relies on an already existing high-energy electron or photon source. The two latter schemes are related, as photons with energy of the order of the electron one are produced via nonlinear Compton scattering of electrons in the intense laser fields [29][30][31][32][33]. Indeed, this configuration was exploited in the seminal E-144 experiment at SLAC [34], where the head-on collision of 46.6 GeV electron beams with TW laser pulses gave the first direct demonstration of nonlinear Breit-Wheeler pair production. However, the laser intensity in this experiment was only marginally relativistic, and very few positrons were produced (∼ 100 positrons over nearly 22 000 electron-laser collisions).
Combining GeV-electron beams with PW and multi-PW lasers will allow to extend this study deep into the relativistic regime, as investigated in [35] that numerically demonstrated the possibility of producing extremely bright high-energy photon sources as well as positrons bunches for laser intensities 10 22 W cm −2 , in conditions relevant to the Apollon facility. This configuration was further studied theoretically and via PIC simulations in [36][37][38], while the possibility to produce pairs via nonlinear Breit-Wheeler with high-energy photons from a Bremsstrahlung source was presented in [39].
In this paper, we focus on the study of nonlinear Breit-Wheeler pair production following from the head-on collision of an extremely intense laser (intensities in the range 10 21 -10 25 W cm −2 ) with a burst of gamma photons (with energies ranging from 100 MeV to few 10 s of GeV). In contrast to [39], where the authors focused on the optimization of the photon source, we consider here the high-energy photon burst as given (only its energy will vary), and do not discuss its origin (various sources of high-energy photons have been proposed [25,[39][40][41][42][43][44][45][46]). Rather, we aim to optimize the conditions of interaction with the colliding high-intensity laser. Motivated by experimental constraints, we investigate the optimal conditions for pair production varying the laser polarization, focusing, spatial or temporal profiles, always considering a fixed laser energy.
In particular, special attention will be paid to the use of Laguerre-Gauss (LG) beams, which have recently attracted the interest of the ultra-high intensity laser-plasma community [38,[47][48][49][50] for their ability to carry orbital angular momentum (OAM) [51]. These beams, that emerge as eigen-modes of the paraxial equation, have unique properties, such as a ring-shaped intensity distribution and the large OAM, that could have an impact on pair production. It was for instance demonstrated in [38] that the collision of an electron beam with a LG beam can lead to efficient production of gamma rays carrying large OAM, to enhanced secondary radiation emission and pair production. However, that study was performed at constant maximum intensity, so it is hard to distinguish the impact of the increased energy (up to 3× higher for the LG beam) from the role of the laser spatial profile itself. We will discuss this aspect in detail in this paper.
Because of the remarkable number of free parameters to study (polarization, focusing, spatial and temporal profile) and the complex spatio-temporal dependence of the electromagnetic fields of LG beams, it is useful to develop and validate a reduced analytical model for pair creation (see reference [52] for an analytical study of the Breit-Wheeler pair production process in a Gaussian beam). Starting from the probability for a given high-energy photon to decay in an electron-positron pair while it crosses an extremely intense laser pulse, we propose a simple, yet accurate, model to compute the number of pairs created during this interaction. The computation relies on three assumptions. (i) The electron-positron pair production rates are obtained in the locally constant cross-field approximation [53]. (ii) The newly created charged particles do not further contribute to pair production. (iii) The potentially complex electromagnetic fields explored by the high-energy photons as they cross the laser pulse are approximated, assuming that each photon sees a succession of single half-periods of monochromatic plane waves, the field amplitude of each half-period being computed so that both the temporal and spatial profiles of the laser pulse at focus are accounted for. As will be further discussed in this work, these assumptions are well satisfied for a wide range of parameters of interest for this study and currently accessible laser parameters. Our model goes further than previous studies [36,53] as it accounts for arbitrary polarization and spatio-temporal profile of the laser beam, and it is not limited to small pair production probability and/or rate.
Three-dimensional PIC simulations performed with the code SMILEI [54] and embarking the relevant QED modules [15] are presented and validate our model. Assumptions (ii) and (iii) are shown to hold under conditions relevant to ultra-short laser pulses (up to few tens laser cycles) typical of the Ti:sapphire technology. Assumption (ii) holding is of particular importance as it demonstrates that, under conditions relevant to forthcoming ultra-short laser facilities, the high-energy-photon-laser interaction proceeds  (19)) on the maximum photon quantum parameter χ 0 (as defined by equation (13)), for a linearly polarized (LP) background field (ε = 0, black line) and a circularly polarized (CP) background field (ε = ±1, blue line). Dashed lines show the asymptotic behavior given by equations (22) and (23) for the LP case. The inset highlights the values of χ 0 for which LP (χ 0 < 2.52) or CP (χ 0 > 2.52) gives higher pair production probability and the relative difference of I ε (in absolute value) between the two cases. essentially in a regime which we refer to as the soft-shower regime [55]. Furthermore, our PIC simulations allow to test our model against complex laser pulse geometries, e.g. considering LG beams, showing that assumption (iii)-which is at the core of the simplification of our model-provides accurate predictions even for non trivial electromagnetic field configurations.
The proposed model and performed study make it possible to clearly identify the laser parameters and high-energy photon energy that optimize pair production for a given laser energy. This work also highlights the conditions in which either stronger focusing or increased gamma photon energy does not improve the production of primary pairs. The paper is structured as follows. In section 2, we introduce the key properties of the nonlinear Breit-Wheeler process, and we propose a reduced model for pair creation in the head-on collision between gamma photons and a plane wave with arbitrary polarization. The possibility to account for a time envelope of the laser field is also considered and the model is benchmarked against 1D PIC simulations. In section 3 we consider the effect of the spatial dependence of the laser field. Both LG and Gaussian beams are considered and we introduce the notion of a total cross-section for which we propose a reduced model, shown to be in very good agreement with 3D PIC simulations. The optimal conditions for pair creation with respect to the laser geometry and intensity are also discussed. In section 4 we then use our model to discuss pair production on upcoming Ti:sapphire laser facilities, before giving our conclusions in section 5.

Reduced model for pair production in a time-varying background field
This section provides the theoretical framework for describing nonlinear Breit-Wheeler pair production, and aims to establish a simple model for pair production in a time-varying background field. Throughout this section, we consider the head-on collision of high-energy gamma photons with a strong, optical, laser pulse (also referred to as the background field). In addition to exact expressions, we provide practical formulae relevant to forthcoming high-power laser facilities, focusing on optical laser systems with micrometric wavelength λ, with relativistic field strength a 0 = eE 0 /(m e cω) 1 (with ω = 2πc/λ the background field angular frequency), corresponding to I 0 λ 2 well beyond 10 18 W cm −2 μm −2 , with I 0 = 0 cE 2 0 /2 the laser intensity. Note that SI units are used throughout this work: e, m e and c denote the elementary charge, electron mass and speed of light in vacuum, respectively, 0 the vacuum permittivity, r e = e 2 /(4π 0 m e c 2 ) the classical radius of the electron, τ e = r e /c the time for light to cross r e , the reduced Planck constant and α = e 2 /(4π 0 c) the fine-structure constant.

Pair production instantaneous rate and probability
Nonlinear Breit-Wheeler is the dominant QED process leading to electron-positron pair production from the collision of a high-energy photon with a strong (optical) laser pulse [56,57]. As discussed in reference [53], for a relativistically intense background field, a 0 1, the local constant field approximation holds and the pair production probability can be computed by integrating the instantaneous pair production rate. For a high-energy photon propagating with velocity c, normalized energy γ γ = ω γ /(m e c 2 ) and instantaneous quantum parameter where E S = m 2 e c 3 /(e ) 1.32 × 10 18 V m −1 is the Schwinger field and E and B denote the electric and magnetic fields at the photon position. The instantaneous pair production rate is given by (see, e.g. references [58,59]): where W 0 = 2α 2 /(3τ e ). For a fixed value of the quantum parameter χ γ , the rate is inversely proportional to the gamma photon energy γ γ , and depends non-trivially on χ γ through the function with μ = 2/[3χ γ ξ(1 − ξ)] and K n (x) the modified Bessel function of the second kind. This function, shown in figure 1(a), can be well approximated (within less than 1% from the numerically computed values for . (4) The denominator in equation (4) provides a simple improvement of the well known Erber approximation [60], which also allows to recover the correct asymptotic behavior (dashed grey lines in figure 1 with c 1 = (3/16) 3/2 3/2 0.344 and c 2 = 3 2/3 45Γ 4 (2/3)/(56π 2 ) 0.569 (Γ(x) denoting the gamma function).
The quantity W BW defined by equation (2) is the instantaneous rate of pair production, so that considering, at a time t 0 , a given number N 0 of identical high-energy photons entering a strong background field, the number of photons N γ remaining at later times t > t 0 satisfies the rate equation: where both N γ and W BW depend on time (the latter when the high-energy photons explore electromagnetic fields with spatial and/or temporal variations). Note that equation (7) does not account for subsequent radiation from the produced pairs, which will be verified to be accurate for the investigated parameter range in section 3. The solution to this equation reads from which one can derive the probability P(t 0 , Δt) = 1 − N γ (t 0 + Δt)/N 0 for a given high-energy photon to decay into an electron-positron pair between time t 0 and t 0 + Δt: When the time-integrated rate (i.e. the integral term in equation (8)) is small, it can be assimilated to the probability P(t 0 , Δt), as P(t 0 , Δt) → t 0 +Δt t 0 W BW t dt 1. This leads some authors to refer to W BW as the pair production probability per unit time. While this is correct in the limit P(t 0 , Δt) 1, it leads to inconsistencies as the probability increases and the time-integrated rate assumes values greater than unity (see, also, references [53,[61][62][63]). Hence, throughout this work, we will distinguish pair production rate (as defined by equation (2)) and probability (as defined by equation (9)).

Pair production probability in an arbitrarily polarized plane wave
Let us now focus on the case of a high-energy photon colliding head-on with an ultra-intense laser pulse described by a plane wave (propagating in theẑ-direction), with electric field: magnetic field B =ẑ × E/c, and polarization ellipticity ε ∈ [−1, 1]. It takes the particular values ε = ±1 for a CP wave (± standing for opposite helicities), and ε = 0 for a LP wave (x denoting by convention the direction of polarization). It is important to stress that our choice of parametrisation (equation (10)) ensures that the laser energy is the same for any value of ε. On the contrary, E 0 / √ 1 + ε 2 measures the maximum laser electric field amplitude and varies with the ellipticity.
When a high-energy photon collides with such a laser field (we define t 0 = 0 as the time at which the photon enters the background field), the value of its quantum parameter evolves in time as the photon explores regions with different electric and magnetic field amplitude (here z γ (t) denotes the photon time-dependent position). In the case of head-on collisions (c = −cẑ), the time-dependent photon quantum parameter reduces to and Note that, with our convention on the background field parametrization (equation (10)), χ 0 denotes the maximum quantum parameter for the LP case only. For arbitrary ε, the maximum achievable quantum parameter χ m is a decreasing function of |ε| Note also that equation (12) takes a simple form when considering either a CP or a LP background field:

Decay probability after propagating through half a period of the background field
For arbitrary values of ε, χ γ (t) is either a constant (for ε = ±1) or a periodic function of time with period τ /4 (we recall τ = 2π/ω) corresponding to the time for the high-energy photon to cross half a wavelength of the background field. It is thus convenient to compute the probability P m of the high-energy photon to decay into a pair during this interval of time. Let us then consider a time t m at which the high-energy photon sees a maximum of the background electric field, and compute the probability P m for the high-energy photon to decay into an electron-positron pair in between the times t 0 = t m − τ/8 and t 0 + τ/4 = t m + τ /8: where R m τ/4 is the time-integrate rate, while R m denotes the average rate. Let us stress that the dimensionless quantity W 0 /(2ωγ γ ) in equation (17) is of order 1 for optical background fields (with micrometric wavelength λ) and GeV-level high-energy photons: The probability P m , that measures the contribution of a single field maximum, depends on the integral quantity which, for a given polarization ellipticity ε, is a function of χ 0 only as shown in figure 1(b) for LP (black line) and CP (blue line). This integral is the key object to compute the pair creation probability and can be calculated either numerically or by means of an approximated formula as described in the following. In appendix A, we generalize the approach proposed in [36] to arbitrary ε and χ, and show that I ε (χ 0 ) can be written in the approximate form: where χ m is defined in equation (14) and erf(x) = 2 √ π x 0 e −t 2 dt is the error function. The function F(s) emerges from the saddle point approximation used to compute the integral when the main contribution comes from the vicinity of the field maximum. It varies slowly with s, as F(s) 2 √ π s for s 1 and F(s) → 1 for s → +∞. In equation (20), F(s) takes for argument s ε (χ m ), with where b 0 (χ) denotes the derivative of b 0 (χ) and c(χ) is a slowly varying function of χ, as c(χ) Moreover it allows to recover the asymptotic behavior of equation (19) in the limiting cases of small and large quantum parameter χ 0 : where we used equations (5) and (6), and c 3 = 27 π/2/64 0.529, c 4 = 9 √ 3π/64 0.765 and c 5 = 3 2/3 45Γ 4 (2/3)/(56π) 1.789.

Total decay probability
As will be further discussed in section 2.3, considering parameters typical of the upcoming generation of laser facilities, the probability for a high-energy photon to decay into a pair after crossing a single half-wavelength of the background field will in general be small compared to 1. This will not however be the case when considering the cumulative effect of crossing several wavelengths.
Let us denote as t 0 = 0 the time at which the photon starts to interact with the background field 6 . The probability for a high-energy photon to decay into a pair after t = nτ /4 of interaction with the background field (or, equivalently, after crossing n local maxima of the field amplitude) P tot (t = nτ /4), can be derived from the number N γ of photons surviving after a time t: Here we have used m as a running index for compactness, and P m (given by equation (17)) andP m denote the probabilities for the photon to decay and not to decay, respectively, during the mth interval. We can then write the total probability of producing a pair as R denoting the average pair production rate. Considering a background field with sin 2 temporal profile with maximum field amplitude (c) a 0 82 corresponding to χ 0 = 0.5, (d) a 0 660 corresponding to χ 0 = 4. Solid lines correspond to the exact probability obtained by integrating numerically equation (9), dashed lines and dots to the theoretical prediction obtained from equation (25), for LP (ε = 0, black) and CP (ε = ±1, blue) cases. The red dash-dotted lines correspond to the theoretical prediction from equation (25) using the approximation of equation (20) for LP.
Let us note that this formula is exact at t = nτ /4 with n ∈ N, and-as will be shown later-can also be used with good approximation at all times t τ/4. Moreover, in the case of a monochromatic plane wave (no temporal envelop), the average pair production rate simply reduces to R = R m . When dealing with a pulse with a finite temporal envelop, however, equation (25) needs to be computed combining the contribution of successive (τ/4)-long intervals, using the local maximum of the background field strength for each time interval. This approach gives very good results even considering ultra-short, few cycles laser pulses. To confirm this, figure 2 shows the temporal evolution of the decay probability for a high-energy photon with γ γ = 10 3 colliding head-on with different background fields. In panels (a) and (b), the high-energy photon interacts with a background field with constant amplitude a 0 82 (χ 0 = 0.5, figure 2(a)) and a 0 660 (χ 0 = 4, figure 2(b)). The interaction lasts for a time τ , i.e. the photon explores two wavelengths of the background field. In panels (c) and (d), the high-energy photon interacts with a background field with a sin 2 temporal profile in intensity, a full-width-half-maximum (FWHM) of 5τ , and maximum field strength a 0 82 and a 0 660, respectively.
In all panels, the solid lines denote the exact probability computed by numerically integrating equation (9), considering either LP (black lines) or CP (blue lines) background field. We can see that, for χ 0 = 0.5, CP produces less pairs than LP, while for χ 0 = 4 CP is slightly more efficient (the difference between the two polarizations will be discussed in details in section 2.3.1). These probabilities are compared against the prediction of equation (25), computed using either the numerically evaluated I ε (χ 0 ) (dashed lines in panels (a) and (b), dots in panels (c) and (d)) or the approximation given by equation (20) to compute R m (red dot-dashed lines, only computed for LP as it is exact for CP). Note that here and in the following, whenever equation (20) is used, the improved Erber approximation equation (4) is exploited. Panels (a) and (b) show an excellent agreement between the exact computations (solid lines) and equation (25) (dashed lines), both using the numerically integrated values of I ε (χ 0 ). A remarkably good agreement is also found when considering the fully analytical approximation of equation (20) for I ε (χ 0 ) (red dot-dashed lines). The probability is slightly overestimated in this case, as expected for χ 0 ∼ 1 for which equation (20) shows the greatest departure from the exact expression of I ε (χ 0 ). Similarly, an excellent agreement between all three approaches is found in panels (c) and (d) considering a short (5τ FWHM in intensity) background field. In this case the value of R m in each τ /4 interval is computed by taking the value of I ε given by the local maximum of the field, i.e. the time variation of the envelope is considered slow during τ /4. The validity of this assumption is confirmed by PIC simulations discussed in the following. In particular, figure 2 shows that the approach leading to the derivation of equation (25), which consists in treating the cumulative contributions of successive maxima of the background field, remains a very good approach even for ultra-short, few cycle, background fields.
Last, we note that the results of 1D PIC simulations performed with Smilei (see appendix B for details) are not distinguishable from the exact computations (solid lines) for all reported cases, and are therefore not shown in figure 2.

Influence of the background field polarization
Interesting insights into the role of the background field polarization on pair production can be gained from the asymptotic behavior provided by equations (22) and (23).
We recall that, motivated by experimental constraints, we compare LP and CP at fixed energy: this condition implies that the CP beam has a lower maximum amplitude. Because of this, for χ 0 1, equation (22) shows that any departure from the LP case leads to a tremendous decrease of I ε (χ 0 ). Indeed the exponential cut-off of the pair production rate is strongly affected by the decrease of χ m with |ε| as shown in equation (14). In contrast, for χ 0 1, f (ε)/(1 + ε 2 ) 1/3 is an increasing function of |ε|: the reduction of the field peak amplitude is compensated by the fact that the absolute value of the amplitude is not time dependent (see equation (15)). Hence, in this range, increasing the background field ellipticity increases the pair production rate. However, the rate increase from the LP to the CP case is small (less than 12%), and it is thus of marginal use to consider CP for boosting pair production 7 .
The differences between the LP and CP cases are clearly shown in the inset in figure 1(b), where I 0 (χ 0 ) is found to be orders of magnitude larger than I ±1 (χ 0 ) for χ 0 1 while I ±1 (χ 0 ) is only about 10% larger than I 0 (χ 0 ) for χ 0 1. Similar conclusions can be drawn from figure 2. For small values of the photon parameter, χ 0 = 0.5 (left panels), the pair production probability is significantly larger considering a LP background field instead of a CP one. For the higher value of χ 0 = 4 (right panels), CP increases, but only marginally, pair production. For this reason, we will from now on focus on LP background fields, unless otherwise specified.

Influence of the background field amplitude and gamma-photon energy
The relative influence on pair production of the photon energy γ γ and the background field amplitude a 0 (or equivalently intensity I 0 , as we consider here a given background field wavelength λ = 0.8 μm), is summarized in figure 3. It shows the probability P m for a high-energy photon to decay in a pair after crossing half a wavelength of a LP background field: (a) integrating numerically equation (17) which provides an exact value of P m ; (b) using the approximate but fully analytical expression given by Let us first discuss panel (a) and note that, in addition to the probability isocontours (solid white lines), the dotted white lines represent the contours of constant quantum parameter, χ 0 = 0.1, 1 and 16.5. As χ 0 ∝ γ γ a 0 , these are straight lines with a −45 • inclination.
In the limit χ 0 1, i.e. the bottom-left of figure 3, P m assumes very small values. In this limit, the time-integrated rate R m τ/4 scales as γ −1 γ χ (22)) and the exponential term, depending only on χ 0 , gives the dominating contribution to P m . As a result, the isocontours of P m (solid white lines) behave as nearly straight lines, roughly parallel to the contours of constant χ 0 . We wish to stress however that the dependence on γ −1 γ cannot be fully ignored: the isocontours of P m are less steep than the isocontours of χ 0 which indicates that P m increases faster with a 0 than with γ γ . In this range of small χ 0 , for which probability and time-integrated rate are equivalent, the importance of increasing a 0 to improve pair production was already pointed out in [36]. However, the probability variation with both a 0 and γ γ needs to be examined in more details at higher χ 0 .
Indeed, the dependence with γ −1 γ plays a more important role as χ 0 increases so that χ 0 cannot be considered as the only relevant parameter in order to optimize the pair production rate. This is clearly seen by considering the isocontours of P m in (a 0 , γ γ )-plane, figure 3(a): the isocontour slope becomes shallow and eventually changes sign. Thus, for any value of P m , a minimum field strength a 0 is needed in order to obtain the desired level of probability. A corollary is that, at constant a 0 , increasing γ γ increases the probability P m up to a maximum value, beyond which a further increase of γ γ would only decrease P m . The minimum of each isocontour can thus be found by solving ∂P m /∂γ γ = 0, which one can recast in the form Interestingly, this equation involves only χ 0 so that the minima of the probability isocontours lie on a straight line of constant χ 0 . This a priori surprising result can be better understood noting that both a 0 and P m are Lorentz invariants, the minimum value of a 0 to reach a given probability P m can depend only on the Lorentz invariant involving the photon energy, i.e. χ 0 . Solving numerically equation (26) for ε = 0 (LP case), we obtain 8 χ 0 16.5, which is reported by a dotted line in figure 3.
To conclude with figure 3, let us now turn to panels (b) and (c). Panel (b) reports the probability P m computed from equation (20) and the corresponding isocontours (red dashed lines). Superimposing the isocontours (white lines) of the probability from panel (a), we can see an excellent agreement between the two approaches. This validates the approximate form, equation (20), which has the advantage to be completely analytical and can now be used to get quick yet precise estimates of the probability P m . Last, panel (c) reports the probability extracted from 1D PIC simulations of the head-on collision of a flash of high-energy photons with a laser beam. The probability was extracted from the depletion of high-energy photons after crossing half a laser wavelength. Here again, a very good agreement is observed with the isocontours (white lines) obtained from the exact integration shown in panel (a). This last panel thus provides a cross-benchmark of our model and the Monte-Carlo module for nonlinear Breit-Wheeler pair production implemented in our PIC code SMILEI.
Our model can be used to predict the effectiveness of pair production for a given set of laser and high-energy photon parameters 9 . Most of the photons crossing the laser will be converted in pairs over a time Δt if the quantity RΔt in equation (25) becomes of order 1, for which 10 P tot (Δt) 0.63. In figure 4, we highlight in yellow the region in the (a 0 , γ γ )-plane where the probability for a high-energy photon to decay into a pair after interacting with a single maximum (half-wavelength) of the laser pulse P tot (τ/4) is equal or larger than 0.63 (the pulse duration is here τ /2, but the time required for the photon to explore half the laser wavelength is Δt = τ /4). The condition P tot (Δt) = 0.63 is also shown for pulses with a step-like time profile (solid lines) and duration 5τ (green, Δt = 2.5τ ) or 100τ (blue, Δt = 50τ ), and pulses with a sin 2 intensity time profile with FWHM 5τ (green dashed line, Δt = 5τ ) and 100τ (blue dashed line, Δt = 100τ ). Note that both step-like and sin 2 time profiles correspond to laser pulses with the same energy and peak power, and lead to similar requirements to reach P tot 0.63. Comparing the blue and green lines, we see that at constant γ γ a longer pulse satisfies the condition of efficient conversion for a lower value of a 0 . Hence, even though the condition Rτ /4 = 1 (yellow solid line) is quite stringent and achieving a high probability level over a time interval τ /4 can be difficult, the condition RΔt 1 is significantly relaxed considering longer interaction time. Indeed, significant (order 1) pair production probability is expected on forthcoming multi-PW laser facilities (see section 4 for details).
Let us also note that, when Rτ /4 1 (yellow area), more than 63% of conversion is achieved after crossing a single wavelength of the background field, and more than 99% of the incident gamma photons are converted into pairs in less than 3 periods. This gives a very robust condition for abundant pair creation, and it determines a limit above which is not useful to further increase the laser amplitude to increase the number of primary pairs. As we will show in the following, this condition (R τ /4 1) can also be invoked to find the optimal laser transverse shape and focal spot for a given laser energy. Indeed 3D PIC simulations, shown by the symbols in figure 4 and discussed in detail in the next section 3, show a different behaviour depending on whether their initial parameters lie in the yellow or green region of figure 4. In the case of the two simulations reported as squares in the figure, we will see that these differences are observed even though the simulation parameters correspond to the same value of χ 0 .
As a final remark, we point out that the probabilities considered so far are computed considering a given high-energy photon energy. In practice, we expect the high-energy photon beam to have a broad energy spectrum or at least a finite energy width (see e.g. references [25,32,[39][40][41][42][43][44][45][46]). Accounting for such energy distributions is however straightforward as there is no coherence effects between the high-energy photons. The total probabilities can thus be computed independently for all relevant high-energy photon energies so that P tot (γ γ ) is actually a function of γ γ . As an example, to compute the total number N pair of produced pairs, one only needs to integrate over the whole high-energy photon spectrum: N pair = dγ γ P tot (γ γ ) dN γ /dγ γ , where dN γ /dγ γ denotes the high-energy photon energy distribution.

Geometrical effects on pair production: the case of LG beams
This section aims to generalize the model developed in section 2.2.2, to overcome the plane wave approximation and account for the laser beam transverse profile. We start by reminding some properties of the LG beams, limiting ourselves to the critical aspects for the understanding of pair production. We then compare the predictions of our theoretical model (generalized to account for the laser transverse profile in section 3.3) with the results of 3D PIC simulations, whose initialization is discussed in section 3.2. A thorough comparison of the efficiency of pair production considering beams with the same energy and different field structures is presented in sections 3.4 and 3.5. We explore first the impact of the beam geometry and amplitude in a regime of efficient pair production (time-integrated rate R m τ /4 > 1) and then in a regime (R m τ /4 < 1) more relevant for current experimental facilities (as will be discussed in the following section 4). Finally in section 3.6, we discuss some properties of the produced pairs, such as the energy spectrum and the propagation direction. 9 In the case of large secondary pair production, our prediction will underestimate the number of pairs as our analysis is limited to the soft-shower regime. 10 Note that in this limit, the probability cannot be assimilated to the time-integrated rate, as in the case of very weak pair production.

Properties of LG beams
The LG modes are solutions of the paraxial equation with cylindrical symmetry and a generalization of the Gaussian beam. Each mode is indexed by two integers: p 0 for the radial order and for the azimuthal one. Hence, the notation LG p will be used in the following to identify a specific mode.
Starting from equation (10), the electric field of an LG beam can be expressed using E 0 = a 0 R u p , where the complex envelope u p in cylindrical coordinates is [64] where C p = p!/(p + | |)! is a normalization constant, L p (x) are the generalized Laguerre polynomials, w(z) = w 0 1 − z 2 /z 2 R is the beam waist with w 0 = w(z = 0) the waist at focus and z R = πw 2 0 /λ the Rayleigh length, and ψ p (z) = 2p + | | + 1 arctan z/z R is the generalized Gouy phase. Note that the LG modes described by equation (27) have all the same energy and that LG 00 corresponds to a Gaussian beam with maximum field amplitude a 0 .
In order to infer how the structure of the LG modes might affect pair production we show in figures 5(a)-(d) the quantum photon parameter χ γ for the head-on collision of photons of energy γ γ = 400 with a LP LG beam at focus (i.e. z = 0) with a 0 = 2000, p = 0 and ∈ [0, 3]. As can be seen in figures 5(a)-(d), for = 0 the maximum of the field, and thus the peak values of the quantum parameter, are no longer along the beam axis (i.e. x = y = 0), in contrast with the Gaussian beam case. The radius at which the field has its maximum amplitude ρ max can be obtained by solving ∂ ρ |u p | 2 = 0, leading to: with δ ij the Kronecker delta. For p = 0, equation (28) can be solved analytically and gives ρ max = w(z) | |/2. While the radial position of the field maximum for p = 0 increases with | |, the amplitude decreases, as shown in figure 5(e), and tends asymptotically to zero for | | 1 as |2π | −1/4 . The maximum amplitude of the envelope u 0 can be computed inserting ρ max into equation (27), which gives Figure 5(e) also shows u max p for p = 0, obtained by inserting the numerical solution of equation (28) into equation (27). The maximum field amplitude for the LG beams decreases with both | | and p, and as a consequence the same behaviour is observed for the maximum value of χ γ (see figures 5(a)-(d)). Therefore, in the following we will consider only | | 5 and p = 0, because of the higher field amplitude (hence higher values of χ γ ) and because of the practical challenges in producing high order modes at relativistic intensity. Note that in the considered cases with = 0, χ γ has 2| | maxima in the transverse plane (see figures 5(b)-(d)) and, looking at the evolution in time, each maximum makes a full rotation around the laser propagation axis in | |τ . This means that, in a fixed transverse plane, a maximum rotates to the location of the next one in half a period.
If we consider the head-on collision of the LG beam with a gamma photon, the photon is crossing half a period of the field in τ/4. Hence, the photon crosses a rotating local maximum in τ/4. This observation is fundamental for generalization of the reduced model presented in the previous section 2.2.2, where we consider subsequent intervals of duration τ /4 to infer the total probability.
As noted in figure 5, the LG beams not only have a non-trivial spatial distribution of the fields in the transverse plane, but also a maximum field amplitude decreasing with | | (following equation (29)). In order to distinguish the impact on pair production of these two aspects, we discuss the interaction with extended Gaussian (EG) beams. For each LG beam (i.e. each value of | |), we define the associated EG beam as the Gaussian beam that has the same maximum field amplitude and carries the same total energy. As u max p decreases with | | (equation (29)) and the energy is kept constant, the waist of the EG beams increases with respect to the fundamental Gaussian (equivalent to LG 00 ) as w = w 0 | |! e | |/2 | | −| |/2 w 0 .

Simulation set-up
In what follows, we present 3D PIC simulations performed with the open-source code SMILEI [54], considering the setup schematically represented in figure 6. Each simulation reproduces a volume of 40λ × 40λ × 14λ (in the x, y and z directions, respectively) with spatial resolution λ/24 and temporal resolution at 95% of the CFL condition [65].
An intense laser pulse (with wavelength λ = 0.8 μm) is injected from the left boundary z = −5λ using the method presented in reference [66] that allows for the Maxwell-consistent injection of an electromagnetic wave with an arbitrary spatio-temporal profile. Here, this method is used to inject the laser pulse by prescribing its spatio-temporal profile in its focal plane [i.e. the (x, y)-plane at z = 0], as given by the real part of equation (27) at z = 0, multiplied by a temporal envelope so that the laser pulse has a sin 2 temporal profile in intensity with FWHM τ FWHM = 5τ , and total duration of 10τ . The laser beam collides head-on with a counter-propagating flash of high-energy photons with a finite longitudinal extension L γ = λ/6, transverse extension equal to the simulation box, and a density n γ equal to the critical density n c = 0 m e ω 2 /e 2 . The duration of the simulation is long enough for the two light beams to stop interacting, but short enough so that none of the photons or secondary particles escape from the simulation box. All results presented hereafter have been extracted at the end of the simulations.
In sections 3.4 and 3.5, we present two series of simulations: the first one with a 0 = 2000 and γ γ = 400, the second with a 0 = 400 and γ γ = 2000, both series corresponding to a reference quantum parameter χ 0 6 × 10 −6 a 0 γ γ ∼ 4.85 (for λ = 0.8 μm). For each series, we varied the transverse spatial profile of the high-intensity laser pulse in order to study the impact of the spatio-temporal field profile, comparing different LG beams (varying , but keeping p = 0) and the corresponding EG beams. Note however that within each series, a 0 is left unchanged and refers to the maximum field amplitude of the reference Gaussian beam (for which we consider w 0 = 3λ). As a result, within each set of simulations, the laser pulse energy is unchanged, equal to 10 kJ for a 0 = 2000 and 410 J for a 0 = 400. Conversely, the photon flash surface energy density m e c 2 γ γ n γ d is 5 mJ/λ 2 for the first series (with γ γ = 400) and 25 mJ/λ 2 for the second one (γ γ = 2000). As highlighted in figure 4 (squares), the reference Gaussian case of the first series of simulations allows us to explore the physics in the regime of high pair production probability after a single τ/4 interval (i.e. R m τ /4 > 1, yellow region). On the contrary, with the second series we investigate a regime where efficient pair production is achieved thanks to the cumulative effects of several wavelength (i.e. R m τ /4 < 1, but the reference Gaussian case is still above the green line in figure 4).

Reduced model for pair creation with arbitrary spatio-temporal laser beam profiles
In this section we generalize the model developed in section 2.2.2 to accurately describe pair production from complex, spatially structured laser beams.
To do so, we exploit the fact that we are here dealing with high-energy photons colliding head-on with the strong background field. Hence, the photon momenta are left unchanged as the photon crosses the background field, and their trajectories are straight lines so that the field seen by a given photon along its trajectory is known and depends only on the photon initial position (x 0 , y 0 ) in the plane transverse to the direction of propagation of both light beams 11 . It follows that the quantum parameter of any photon is known throughout the photon interaction with the background field, and for a given background field, is determined uniquely by the photon initial position (x 0 , y 0 ). As a result, the total probability P tot for any photon to decay into a pair after interacting with the laser pulse is also determined by the photon initial position, and the total number of produced pairs at the end of the interaction between the flash of high-energy photons (with density n γ and longitudinal width L γ ) simply reads N pair = n γ L γ σ tot with σ tot = dx 0 dy 0 P tot (x 0 , y 0 ). (30) A key quantity to model pair production is thus the total cross-section σ tot , here defined as the integral of the total probability over the transverse plane. Notice that, as discussed in the last paragraph of section 2, for a given gamma photon spectrum σ tot can be considered a function of γ γ , and the total number of produced pairs can be computed integrating over the incident high-energy photon spectrum. Moreover, equation (30) can be easily generalised to account for a transverse profile of the incident high-energy photon beam. In particular, in the limit where this beam is much narrower than the laser beam, equation (30) reduces to σ tot σ γ P tot , where σ γ is the transverse area of the high-energy photon beam, and P tot is computed at the (transverse) location of the high-energy photon beam.
To compute the total probability at any position (x 0 , y 0 ) we use equation (25) and proceed as described in section 2.2.2, reconstructing the total probability from successive time intervals of duration τ /4 (as discussed in section 3.1, a photon at any position in the transverse plane explores in τ/4 a local maximum of the field amplitude whatever value assumes). One then only needs to compute the maximum field amplitude satisfactorily. For simplicity, we do so measuring the field amplitude at focus, taking the absolute value of the field envelope from equation (27) at z = 0. This approach, which neglects the effects of diffraction, is justified whenever the duration of the interaction (∼ τ FWHM ) does not exceed the time z R /c To test the impact of the approximations used in the model (phase and diffraction effects are neglected), we show in figure 7(b) the probability map extracted from a 3D PIC simulation, where all effects neglected in our model are taken into account. Here the probability is extracted at the end of the interaction as the ratio of the remaining photon linear density by the initial one. An excellent agreement is found with the theoretical prediction of panel (b), with a measured total cross-section σ tot 91.0 λ 2 , which confirms the predictive capability of our simple model for the case of spatially structured beams.

Simulations in the regime of high pair production probability (R m τ /4 > 1)
In order to investigate the dependence on the laser structure of the number of produced pairs, we present here a study in the regime satisfying R m τ /4 > 1 for the reference Gaussian beam, in which we expect efficient pair production. We remind that we consider here photons with normalized energy γ γ = 400 and a reference Gaussian beam with maximum amplitude a 0 = 2000, leading to χ 0 4.85 (all other parameters have been specified in section 3.2).
To test our analytical predictions and the assumptions at the base of our model, we performed simulations (hereafter referred to as frozen cases) in which the produced pairs are prevented to further radiate. In this condition, no secondary pairs are produced and the simulations allow to directly test the proposed analytical model. We also performed simulations including the full system dynamics (hereafter referred to as full cases), meaning that the produced pairs interact with the laser fields and potentially radiate. This allowed us to verify that we are still in the soft-shower regime, and that only a small departure from predicted number of produced pairs is indeed observed. Figure 8(a) shows the number of produced pairs at the end of the interaction for all the tested cases, LP LG (in black), CP LG (blue) and EG (green), dots refer to the frozen simulations and squares to the full ones. Considering a LP LG laser beam, pair production efficiency increases with . This means that for the high value of a 0 considered here, it is preferable to increase , even if it is decreasing the maximum amplitude (as shown in figure 5(e)), as it increases the characteristic transverse size of the beam, and so the interaction area where the probability is high (see inset in figure 8(a)). The EG beams (which have by construction the same maximum amplitude of the LG beams but Gaussian profile) are more efficient than the corresponding LG, up to = 4, as they have a bigger region with substantial probability in their probability maps. However, a slightly greater number of pairs are produced in the simulation with LG 05 than EG 05 beam. This because the maximum field amplitude has dropped substantially and a large part of the interaction region in the EG case has low probability of pair production (see inset in figure 8(a)). Predictions form our model discussed in section 3.3 and based on equation (30), shown in dashed lines, are in very good agreement with both EG and LG simulations.
Based on the prediction from figure 1, for the considered χ 0 4.85 we would expect at first look a higher efficiency in the CP case than with LP, contrary to what is obtained from PIC simulations (blue symbols in figure 8(a)). However, given the spatial structure of the laser field in the transverse plane, in a large part of the interaction region the quantum number is χ < 2.5, i.e. the value above which CP should be favoured based on the analysis shown in figure 1(b) obtained within the plane wave approximation. Hence, the slightly higher number of pairs produced in the CP case with respect to LP in the region where χ > 2.5, cannot compensate the higher efficiency of LP in the region of χ < 2.5. Note that the results using a CP beam are independent from the sign of (i.e. from the total angular momentum). This suggests that the total angular momentum of the laser has no effect on the number of pair produced. Note however that we do not account for the spin polarization effects. Since we are focusing on the soft shower regime where there are few secondary pairs, the dynamics of the pairs affected by spin effects is expected to play a minor role for pair production in our case. Moreover the dynamics of the produced pairs is not the main concern of this study. Nevertheless, in other configurations spin effects can impact the created pairs and gamma photons dynamics as shown in references [67,68].
For completeness, we compare the results of the frozen cases with simulations reproducing the full dynamics of the system (squares in figure 8(a)). As highlighted in figure 8(b), the difference in the produced number of pairs is always below 25% of its value in the frozen case, and it decreases with | |. This confirms that the majority of pairs comes from the conversion of primary gamma photons and not from subsequent radiation, as expected for a soft-shower. In this regime, our model, which accurately captures the physics of pair production for the frozen simulations (dashed lines), can be therefore used to predict pair production in realistic conditions.

Simulations in the regime of moderate pair production probability (R m τ /4 < 1)
In this second series of simulations we explore the regime of R m τ /4 < 1, maintaining the same maximum χ 0 as in the previous section by exchanging the values of field amplitude (now equal to a 0 = 400) and photons energy (now equal to γ γ = 2000). Note that considering a laser duration (τ FWHM = 5τ ), we still expect efficient pair production (the reference plane wave simulation lies above the green dashed line in figure 4). This can be seen in the inset of figure 9(a) for the Gaussian case LG 00 where a significant region exhibits a probability of order 1.
As figure 9 shows, for the simulations with LG beams the number of produced pairs is optimized for = 3. The increase for | | 3 can be explained as in the previous section by geometrical considerations, i.e. the area with the field being large enough to reach P tot of order one is getting wider with | |. In the contrary for | | > 3, the laser field amplitude becomes too low and having a laser beam with a larger transverse extension does not improve pair production. This is due to the fact that the field is becoming too low, so that the effect of decreasing P tot on pair production is more important, and is not compensated anymore by the increase in transverse size.
Moreover, in this regime, all LG beams perform better than the corresponding EG, given that the region with high (even though smaller than in the previous section 3.4) probability is larger for the LG than EG, as shown in the inset of figure 9 for = 5. Also, a similar behaviour as discussed in the previous section 3.4 is observed for the CP simulations and can be explained with the same reasoning. Last but not least, our model is reproducing with high accuracy the observed number of pair in this case too.

Energy and angular distribution of the produced pairs
The energy-angular distribution of the produced pairs in the z-x plane (i.e. the plane formed by the laser propagation directionẑ and the polarization directionx for the LP case) recorded at the end of the interaction is shown in figure 10 considering the reference LP and CP Gaussian beams. The top row gives the positron distribution in θ = arctan(p x /p z ) and in energy for the LP beams and the bottom one for the CP beams. The insets show the energy spectra of particles moving in the positive (same as the laser, black line) and negative (same as the gamma flash, red line) z-directions.
At high a 0 (panels (a) and (c)) the produced pairs, initially generated with momentum in the negative z-direction, are predominantly moving in the same direction as the laser (corresponding to θ = 0 • ) within a cone of aperture 40 • , meaning that they have been slowed down and then pushed by the laser. They reach a maximum energy of 10 4 m e c 2 for the LP case and 8 × 10 3 m e c 2 for the CP beam, consistent with the decrease of the peak field amplitude in CP with respect to LP at constant energy.
The spatial distribution is drastically modified when considering a 0 = 400. Even if there still is a substantial amount of pairs propagating in the laser direction, the majority of them keep their original direction, along the initial gamma photons propagation axis and anti-parallel to the laser pulse one. For these pairs the maximum energy is independent from the polarization, as shown in figures 10(b)-(d), contrarily to the pairs that are slowed down and pushed by the laser (right side of the quadrants or black line in the insets).
This study, supported by complementary simulations (not shown), suggests that for a 0 γ γ most of the created particles will end up propagating along the laser direction, while for smaller values of the laser maximum amplitude a large fraction of pairs will keep their original direction. A complete characterisation of the pairs spectrum and directionality is beyond the scope of this work, and will be investigated in more details in the future, as it can give important information for the planning of upcoming experimental campaigns.

Discussion
The model developed in this work allows to predict the pair production capability of upcoming facilities. In figure 11, we discuss the pair production maximum probability and total cross-section as a function of the photon energy and laser amplitude. The considered parameter range is relevant to ultra-short Ti:sapphire (λ = 0.8 μm) lasers, where studies of high-field physics and QED are already planned [69][70][71][72]. Two types of facilities can be distinguished. First, the Apollon 12 , CoReLS 13 and ZEUS 14 facilities (shown in black in figure 11) are designed to deliver multiple light beams with duration in between 15 and 30 fs FWHM and peak power from 1 to 10 PW. Second, the FACET-II 15 and LUXE 16 experiments will couple conventional electron accelerators with high-intensity (100 TW-class) ultra-short laser pulses (in white in figure 11). In Note in addition that, as most of these facilities will operate in the χ 0 1 regime, we report for the photon energy that of the electron beams either expected from laser-driven wakefield acceleration (for Apollon, CoReLS and ZEUS) or emerging from the linear accelerators (for FACET-II and LUXE). In the case of a broad (e.g. synchrotron-like) energy distribution for the high-energy photon source, the expected number of produced pairs can be computed integrating the reported values of P tot or σ tot over the full high-energy photon spectrum. As detailed below, this integration can be carried out in two rough but straightforward ways depending on the kind of facility (multi-PW or laser combined with a conventional accelerator) that is considered.
Let us now focus on panel (a), which shows the maximum total decay probability (from equation (25)) for a photon with energy γ γ colliding head-on with a LP Gaussian laser pulse with maximum field strength a 0 and duration τ FWHM = 25 fs (typical for the facilities listed above). This panel suggests the strategy to optimize pair production depending on the type of facility.
First, despite the extremely high electron beam energies, facilities such as LUXE and FACET-II (in white in figure 11) offer the possibility to probe the regime of moderate quantum parameter (χ 0 1) for which pair production may be observed (P tot 10 −2 ) but will not be abundant. This limit might be overcome by increasing the laser power up to few 100's of TW (300 TW are envisioned at LUXE), which allows to enter the quantum regime (χ 0 > 1) and increase P tot to values exceeding 0.1. On these facilities however, because with focusing aperture f/3 (black lines) and f/1.5 (red lines). In this last panel, at a given power, the laser energy is the same for the different focusing aperture f/N. P tot will not assume large values, increasing a 0 (e.g. by focusing the laser pulse as much as technically possible) is one of the most promising path to achieve abundant pair production.
In contrast, provided that multi-GeV electron beams can be obtained from laser wakefield acceleration, abundant pair production (P tot 0.63) is expected on all multi-PW laser facilities with the focalisation technique already in place (typically an aperture in between f/3 and f/4 was considered in figure 11).
Similar conclusions can be drawn from panel (b) where we examine in more details the influence of the laser pulse spatial profile, by looking at the total cross-section 17 σ tot (as defined by equation (30)) normalized to w 2 0 , considering a Gaussian laser beam with the same temporal properties as in panel (a). For the FACET-II and LUXE facilities, the total cross-section assumes very small values but increases fast as a 0 is increased: e.g. for LUXE, σ tot increases by 3 orders of magnitude increasing a 0 by a factor 10. Clearly, operating at large field strength will be a bottleneck for achieving abundant pair production on these facilities.
In contrast, in the range of parameters covered by multi-PW facilities, the dependence of σ tot with a 0 is much weaker. As in addition, σ tot w 2 0 (consistent with P tot 1 at the center of the beam), it becomes interesting for these facilities to increase the laser transverse size rather than opt for tight focusing.
Let us then discuss in more details the impact of focusing for the three upgrades of the Apollon facility (at 1, 3 and 7.5 PW). In panel (c), we present the total cross section (in units of λ −2 ) as a function of γ γ , considering either the standard f/3 aperture (black lines) or the more challenging f/1.5 aperture (red lines). For each laser power, a tighter focusing reduces the beam waist but increases the maximum a 0 .
In all cases, the cross section is rapidly increasing with the photon energy for γ γ 2 × 10 3 ( 1 GeV). It is almost flat as γ γ reaches 5 × 10 3 ( 2.5 GeV), which suggests that it is not worth increasing the photon energy above this value in forthcoming experiments to optimize pair production in the soft-shower regime. 17 Since the Gaussian field profiles depend on the transverse coordinate through the ratio ρ/w 0 (equation (27)), σ tot /w 2 0 is independent of w 0 .
As shown in panel (c), for a given laser power the cross section for f/1.5 is larger than for f/3 only for small values of γ γ , while the opposite behaviour is observed at large γ γ (of the order of a few GeV). This means that such tight focus, which is technologically very challenging, is not necessary for pair creation in forthcoming experiments on Apollon.
To test this prediction, we have performed complementary 3D PIC simulations of the 3 PW case with f/3 and f/1.5 (a 0 = 138.8 and 276.7, respectively) interacting with a gamma flash with energy γ γ = 10 4 . These simulations give σ tot 26.8 λ 2 for f/3, which is indeed larger than the value σ tot 19.5 λ 2 obtained for the tight-focusing f/1.5 case. These simulation results are consistent with our reduced model predictions, but also noticeably above its predictions, which are σ tot 16.2 λ 2 for f/3 and σ tot 6.5 λ 2 for f/1.5 (as given by figure 11(c)).
One of the reason of the discrepancy is that our calculations are performed at focus, neglecting diffraction. While this was a very good approximation for the parameters considered earlier (see section 3.3), it is not correct when considering tightly focused intense background fields and an interaction region larger than the Rayleigh length. The model can be easily generalized to include diffraction. Doing so improves our predictions to σ tot 22.2 λ 2 for f/3 and σ tot 9.8 λ 2 for f/1.5. These values are now significantly closer to our 3D PIC simulation results, and are also in excellent agreement (within 1%) with additional 3D PIC simulations in which secondary pair production is not accounted (frozen case). Secondary pairs, responsible for a further increase of about 20% for the f/3 case and by a factor ×2 for the tightly-focused f/1.5 case, are not included in the model and are the source of the remaining discrepancy. This suggests that for the highest peak intensity considered here we are at the limit of the soft-shower regime and the model should be considered as providing a lower bound estimate for pair production.
Let us finally briefly comment on how the effect of a broad energy spectrum and transverse extension of the high-energy photon beam can be accounted for depending on the type of facility considered.
In the case of FACET-II and LUXE facilities, both P tot and σ tot are very steep functions of γ γ , e.g. increasing by several (more than 4) orders of magnitude for the FACET-II configuration as γ γ is increased from 5 × 10 3 to 3 × 10 4 . One can thus expect that, for these facilities, only the photons with the highest energy contribute to pair production. The total number of produced pairs can then be estimated considering only those photons with the highest energy (maximum γ γ as reported on panels (a) and (b)). If the high-energy photon beam is much narrower that the laser beam, the total number of produced pair can be computed from the number of photons and P tot evaluated at this maximum γ γ . If instead the high-energy photon beam has a large transverse distribution, the final number of pairs can be calculated by taking the high-energy photon surface density and σ tot computed at the maximum γ γ .
For the multi-PW laser facilities, and large background field strength a 0 , we have seen that both P tot and σ tot (for the latest, see e.g. panel (c) for the Apollon facility) reach a plateau above a minimum value of the photon energy (e.g. γ γ 5 × 10 3 in the case of panel (c) discussed above). Hence, pair production will follow primarily from the decay of the high-energy photons in this plateau. A simple estimate of the number of produced pairs will then be obtained from either the number of all photons in the plateau region and a typical value of P tot (if the photon beam transverse size is small with respect to the laser beam waist) or from their surface density and σ tot obtained for this energy range (if the photon beam transverse size is large).

Conclusion
In summary, we have presented a systematic study of nonlinear Breit-Wheeler pair production in the head-on collision of high-energy gamma photons with an ultra-high intensity laser beam. Combining analytical modeling and PIC simulations embarking the relevant QED modules, we have evaluated the impact on the efficiency of pair production of the laser spatio-temporal profile (comparing e.g. Gaussian with different focal spots and LG pulses) and parameters such as polarization, intensity and duration. Motivated by experimental constraints, we have considered fixed laser energies while changing these parameters.
We have explored laser field strength and photon energy parameter ranges relevant for currently and/or upcoming high-power laser facilities, focusing in particular on ultra-short Ti:sapphire laser facilities. A reduced model was proposed that allows to describe pair production in the regime of soft-shower, where secondary pair production can be neglected. This model, that would in principle allow to predict a minimum number of produced pairs, is found to agree remarkably well with 3D PIC simulations over a broad range of parameters, highlighting the importance of the soft-shower regime for experiments on the forthcoming laser facilities.
The model also allows to distinguish two regimes of interaction depending on whether the probability for a photon to decay into a pair as it crosses a single half-wavelength of the laser pulse is large (of order 1) or not. The two regimes are not however fully determined by the photon maximum quantum parameter χ 0 , so that the laser field strength parameter a 0 and photon energy γ γ do not play a symmetric role in the system dynamics. This was confirmed in 3D PIC simulations where the two regimes were investigated at χ 0 = 4.85 considering laser beams with complex spatio-temporal profiles, LG beams in particular, and different polarizations. It was found that, for a fixed laser energy, using circular polarization or LG does not improve significantly pair production. As they are also quite difficult to achieve experimentally, their effect on pair production is too marginal to be interesting for applications. It was found however that, in the regime of high pair production probability, it can be preferable not to focus much the laser beam (keeping the same total energy with lower peak intensity) and maintain a Gaussian shape in order to maximise the area with high-enough fields for efficient pair production.
These findings help to draw guidelines for future experiments on ultra-high intensity facilities. We show in particular that the path to abundant pair production is different whether one considers 100 TW-class or multi-PW laser systems. In particular, for facilities such as FACET-II or LUXE, 10 to 100 TW-class lasers will be coupled to ultra-relativistic 10-20 GeV electron beams emerging from conventional electron accelerators. On these facilities, pair production will develop in the low probability regime and operating in a tightly-focused configuration to access the highest possible laser field strength will be a major experimental challenge. In contrast, using multiple laser beam multi-PW facilities will operate in the regime of high pair production probability. In this regime, the pair production efficiency (measured in our work in terms of a total cross section) depends more strongly on the laser pulse transverse size than on its maximal field strength so that operating with standard (not too tight) focusing aperture increases the number of produced pairs. Interestingly, producing high-energy photons (or electrons, e.g. through laser wakefield acceleration) of a few (∼ 5) GeV shall be enough for these facilities as higher energies do no help increase significantly the pair production efficiency in this high probability regime.
To conclude, this work provides a deeper understanding of the optimal conditions for pair production in upcoming experimental campaigns exploring nonlinear Breit-Wheeler pair production. It focused on pair production seeded by high-energy photons, and it will be interesting for future works to investigate the differences with the case where ultra-relativistic electrons are used to seed the pair production process. Studying the transition from the soft-shower to the cascade/avalanche regime in which secondary pairs and radiations play a dominant role is another interesting perspective.
With the ansatz equation (A.1), the integral can be performed analytically, leading to I ε (χ 0 ) π b 0 (χ m ) F(s ε (χ m )) with F(s) = 2/πs erf π √ 2/(4s) . (A. 3) The ansatz equation (A.1) gives a good approximation of the integrand for arbitrary values of ε as long as χ m is small enough. It is also exact (for all χ 0 ) in the cases ε = ±1 for which s ±1 (χ m ) → +∞. Hence, equation (A.3) provides a very good, fully analytical approximation of equation (19) for a wide range of ε and χ 0 . The approximation however needs to be corrected for χ 0 1 and |ε| < 1. In this limit, the integrand in equation (19) can be approximated using equation (6) (20). This form ensures the correct asymptotic behaviour of the integral equation (19) in both small and large χ 0 limits. It departs from the exact expression (integrated numerically) by less than 20%, the error being maximum for intermediate values of