Dominance of γ-γ electron-positron pair creation in a plasma driven by high-intensity lasers

Creation of electrons and positrons from light alone is a basic prediction of quantum electrodynamics, but yet to be observed. Our simulations show that the required conditions are achievable using a high-intensity two-beam laser facility and an advanced target design. Dual laser irradiation of a structured target produces high-density γ rays that then create >108 positrons at intensities of 2× 1022 Wcm−2. The unique feature of this setup is that the pair creation is primarily driven by the linear Breit-Wheeler process (γγ → e+e−), which dominates over the nonlinear Breit-Wheeler and Bethe-Heitler processes. The favorable scaling with laser intensity of the linear process prompts reconsideration of its neglect in simulation studies and also permits positron jet formation at experimentally feasible intensities. Simulations show that the positrons, confined by a quasistatic plasma magnetic field, may be accelerated by the lasers to energies > 200 MeV.


Introduction
High-power lasers, focused close to the diffraction limit, create ultrastrong electromagnetic fields that can be harnessed to drive high fluxes of energetic particles and to study fundamental physical phenomena 1 . At intensities exceeding 10 23 Wcm −2 , those energetic particles can drive nonlinear quantum-electrodynamical (QED) processes 2, 3 otherwise only found in extreme astrophysical environments 4,5 . One such process is the creation of electron-positron pairs from light alone. Whereas multiphoton (nonlinear) pair creation has been measured once, using an intense laser 6 , the two-photon process (γγ → e + e − , referred to here as the linear Breit-Wheeler process 7 ) has yet to be observed in the laboratory with real photons. As the probability of the nonlinear process grows nonperturbatively with increasing field strength 8,9 , it is expected to provide the dominant contribution to pair cascades in high-field environments, including laser-matter interactions beyond the current intensity frontier 10,11 and pulsar magnetospheres 12 .
The small size of the linear Breit-Wheeler cross section means that high photon flux is necessary for its observation. Achieving the necessary flux requires specialized experimental configurations 13,14 and therefore its possible contribution to in situ electron-positron pair creation has hitherto been neglected in studies of high-intensity laser-matter interactions. However, these interactions create not only regions of ultrastrong electromagnetic field, but also high fluxes of accelerated particles, because relativistic effects mean that even a solid-density target can become transparent to intense laser light 15,16 . In the situation of multiple colliding laser pulses, which is the most advantageous geometry for driving nonlinear QED cascades 10,[17][18][19] , there are, as a consequence, dense, counterpropagating flashes of γ rays, and so the neglect of linear pair creation may not be appropriate.
Recent construction of multi-beam high-intensity laser facilities, such as ELI-Beamlines 20 , ELI-Nuclear Physics 21,22 , and Apollon 23 , and a significant progress in fabrication of µm-scale structured targets 24,25 open up qualitatively novel regimes of pair production for exploration. Specifically, we show that a structured plasma target irradiated by two laser beams creates an environment where the linear process dominates over the nonlinear and over the Bethe-Heitler process. Remarkably, this regime does not require laser intensities beyond than what is currently available. At I 0 < 5 × 10 22 Wcm −2 , the positron yield from the linear process is ∼ 10 9 , which is four orders of magnitude greater than that envisaged in Refs. [13] and [14]. These positrons are generated when two high-energy electron beams, accelerated by and copropagating with laser pulses that are guided along a plasma channel, collide head-on, emitting synchrotron photons that collide with each other and the respective oncoming laser. Not only does this provide an opportunity to study the linear Breit-Wheeler process itself, which is of interest because of its role in astrophysics [26][27][28] , but also the transition between linear and nonlinear-dominated pair cascades. In an astrophysical context, the balance between these two determines how a pulsar magnetosphere is filled with plasma; as in the laser-plasma scenario, the controlling factors are the field strength and photon flux [29][30][31][32] . We also show that the positrons, created inside the plasma channel coterminously with the laser pulses, may be confined and accelerated to energies of hundreds of MeV, which raises the possibility of generating positron jets. The transverse confinement needed to accelerate positrons is provided by a slowly evolving plasma magnetic field. Crucially, it is the same field that enables acceleration of the ultra-relativistic electrons prior to the collision of the two laser pulses.

Results
The target configuration considered in this work is shown in Fig. 1(a). A structured plastic target with a pre-filled channel is irradiated from both sides by two 50-fs, high-intensity laser pulses that have the same peak normalized laser amplitude a 0 , in the range 100 ≤ a 0 ≤ 190. Here a 0 = 0.85I , where I 0 is the peak intensity of the laser and λ 0 = 1 µm its wavelength in vacuum. The target structure, where a channel of width d ch = 5 µm and electron density n e = (a 0 /100)3.8n c is embedded in a bulk with higher density n e = 100n c , enables stable propagation 33 and alignment of the two lasers. Here n c = πmc 2 /(eλ 0 ) 2 is the so-called critical density, where e is the elementary charge, m is the electron mass, and c is the speed of light. At relativistic laser intensities (a 0 1), the cutoff density for the laser increases roughly linearly with a 0 due to relativistically induced transparency. Scaling the channel density with a 0 ensures that the optical properties of the channel and thus the phase velocity of the laser wave-fronts are approximately unchanged with increase of a 0 . Structured targets with empty channels have successfully been used in experiments 24,25 and it is now possible to fabricate targets with prefilled channels, similar to those considered in this work 1 .
The interaction is simulated in 2D with the fully relativistic particle-in-cell (PIC) code EPOCH 34 , which includes Monte Carlo modules for quantum synchrotron radiation and nonlinear pair creation 35 . At each time-step, the quantum synchrotron radiation module computes the quantum nonlinearity parameter, for each charged macro-particle using the electric and magnetic fields (E and B) at the particle location, as well as the particle relativistic factor γ and velocity v. Here E S ≈ 1.3 × 10 18 V/m is the Schwinger field [36][37][38] . The parameter χ controls the total radiation power and the energy spectrum of the emitted photons. In the quantum regime χ 1, which is reached in this work, it is necessary to take into account the recoil experienced by the particle when emitting individual photons. This is done self-consistently by the PIC simulation, which uses the Monte Carlo algorithm described in Refs. [35,39]. Note that, since the ion species is fully ionized carbon, Bethe-Heitler pair creation, already demonstrated in laser-driven experiments 40,41 , may be neglected. Detailed simulation and target parameters are provided in the Methods section. All the results presented in this paper have been appropriately normalized by taking the size of the ignored dimension to be equal to the channel width d ch , i.e. 5 µm.
The plasma channel, being relativistically transparent to the intense laser light 15,16 , acts as an optical waveguide. The laser pulses propagate with nearly constant transverse size through the channel, pushing plasma electrons forward. This longitudinal current generates a slowly evolving, azimuthal magnetic field with peak magnitude 0.6 MT (30% of the laser magnetic field strength) at a 0 = 190, as shown in Fig. 1(c). The magnetic field enables confinement and direct laser acceleration of the electrons 33,42 . After propagating for ∼30 µm along the channel, laser #2 in Fig. 1(a) has accelerated a left-moving, high-energy, high-charge electron beam that performs transverse oscillations of amplitude ∼2 µm: the number of electrons with relativistic factor γ > 800 is 4 × 10 11 , which is equivalent to a charge of 64 nC. Laser #1 generates a similar population of electrons moving to the right, with a representative electron trajectory shown in Fig. 3(c).
The plasma magnetic field has an essential role in enabling generation of ultrarelativistic electrons. Transverse deflections by the magnetic field keep p y antiparallel to the transverse electric field E y of the laser, despite the oscillation of the latter. As a result, the electron continues to gain energy while moving along the channel and performing transverse oscillations, as may be seen in Fig. 3(a) and Fig. 3(c). In the absence of the magnetic field, the oscillations of E y would terminate the energy gain prematurely. The magnetic field of the plasma has to be sufficiently strong to ensure that the electron deflections occur on the same time scale as the oscillations of E y . This criterion can be formulated in terms of the longitudinal plasma current and is provided in Ref. [42]. The time evolution of the electron spectrum, which shows this mechanism in action, is shown in Fig. 2(a).
The target length is such that no appreciable depletion of the laser pulses occurs by the time they reach the midplane (x = 0), t = 0. Here, the high-energy electron beams collide head-on with the respective oncoming laser pulse, each of which has an intensity at least as large as its initial value (the magnitude can increase slightly due to pulse shaping during propagation along the channel). This configuration maximizes the quantum nonlinearity parameter χ for the electrons, as the two terms under the square root in Eq. (1) are additive for counterpropagation. In copropagation, by contrast, they almost cancel each other. This is why radiation prior to the collision, when electrons propagate in the same direction as the accelerating laser pulse, is driven primarily by the plasma magnetic field. As is shown in Fig. 2(b), the largest value of χ ≈ 0.25 until the collision occurs; immediately thereafter, the cancellation is eliminated and χ increases rapidly to approximately 1.25. Figure 1(b) shows the impact of the collision on the energetic electrons from Fig. 1(a): they radiate away a substantial fraction of the energy they gained during the acceleration phase and are scattered out of the channel. Similar behavior is shown in Fig. 3(c): the electron encounters the counterpropagating laser beam at about t = 10 fs and then its energy decreases rapidly.
The observed increase in χ during the electron-laser collision increases the radiation power of the individual electrons. This emission occurs in a highly localized region, which leads to the marked increase in photon density shown in Fig. 4 for the case where a 0 = 190. The conversion efficiency of the laser energy into photons with energies 100 keV ≤ ε γ ≤ 10 MeV is shown in Fig. 5 over a wide range of a 0 . We are interested in the photons in this energy range because these are the photons that Original electron Generated positron 3 participate in the linear Breit-Wheeler process in our setup, as shown in Fig. 11 of the Methods section. As expected, there is a significant increase in the conversion rate caused by the electron-laser collision. The configuration under consideration here therefore represents a micron-scale, plasma-based realization of an all-optical laser-electron-beam collision (see also Ref. [17]). This geometry is the subject of theoretical [43][44][45] and experimental 46,47 investigation into radiative energy loss in the quantum regime, as well as nonlinear pair creation 48 . It is worth emphasizing that the use of the structured target has two key benefits compared to the commonly used gas targets: automatic alignment of the colliding electrons with an oncoming laser beam and a considerably higher density of colliding electrons. The angularly resolved spectrum of the emitted photons is shown in Fig. 6. There are approximately 2 × 10 14 photons with energies between 100 keV and 10 MeV and with 90 • ≤ θ ≤ 180 • . This is essentially half of the energetic photons emitted by the left-moving electrons (the other half is emitted with −180 • ≤ θ ≤ −90 • and has a similar spectrum). The photons emitted by one electron beam collide with both the oncoming laser and the photons emitted by the other electron beam. The former drives electron-positron pair creation by the nonlinear Breit-Wheeler process, γ EM field −−−−→ e + e −2, 9 : at a 0 = 190, our simulations predict a yield of 5 × 10 8 pairs.
The positrons subsequently undergo direct laser acceleration in much the way as the electrons: PIC simulations show

Initial channel boundary
Initial target boundary

Energy conversion efficiency [%]
! " Prior to the collision Entire interaction Figure 5. Conversion efficiency of the laser energy into γ rays with energies 100 keV ≤ ε γ ≤ 10 MeV, in 2D PIC simulations with two counterpropagating laser pulses that have the same peak normalized laser amplitude a 0 . The blue markers show the conversion efficiency recorded before the two pulses collide, whereas the red markers show the conversion efficiency computed over the duration of the entire laser-target interaction.
that the typical relativistic factor of a right-moving positron increases to γ ≈ 1000 as it propagates from x ≈ 0 to x ≈ 20 µm. This is illustrated in Fig. 1(c) and corroborated by the time evolution of the positron energy spectrum shown in Fig. 7(a). A representative trajectory for a positron moving from the central region towards the left target boundary is shown in Fig. 3(d). Acceleration is made possible by the plasma magnetic field, which is confining (on the left-hand side of the target) for electrons moving to the right, or equivalently, positrons moving to the left [compare Fig. 3(c) and Fig. 3(d)]. Crucially, Fig. 1(c) shows that this magnetic field polarity is preserved well after the lasers and electron beams collide. This is why, after the two laser pulses collide and pass through each other, they can accelerate the positrons, but not the electrons, created in by photon-photon collisions, as seen in Fig. 7. The generated electrons are not transversely confined in our magnetic field configuration when moving from the center towards either of the channel openings. However, the continued propagation of the lasers along the channel raises the possibility of accelerating positron jets, if there is sufficient pair creation in the channel center.
, of the photons emitted inside the channel in the 2D PIC simulation from Fig. 1, where θ is the angle defined in Fig. 1 We now show that this is the case, and furthermore that the pair creation is dominated by the linear Breit-Wheeler process. The cross section is 7 : where r e = e 2 /(mc 2 ) is the classical electron radius, β = 1 − 1/ς , and √ ς is the normalized center-of-mass energy, ς = ε 1 ε 2 (1 − cos ψ)/(2m 2 c 4 ), for two photons with energy ε 1,2 colliding at angle ψ. Equation (2) is the cross section for Laser collision Laser collision two-photon pair creation in vacuum: while it is modified by a strong EM field [49][50][51][52] , these corrections, which scale as (χ γ /ς ) 2 for photon quantum nonlinearity parameter χ γ 53 , are negligible for the scenario under consideration here (see the Supplemental Material for details).
We take as a representative value σ γγ ≈ 2r 2 e (approximately its maximum, at ς ≈ 2) and assume that we have two photon populations of number density n γ , colliding head-on in a volume of length cτ (the laser pulse length) and width d ch (the width of the channel). The number of photons (in each beam) is N γ ≈ 10 9 λ 0 [µm]P γ (n e /n c )(cτ/λ 0 )(d ch /λ 0 ) 2 , where P γ is the number of photons emitted per electron, n e is the electron number density and λ 0 is the laser wavelength. The number of positrons produced, The physical parameters are n e = 7n c , τ = 50 fs, d ch = 5 µm, and λ 0 = 1 µm. The number of photons emitted per laser period by a counterpropagating electron is P γ ≈ 18αa 0 , where α 1/137 is the fine-structure constant. By setting P γ = 20, we obtain a total number of photons, 2N γ ≈ 1.3 × 10 14 , which is approximately consistent with the simulation result. As a consequence, we predict that N BW lin ≈ 7 × 10 9 . Given that P γ ∝ a 0 and n e ∝ a 0 , we predict a scaling of N BW lin ∝ a 4 0 .

0.07
Laser #1 Laser #2 Initial channel boundary This is considerably larger than the number of pairs expected from the nonlinear Breit-Wheeler process; moreover, as the probability rate for the latter is exponentially suppressed with decreasing a 0 , we expect the yield to be much more sensitive to reductions in laser intensity. The potential dominance of the linear process motivates a precise computation, which takes into the account the energy, angle and temporal dependence of the photon emission.
However, direct implementation of the linear Breit-Wheeler process in a PIC code is a significant computational challenge, as it involves binary collisions of macroparticles and the interaction must be simulated in at least 2D. The simulation at a 0 = 190 6/16 generates ∼10 8 macrophotons in the energy range relevant for linear Breit-Wheeler pair creation and therefore ∼10 16 possible pairings. This can be reduced by using bounding volume hierarchies 54 , which is effective if the photon emission and the pair creation are well-separated in time and space. In our case, there is no such separation. As such, we postprocess the simulation output to obtain the yield of linear Breit-Wheeler pairs, using the algorithm described in Methods. Note that the photons used to compute this yield are the same photons used by the simulation to compute the yield of nonlinear Breit-Wheeler pairs. As such, while the photon number would change if the simulation were performed in 3D rather than 2D, the yield of both processes would be affected in a similar way. The location and time that pairs are created by the linear process, as determined by this algorithm for the case that a 0 = 190, are shown in Fig. 8 and Fig. 9 respectively. Approximately 59% of the pairs are created inside the original channel boundary. The majority (74%) of pairs are created by photons emitted after t = 0, when the high-energy electrons collide with the respective counterpropagating laser. There is a smaller contribution from photons that are emitted during the acceleration phase, t emit < 0; radiation in this case is driven by the plasma magnetic field, because the energetic electrons are moving in the same direction as the laser 33,55 . The dominance of the post-collision contribution is caused by the increase in the quantum parameter χ e for counterpropagation. The fact the pair creation overlaps with the laser pulses (in both time and space) indicates that the positrons could be accelerated out of the channel, as the magnetic field, shown in Fig. 1(c), has the correct orientation to confine them.
The pair yields for the linear and nonlinear processes are compared in Fig. 10. The results for the latter are obtained by performing four simulation runs for each value of a 0 with different random seeds: points and error bars give the mean and standard deviation obtained, respectively. At a 0 < 145, fewer than ten macropositrons are generated per run, so the corresponding data points are not shown. Our analytical estimates for linear Breit-Wheeler pair creation lead us to expect a yield that scales as a 4 0 : this is consistent with a power-law fit to the data in Fig. 10, which gives a scaling ∝ a m 0 , where m ≈ 3.93. We find that the linear pair yield is significantly larger for a 0 < 190.
The number of positrons produced by the linear Breit-Wheeler process exceeds 10 6 even for a 0 = 50, equivalent to I 0 = 3.4 × 10 21 Wcm −2 , which is well in reach of today's high-power laser facilities. In order to determine whether this is sufficient to be observed, we estimate the number of pairs produced by the Bethe-Heitler process, which is the principal source of background. In this process, a γ ray with energyhω > 2mc 2 creates an electron-positron pair by interacting with the Coulomb field of an atomic nucleus. The calculation is described in detail in the Methods section. We sum the pair creation probabilities for each simulated photon, taking into account the distance each photon travels in the plasma channel, to obtain the blue circles in Fig. 10. The Bethe-Heitler background is smaller than the linear Breit-Wheeler signal by approximately two orders of magnitude, which supports the feasibility of using a plasma channel as a platform for investigating fundamental QED effects.

Discussion
We have shown that laser-plasma interactions provide a platform to generate and accelerate positrons, created entirely by light and light, at intensities that are within the reach of current high-power laser facilities. While previous research into pair creation at high intensity has focused largely on the nonlinear Breit-Wheeler process, we show that the high density of photons afforded by a laser-plasma interaction can make the linear process dominant instead. As such, the geometry we consider has the potential to enable the first experimental measurement of two-photon pair creation, driven entirely by real photons. More broadly, it motivates reconsideration of the neglect of two-particle interactions in simulations of dense, laser-irradiated plasmas. Such interactions will form a major component of the physics investigated in upcoming high-power laser facilities. From the theory perspective, our results also motivate investigation of field-driven corrections to the two-photon cross section. The theory for the inverse process, pair annihiliation to two photons, has recently been revisited 56 .
One of our surprising findings, besides the dominance of the linear Breit-Wheeler process, is that the plasma magnetic field preserves its polarity after the two laser pulses collide and pass through each other. The polarity of the magnetic field enables transverse confinement of the positrons within the channel and their acceleration by one of the laser pulses to energies approaching 1 GeV. We have confirmed this directly for the positrons generated via the nonlinear Breit-Wheeler process. This should also be the case for the positrons generated via the dominant linear Breit-Wheeler process, because the particles are created inside the channel magnetic field in the presence of a laser pulse, which are the prerequisites for the direct laser acceleration. We therefore expect the positrons to be ejected from the target in the form of collimated jets. The collimation should aid positron detection outside of the target. Moreover, their detection at lower values of a 0 should be a clear indicator of the linear Breit-Wheeler process being the source, as the nonlinear process is heavily suppressed for a 0 150.
Finally, we point out that our observations regarding the dominance of the linear Breit-Wheeler process apply to a range of channel densities. In our simulations, the electron density in the channel is set at n ch = (a 0 /100)3.8n c , such that it increases linearly with a 0 during the intensity scan. Two channel density scans provided in the Supplemental Material show that our observations hold for channel densities that are within a ±20% window of n ch .  Table 1 provides detailed parameters for the simulations presented in the manuscript. Simulations were carried out using the fully relativistic particle-in-cell code EPOCH 34 . All our simulations are 2D-3V. The axis of the structured target is aligned with the axis of the counterpropagating lasers (laser #1 and laser #2) at y = 0. The target is initialized as a fully-ionized plasma with carbon ions. The bulk electron density is constant during the intensity scan while the electron density in the channel is set at n e = (a 0 /100)3.8n c . Each laser is focused at the corresponding channel opening. The lasers are linearly polarized with the electric field being in the plane of the simulation. In the absence of the target, the lasers have the same Gaussian profile in the focal spot with the same Gaussian temporal profile.

Particle-in-cell simulations
We performed additional runs at a 0 = 190 with higher spatial resolutions (40 by 40 cells per µm and 80 by 80 cells per µm). There are no significant variations in the photon spectra for multi-MeV photons and for photons with energies above 50 keV. The electrons that emit energetic photons, as the one whose trajectories in physical and momentum space are shown in Fig. 3(a) and Fig. 3(c), undergo their energy gain without alternating deceleration to non-relativistic energies and re-acceleration. This is likely the reason why they are not subject to a more severe constraint discussed in Refs. [57][58][59] that requires for the cell-size/time-step to be reduced according to the 1/a 0 scaling in order to achieve convergence.

Postprocessing algorithm for determination of the linear Breit-Wheeler pair yield
In order to compute the yield and spatial distribution of the linear Breit-Wheeler pairs, we approximate the photon population as a collection of collimated, monoenergetic beamlets. Discretization into beamlets is achieved by recording the location (x 0 , y 0 ), energy ε γ , and angle θ of each photon macroparticle at the time of the emission t emit . The photon emission pattern suggests that the emission profile across the channel can be approximated as uniform. We thus represent the emitted photons by a time-dependent distribution function f = f (x 0 , s γ , θ ;t emit ), where s γ ≡ log 10 (ε γ /MeV). It is sufficient to limit our analysis to −40 µm ≤ x 0 ≤ 40 µm, −3 ≤ s γ < 3, and 0 • ≤ θ ≤ 180 • . We split each interval into 70 equal segments to obtain 2.6 × 10 5 beamlets. We only check for collisions of beamlets propagating to the right with beamlets propagating to the left. The yield is multiplied by a factor of two to account for beamlets with −180 • ≤ θ ≤ 0 • .
The temporal dependence of a beamlet is represented by slices of given density and fixed thickness. For each beamlet pairing, our algorithm finds the interaction volume V , the intersections of the beamlet axes and the crossing angle ψ. The pair yield is given by ∆N BW lin = σ γγ c(1 − cos ψ)V n 1 n 2 dt, where n 1 and n 2 are the photon densities in two overlapping slices at the intersection point. In general, the shape of the overlapping region is not rectangular, so the pair creation is visualized by depositing ∆N BW lin onto a rectangular grid, into cells with centers inside volume V . The procedure is repeated for each beamlet pairing to obtain the density of generated pairs.
To show that the limitation −3 ≤ s γ < 3 is justified, we plot the distribution of linear Breit-Wheeler pairs as a function of photon energies. We use s γ rather than ε γ to capture a wide range of energies. Figure 11 confirms that the pair yield drops off for |s γ | > 2, which justifies the energy range selected in the manuscript. The algorithm is a simplification that replaces a direct approach of evaluating all possible collisions of beamlet slices. In a head-on collision, each slice collides with many counterpropagating slices within the interaction volume, which makes the calculation computationally intensive. Our algorithm takes advantage of the fact that the typical duration of beamlet emission, τ, is much longer than the time it takes for photons to travel between the sources emitting the two beamlets, /c, where is the distance between the sources in the case of near head-on collision. As shown in the Supplemental Material, our approach is a good approximation as long as /cτ < 1, with the error scaling as ( /cτ) 2 .
The postprocessing algorithm neglects the depletion of the photon population due to the linear Breit-Wheeler process. This is justifiable, because only a small fraction of the considered photons actually pair-create (and would therefore be lost). Using the maximum photon density of n γ ≈ 600n c from Fig. 4, we obtain a mean free path with respect to the linear Breit-Wheeler process, that is much larger than the characteristic size of the photon cloud of 10 µm. We estimate depletion of the photon population due to the linear Breit-Wheeler process (as a fraction of the initial size) to be smaller than 2 × 10 −4 .

Estimated background from Bethe-Heitler pair creation
The principal source of background in a prospective measurement of linear (or nonlinear) Breit-Wheeler pair creation is the Bethe-Heitler process, wherein a photon with energyhω > 2mc 2 produces an electron-positron pair on collision with an atomic nucleus 60 . In order to estimate the contribution from this process, we sum the pair creation probability for each macrophoton in the simulation: N BH = ∑ k w k P k,BH , where w k is the weight of the kth photon, scaled assuming that the third dimension has size 5 µm. The probability P k,BH = n i k σ BH , where n i is the density of carbon ions in the channel and k is the distance the photon travels before it leaves the channel. We estimate n i and k using the unperturbed properties of the channel, i.e. those at the start of the simulation, and taking into account the photon's point of emission and direction of propagation. Thus n i = 3.8a 0 n c /(100Z). We approximate the cross section σ BH by that for an unscreened, fully ionized, point carbon nucleus (formula 3D-0000 in Ref. [61] with Z = 6). The functional dependence of the cross section on the normalized photon energy γ =hω/(mc 2 ) is given by σ BH (γ) αr 2 e Z 2 (2π/3)[(γ − 2)/γ] 3 for γ − 2 1 and σ BH (γ) αr 2 e Z 2 [28 ln(2γ)/9 − 218/27] for γ 1, where r e is the classical electron radius 61 . Our results are shown as blue circles in Fig. 10. This estimate neglects contributions from pair creation in the plasma bulk, which can be controlled by reducing the thickness of the channel walls. Furthermore, the difference in magnitude between background and signal is sufficiently large that it provides a margin of safety.

Supplemental material
Channel density scan Table 2. Number of pairs at a 0 = 160 for different electron densities, n e , in the channel. The density n e is given in terms of n ch = (a 0 /100)3.8n c . n e /n ch 0.8 In our simulations, the electron density in the channel is set at n ch = (a 0 /100)3.8n c . This value was chosen to achieve prolonged laser propagation inside the channel, such that each laser pulse can generate ultrarelativistic electrons without becoming significantly depleted prior to the collision. In order to show that the observed trend is valid for a range of channel densities, we have performed channel density scans for a 0 = 160 and a 0 = 190. The pair yield for the linear and nonlinear Breit-Wheeler processes is given in Table 2 for a 0 = 160 and in Table 3 for a 0 = 190. The ratio of the linear to nonlinear pair yield increases for 0.8 ≤ n e /n ch ≤ 1.2 as we reduce a 0 from 190 to 160. This is the trend that is reported in the main text, which indicates that our observations apply to a range of channel densities and no fine-tuning is necessary.  Figure 12 provides additional information on the production of nonlinear Breit-Wheeler pairs. The horizontal scale shows the emission time, t emit , of energetic photons that go on to produce pairs by interacting with the laser photons. The energetic photons are generated and subsequently propagated as particles by the PIC code. The time of the pair production, t BW nonlin , is shown along the vertical scale. The number of pairs for each pairing of t emit and t BW nonlin is color-coded. The pair production is directly computed by the PIC code. The numbers shown in the figure are obtained by assuming that the spatial scale along the z axis is equal to the initial width of the channel.

11/16
The vast majority of the nonlinear pairs are produced by photons emitted after the two lasers collide. Indeed, t emit = 0 is the time when the two laser beams begin to collide. Most of the pairs are located at t emit > 0, which indicates that the parent photons were generated after the collision.

Strong-field modifications to the linear Breit-Wheeler cross section
The cross section we use to determine the number of electron-positron pairs created by the linear Breit-Wheeler process, Eq. (1), is calculated assuming that the photons are in vacuum. However, our calculations show that the pairs are created within the plasma channel as the laser pulses overlap, where the background electromagnetic field is strong. The presence of a strong magnetic field 49,50 or plane EM wave is known to modify the cross section for this process. Calculations of the cross section in the latter case have focused on changes at moderate a 0 62 or on resonance features 51,52 . Resonances occur where the intermediate fermion is on mass shell, which corresponds to the incoherent combination of nonlinear pair creation, followed by photon absorption by one of the daughter fermions, in the high-intensity regime. This contribution, which is effectively driven by a single photon, is already counted as our simulations include nonlinear pair creation; photon absorption by an electron in a strong field 9, 63 is suppressed unless the photon and electron are aligned within electron's emission cone. (See supporting simulations in Ref. [64].) As such, in order to estimate the effect of the strong field at a 0 1, we use the cross section for two-photon pair creation in a constant, crossed field given in Sec. 5.7 of Ref. [53]. If two-photon pair creation is kinematically allowed, i.e. ς > 1, corrections to the cross section scale as (χ γ /ς ) 2 , where χ γ is the quantum nonlinearity parameter 53 : in the scenario under consideration here, the photons which undergo linear pair creation have MeV energies and χ γ 1, and therefore these corrections can be neglected. On the other hand, the fact the laser-plasma interactions provide a platform for prolific two-photon pair creation in a region of strong EM field, as our results show, motivates a more general treatment that can investigate where these field-driven corrections become substantial. We note that a related process, pair annihiliation to two photons in a pulsed, plane EM wave, has recently been revisited 56 .

Pair production algorithm for head-on collisions of beamlets
Our algorithm for computing the pair production yield due to the linear Breit-Wheeler process leverages the fact that the colliding photons are emitted over an extended period of time compared to the characteristic travel time between the emission locations. In what follows, we illustrate its implementation for head-on collisions of two beamlets. Near head-on collisions are the biggest contributor to the pair yield and this is also the regime that greatly benefits from our simplified approach in terms of computational efficiency.
We are considering a head-on collision of two counterpropagating beamlets that have the same transverse area S. The first beamlet is being emitted at x = x 1 and it propagates in the positive direction. It is convenient to use the emission time τ as a marker for the photons in each beamlet. The corresponding photon density in the first beamlet is then n 1 (τ 1 ). The counterpropagating beamlet is emitted at x = x 2 > x 1 and its photon density is n 2 (τ 2 ). The total pair yield by these two beamlets interacting with each other is where l ≡ x 2 − x 1 .
Equation (4) presents a direct approach to calculating the pair yield and it involves a double integral. In our case, the typical beamlet emission lasts longer than l/c, which suggests a possible simplification of replacing n 2 (τ 2 ) with n 2 (τ 1 ). The inner integral in Eq. (4) can then be directly evaluated and we find that where V ≡ lS is the interaction volume. Note that the two beamlets are interchangeable in this expression. In order to estimate the error, we expand n 2 (τ 2 ) around τ 2 = τ 1 and retain linear and quadratic terms in the expansion. The time integral of the linear term in Eq. (4) is equal to zero, which means that the error in our approach is determined by the quadratic term. We thus estimate that the relative error in the number of produced pairs scales as (l/c∆τ) 2 , where ∆τ is the characteristic duration of beamlet emission. We compute the spatial distribution of pairs by simply assigning the pair density ∆n pairs = ∆N pairs /V

12/16
to each position within the interaction volume. We call this the average density method. A direct approach would however require us to compute the following integral at each position along the x axis: where The direct approach is much more demanding computationally. Even though our approach for calculating the spatial distribution of pairs is deliberately crude for a single pair of beamlets, it is effective when applied to a large ensemble of spatially distributed beamlets. As an example, we have carried out calculations for approximately 3000 spatially distributed beamlets, where 1573 beamlets are directed to the right and 1590 beamlets are directed to the left. The photon energy range is 0.19 MeV < ε γ < 2.05 MeV. We used our PIC simulation to generate these beamlets by selecting photons emitted with |θ | < 5.15 • or |θ − π| < 5.15 • . The average density method that uses ∆n pairs from Eq. (7) for each beamlet-beamlet collision gives the curve shown with small solid circles. The direct approach detailed by Eq. (8) gives the curve shown with open circles. The two curves have a similar shape, but the direct approach took almost two orders of magnitude longer in terms of computational time. The relative difference in the total number of pairs between the two methods is less than 17%.

! [µm]
Our method evenly distributes the generated pairs over the interaction volume, which is the key to achieving a good agreement with the direct but more computationally expensive approach. In order to illustrate this aspect, we performed another calculation. In this case, all of the pairs produced by two beamlets are placed into the center of the interaction volume without being evenly distributed. We call this the point density method. The density is calculated as ∆N pairs /S∆x, where ∆N pairs is given by Eq. (6) and ∆x is the thickness of slices that we use for spatial discretization into beamlets. The result of this procedure is shown with star markers in Fig. 13. There are no significant savings in terms of computational costs compared to the average density method. The characteristic width of the spatial distribution shows considerable deviation from that for the direct approach.