Optimal laser focusing for positron production in laser–electron scattering

In an intense electromagnetic background, charged particles obtain relativistic velocities and emit energetic photons. A fraction of these photons decays into electron–positron pairs, which can themselves be accelerated by the fields and radiate new photons [1, 2]. The repeated recurrence of photon emission and pair creation leads to the formation of the so-called QED cascades, where the number of particles in a plasma grows exponentially with time.

Pair cascades are believed to affect plasma dynamics in extreme astrophysical environments, e.g., in pulsar magnetospheres and polar caps [3–5]. It was recently proposed that one could re-create a comparable energy density in the laboratory using counter-propagating intense laser pulses [6], which has prompted many scientists to study related configurations using kinetic particle-in-cell simulations [7–27]. As the required laser intensity (I ∼ 10²⁴ W cm⁻²) is still beyond the extent of the current laser technology, there are many unknowns about the highly nonlinear dynamics associated with plasmas in extreme conditions. Previous proposals [28, 29] have shown that the required laser intensity for copious GeV pair production in near-critical-density plasmas can be reduced to the order of I ∼ 10²² W cm⁻².

Before lasers become sufficiently intense to generate dense e⁺ e⁻ pair plasmas from light, a head-on collision between a pulsed laser and a very energetic electron beam can allow us to generate dilute e⁺ e⁻ beams by applying the currently available technology [30]. The famous SLAC E-144 experiment has shown the onset of nonlinear Compton scattering and Breit–Wheeler pair production combining a 10¹⁸ W cm⁻² laser with a $\sim 50\enspace \mathrm{G}\mathrm{e}\mathrm{V}$ electron beam [31, 32], albeit with a very low positron yield. Experiments are underway, expected to produce more pairs per shot [33]. Research on stochastic effects in radiation reaction is also expected to benefit from the laser–electron scattering experiments [34–38], with new ways to infer the peak laser intensity at the interaction point [39, 40] and probing the transition from the classical to the quantum-dominated laser–electron interaction. Two all-optical experiments have shown the electron slowdown due to radiation reaction [41, 42], but were not able to discriminate between different theoretical descriptions of radiation reaction. We anticipate the near-future facilities (e.g. ELI [43], Apollon [44], CoReLS [45], FACET-II [46, 47], LUXE [48], EXCELS [49], ZEUS [50]) are to probe the electron–positron pair production covering several different regimes of interaction. This manuscript focuses on head-on laser–electron scattering that maximizes the strength of the electric field in the electron rest frame. This is the first experiment planned in most of the aforementioned facilities, and we aim to improve the current predictive capabilities for positron creation.

Due to the inherent non-linearity of the Breit–Wheeler pair production, there is no general roadmap on what would be an optimal strategy to obtain the highest possible positron yield using any given laser system. If the laser is assumed to be a plane wave (adequate when the laser is much wider than the interacting beam), the analytical predictions state that the best strategy would be to use the highest conceivable laser intensity. Therefore it is tempting to conclude that the laser should be focused on the smallest attainable focal spot. Our work shows that this strategy may not always be optimal, as there is a trade-off between the high laser intensity and the size of the interaction volume. With a short focus, the highest intensity region becomes small both transversely and longitudinally, which can reduce the number of seed electrons that interact with the close-to-the-maximum intensity, as well as the duration of this interaction. Using tight focusing also increases the number of particles that are not temporally synchronized with the peak of the laser pulse at the focal plane (these electrons never get to interact with the peak laser intensity, regardless of the transverse position). Finally, the wavefront curvature can also change the effective angle of laser–particle interaction.

Each of the mentioned factors affects the resulting number of positrons; hence a correct optimization strategy would have to take all of them into account at the same time. This can be achieved by resorting to full-scale three-dimensional particle-in-cell simulations, making sure enough statistics is used to represent the interacting electron beam with all its features, as well as high grid resolution in all spatial directions to correctly describe the laser dynamics. This approach requires a lot of computational resources (several million CPU-hours for each parameter set) and can be justified for the support of a specific ongoing experiment where most parameters are not free. However, for future experiments, there are many possible choices. It would thus be practical to devise a simple and cost-effective way for their consideration. This would ensure that the best possible strategies are applied when constructing new laser facilities.

This manuscript aims to simplify the pair production optimization in multi-variable parameter space. The goal is accomplished by extending the analytical scaling laws previously developed for electron interaction with a plane wave to more realistic geometries, taking into account the finite electron beam size and the laser focusing. Our method allows for predicting the number of positrons created in a laser–electron collision with a temporal (longitudinal) or a perpendicular offset. With the model proposed in this paper, predictions can be obtained without simulations. The most demanding calculation required is a numerical integration of an analytical function that can be performed on a single CPU, and the results obtained within minutes. This paper is structured as follows. In section 2, we revisit the pair production in a plane wave. Section 3 covers different methods for calculating the overall positron yield. In section 4, we give predictions for a special case where the electron beam is long and wide compared to the laser beam at the focus. Section 5 discusses a thin electron beam, while section 6 features a short electron beam. Finally, in section 7, we show the optimization of the positron number expected as a function of the available laser and electron beam parameters.

The simplest description of a laser pulse is a plane wave with a temporal envelope. Such a wave is fully described by the wavelength λ, a pulse duration τ (that defines the extension of the pulse's temporal envelope), and the normalized vector potential a₀, which relates to the intensity through ${a}_{0}=0.855\sqrt{I[1{0}^{18}\enspace \mathrm{W}\enspace {\mathrm{c}\mathrm{m}}^{-2}]}\lambda \enspace (\mu \mathrm{m})$ (for linearly polarized lasers). As a relativistic electron interacts with the strong electromagnetic wavepacket, it emits high-energy photons that themselves interact with the laser field and can decay into electron–positron pairs. The process of electron–positron generation is mediated by a real photon through the Breit–Wheeler mechanism [51].

In the plane wave approximation, the total number of new pairs per interacting electron can be estimated if we know the electron energy γ₀ m_e c² (where γ₀ is the initial electron Lorentz factor, m_e the electron mass and c is the speed of light), and the laser parameters (peak a₀, central wavelength λ and pulse duration τ, which is defined as the full width at half maximum of the laser intensity).

The total number of pairs is then given by [52]:

$$\begin{equation}{\left.{N}_{+}^{\mathrm{P}\mathrm{W}}({\gamma }_{0},{a}_{0},\lambda ,\tau )\simeq 3\sqrt{\frac{\pi }{2}}{P}_{\pm }({\omega }_{\mathrm{c}})\enspace {\chi }_{\mathrm{c},\mathrm{r}\mathrm{r}}\frac{{({\gamma }_{0}m{c}^{2}-\hslash {\omega }_{\mathrm{c}})}^{2}}{\hslash {\gamma }_{0}m{c}^{2}}\frac{\mathrm{d}{N}_{\gamma }}{\mathrm{d}\omega }\right\vert }_{\omega ={\omega }_{\mathrm{c}}}\end{equation} \tag{ 1 }$$

The first term P_±(ω_c) represents the probability of emitting a photon of frequency ω_c; the second is the recoiled χ_c,rr which accounts for the radiation reaction on the beam electrons and the final term dN_γ /dω is the value of the emitted photon distribution at ω = ω_c. This approximation underestimates the number of low-energy photons, which does not significantly affect the positron production calculation. According to this model, all positrons are generated from photons with a critical frequency ω_c, and there is no feedback by the produced pairs on the photon spectra (in other words, there is only one generation of secondary particles). Furthermore, the model assumes a semi-classical equation of motion of the electrons as they lose energy through the emission of radiation and uses the locally constant field and rigid-beam approximations. This allows for an implicit calculation of the laser phase ϕ_c at which χ is maximized, which is later used to estimate the emitted photon spectrum and consequently the number of created pairs (figure 1).

Zoom In Zoom Out Reset image size

Let us now consider a diffraction-limited laser pulse illustrated in figure 3. The maximum laser intensity an individual particle within the electron beam interacts with depends on two geometrical factors: the transverse offset from the laser axis compared to the laser spotsize and the initial longitudinal position that affects the temporal synchronization of the interaction. In other words, while interacting with a Gaussian laser pulse, electrons far from the focus interact with a lower average (and peak) field, which must be taken into account. The electron encounters the peak of the laser pulse at time t in an (x, y, z) point of configuration space which defines the maximum field felt by this particle. We can therefore assign an effective vector potential a_0,eff(t, x, y, z) that corresponds to the maximum laser intensity the particle experiences during the interaction.

We define an equivalent distribution of beam particles according to the maximum intensity they interact with during the scattering. This intensity is identified through the maximum instantaneous vector potential associated with an individual beam particle as a_0,eff. In the case of a plane wave interaction, there is no defocusing, and particles always interact with the same intensity, regardless of where or when they overlap with the peak of the laser (a_0,eff ≡ a₀, and the equivalent distribution, in this case, would be a dirac Delta function). For a more general case, by considering a corrected a_0,eff for each particle, we can apply the equations already derived for a plane wave (equation (1)), and then integrate over the distribution function in a_0,eff to obtain the total yield of positrons in the laser–electron scattering. The integration can be performed by sampling the distribution function numerically or performing an analytical integration over the configuration space in some special cases where this is possible. In reference [52] the authors calculate the total number of positrons by analytically and numerically integrating the scaling law for the plane wave in full coordinate space x, y, z (for comparison with our predictions cf section 6). The approach in the present manuscript allows to make this procedure shorter and applicable to different beam shapes and sizes. After deriving the distribution of particles dN/da_0,eff a one-dimensional integration in a_0,eff-space is enough to make predictions on the total number of expected positrons (or some other relevant quantity).

We first obtain a distribution of the interacting particles according to their a_0,eff. For every bin in this distribution dN(a_0,eff)/da_0,eff one can calculate the contribution for pair production using ${N}_{+}^{\mathrm{P}\mathrm{W}}({\gamma }_{0},{a}_{0,\mathrm{e}\mathrm{ff}},\lambda ,\tau )$. This method is cost-effective because it eliminates the need to perform multiple variable integration in the configuration space.

The problem can be addressed using cylindrical coordinate system (ρ, ϕ, z), centered at the laser focus. For a Gaussian laser beam (in the paraxial approximation), the configuration space can be mapped according to the laser intensity isosurfaces shown in figure 2, that do not depend on the coordinate ϕ. For simplicity, let us first assume that the electron beam is a cylinder with a constant density n_b. Each particle meets the laser beam at a different point of space, and is assigned a_0,eff(ρ, z), where ρ and z are its coordinates at the instant of time when it is synchronised with the peak of the laser. Performing a one-to-one mapping to the new coordinates of a flat-top relativistic beam in counter-propagation with the laser, the beam density in the new coordinates doubles and the length shrinks by two because the laser–electron crossing happens at twice the speed of light. The number of particles dN_b(a_0,eff) with a_0,eff that falls in the interval (a_0,eff, a_0,eff + da_0,eff) can then be estimated to be dN_b(a_0,eff)/da_0,eff = 2n_bdV/da_0,eff, where dV is the volume between two adjacent isosurfaces associated with a_0,eff and a_0,eff + da_0,eff. Due to the geometry of the problem, this expression can be transformed to the following:

$$\begin{equation}\frac{\mathrm{d}{N}_{\mathrm{b}}({a}_{0,\mathrm{e}\mathrm{ff}})}{\mathrm{d}{a}_{0,\text{eff}}}={\int }_{S}\frac{2{n}_{\mathrm{b}}\enspace dS}{{\Vert}\nabla {a}_{0,\text{eff}}{\Vert}}\end{equation} \tag{ 2 }$$

where the surface element $dS=\rho \sqrt{\mathrm{d}{\rho }^{2}+\mathrm{d}{z}^{2}}\mathrm{d}\phi =\rho \sqrt{1+{(\partial \rho /\partial z)}^{2}}\mathrm{d}z\enspace \mathrm{d}\phi$ is calculated at the isosurface that is by definition perpendicular to the gradient of the vector potential given by ${\Vert}\nabla {a}_{0,\text{eff}}{\Vert}=\sqrt{{\left(\partial {a}_{0,\text{eff}}/\partial \rho \right)}^{2}+{\left(\partial {a}_{0,\text{eff}}/\partial z\right)}^{2}}$.

Zoom In Zoom Out Reset image size

Letting beam plasma density vary in space ${n}_{\mathrm{b}}(\overrightarrow{r})$ allows considering cases of short or long, wide or narrow beams, including non-ideal spatio-temporal synchronization with the laser (as discussed later for cases illustrated in figure 3). It is worth noting that even a point-particle interaction with a Gaussian beam is not equivalent to a plane wave approximation unless the particle is in perfect temporal synchronization with the laser pulse.

Zoom In Zoom Out Reset image size

An alternative way to obtain the distribution of particles in a_0,eff is to numerically sample the electron beam density n_b(r, θ, ϕ) in space (the 'sampling' method). The sampled distribution function can be directly binned into a histogram according to the maximum value of a₀ each section interacts with. This allows to take into account arbitrary beam shapes, as well as include spatial structure like correlated energy chirp that can be present in LWFA beams. As long as the equation (1) for the number of expected pairs in a plane wave interaction is correct, the 'sampling' method is expected to predict a correct result for a focused Gaussian laser beam interacting with an electron beam of any shape and size. In the rest of the manuscript, we use it as a verification of the explicit distribution functions obtained analytically for several typical cases expected in an experiment. A good agreement between the predictions of the two methods confirms that the analytical distribution functions are correct and that the approximations taken in the electron beam description are valid.

Once the particle distribution in equation (2) is calculated, we can extract field moments $\langle {{a}_{0,\text{eff}}}^{k}\rangle =\int {{a}_{0,\text{eff}}}^{k}\enspace \mathrm{d}N/\mathrm{d}{a}_{0,\text{eff}}\enspace \mathrm{d}{a}_{0,\text{eff}}$, which can for example be used to calculate the average laser intensity [40].

As a first application of the ideas presented in the last section, let us consider the case of the scattering between a focused Gaussian laser pulse and a wide electron beam. The spatio-temporal intensity distribution of a Gaussian laser is characterized by the peak vector potential a₀, the laser wavelength λ, and the Rayleigh-range ${z}_{\mathrm{R}}=\pi {W}_{0}^{2}/\lambda$, where W₀ is the transverse spot size. The effective vector potential has the following spatial dependence ${a}_{0,\mathrm{e}\mathrm{ff}}=\left({a}_{0}/\sqrt{1+{(z/{z}_{\mathrm{R}})}^{2}}\right)\mathrm{exp}\left(-({\rho }^{2}/{W}_{0}^{2})/(1+{(z/{z}_{\mathrm{R}})}^{2})\right)$, where z is the distance from the focal plane and ρ is the distance from the laser propagation axis. Our definition of 'a wide beam' is that the beam radius is much larger than the laser focal spot W₀. The gradient of a_0,eff can be written as ${\Vert}\nabla {a}_{0,\text{eff}}{\Vert}=\vert \partial {a}_{0,\text{eff}}/\partial \rho \vert \enspace \sqrt{1+{(\partial \rho /\partial z)}^{2}}$, where $\vert \partial {a}_{0,\text{eff}}/\partial \rho \vert =2\rho {a}_{0,\text{eff}}/({W}_{0}^{2}(1+{(z/{z}_{\mathrm{R}})}^{2}))$. This simplifies the particle distribution in a_0,eff according to equation (2):

$$\begin{equation}\frac{\mathrm{d}{N}_{\mathrm{b}}({a}_{0,\mathrm{e}\mathrm{ff}})}{\mathrm{d}{a}_{0,\text{eff}}}=\frac{2\pi {n}_{\mathrm{b}}{W}_{0}^{2}}{{a}_{0,\mathrm{e}\mathrm{ff}}}{\int }_{{z}_{\mathrm{min}}}^{{z}_{\mathrm{max}}}1+{\left(\frac{z}{{z}_{\mathrm{R}}}\right)}^{2}\enspace \mathrm{d}z\end{equation} \tag{ 3 }$$

where the limits of integration in z direction will depend on the beam length and its temporal synchronization with the laser pulse. If the entire isosurface associated with a specific a_0,eff is covered with interacting particles, ${z}_{\mathrm{max}}=-{z}_{\mathrm{min}}={z}_{\mathrm{R}}\sqrt{{({a}_{0}/{a}_{0,\text{eff}})}^{2}-1}$. Otherwise, a portion of the volume associated with a specific laser intensity may be empty due to the finite beam length and temporal synchronization. The interaction limits imposed by the beam are ${z}_{\mathrm{min}}^{\ast }={{\Delta}}_{{\Vert}}-L/4$ and ${z}_{\mathrm{max}}^{\ast }={{\Delta}}_{{\Vert}}+L/4$, where Δ_∥ is the longitudinal displacement of the electron beam centre from the focal plane when the laser is at the focus. For every a_0,eff, one has to evaluate what the appropriate integration limits are on each side, and there is a transition in the distribution function at values of a_0,eff that corresponds to the beam edge on axis. In the case where Δ_∥ = 0, the a_0,eff above which the iso-volumes are completely full is ${a}_{z}\equiv {a}_{0}/\sqrt{1+{(L/4{z}_{\mathrm{R}})}^{2}}$, as shown by the examples in figure 4(a). The distribution function for a centered, wide, flat-top electron beam is given explicitly in table 1. For clarity, the table summarizes the results only for temporally synchronised beams, but the model allows considering more general cases, that are covered in the appendix A.

We now illustrate the obtained particle distributions according to the effective laser intensity on an example. The SFQED experiment [46] will study pair-production using a 0.61 J laser pulse (a₀ = 7.3, λ = 0.8 μm, W₀ = 3 μm, τ = 35 fs) and a 13 GeV, 2 nC electron beam. The electron beam follows a non-symmetric Gaussian density distribution transversely with σ_x = 24.4 μm, σ_y = 29.6 μm, and has a $\sim 250\enspace \mu \mathrm{m}$ long flat-top longitudinal profile.

Zoom In Zoom Out Reset image size

Table 1. Particle distributions according to the effective vector potential for different beam geometries. Here, ${a}_{z}\equiv {a}_{0}/\sqrt{1+{(L/4{z}_{\mathrm{R}})}^{2}}$ is the a_0,eff associated with the integration limits imposed by the longitudinal size of the electron beam, N_b represents the total number of particles in the beam, n_b is the beam density, R and L are the beam radius and length respectively. The laser spot size is W₀, ${z}_{\mathrm{R}}\equiv \pi {W}_{0}^{2}/\lambda$ is the Rayleigh length, and Δ_⊥ is the perpendicular displacement of the beam centre from the laser propagation axis.

Setup	Particle distribution for temporally centered beams
Wide beam	$\frac{\mathrm{d}{N}_{\mathrm{b}}}{\mathrm{d}{a}_{0,\text{eff}}}=\begin{cases}\frac{4\pi \enspace {n}_{\mathrm{b}}\enspace {W}_{0}^{2}\enspace {z}_{\mathrm{R}}}{{a}_{0,\text{eff}}}\frac{\sqrt{{a}_{0}^{2}-{a}_{0,\text{eff}}^{2}}}{3{a}_{0,\text{eff}}}\left(2+{\left(\frac{{a}_{0}}{{a}_{0,\text{eff}}}\right)}^{2}\right),{a}_{0,\text{eff}}\geqslant {a}_{z}\quad \\ \frac{4\pi \enspace {n}_{\mathrm{b}}\enspace {W}_{0}^{2}\enspace {z}_{\mathrm{R}}}{{a}_{0,\text{eff}}}\frac{L}{4{z}_{\mathrm{R}}}\left(1+{\left(\frac{L}{4{z}_{\mathrm{R}}}\right)}^{2}\right),\enspace {a}_{0,\text{eff}}< {a}_{z}\quad \end{cases}$
Thin beam	$\frac{\mathrm{d}{N}_{\mathrm{b}}}{\mathrm{d}{a}_{0,\text{eff}}}=\begin{cases}\frac{4{N}_{\mathrm{b}}{z}_{\mathrm{R}}}{L}\frac{{a}_{0}^{2}}{{a}_{0,\text{eff}}^{2}}\frac{1}{\sqrt{{a}_{0}^{2}-{a}_{0,\text{eff}}^{2}}},{a}_{0,\text{eff}}\geqslant {a}_{z}\quad \\ \quad 0,\enspace {a}_{0,\text{eff}}< {a}_{z}\quad \end{cases}$
Short beam	$\frac{\mathrm{d}{N}_{\mathrm{b}}}{\mathrm{d}{a}_{0,\text{eff}}}=\frac{{N}_{\mathrm{b}}\enspace {W}_{0}^{2}}{{R}^{2}\enspace {a}_{0,\mathrm{e}\mathrm{ff}}}{\left(\frac{{a}_{0,\text{eff}}}{{a}_{0}}\right)}^{{({W}_{0}/R)}^{2}}{\mathrm{e}}^{-{{\Delta}}_{\perp }^{2}/{R}^{2}}\enspace {I}_{0}\left(\frac{2{{\Delta}}_{\perp }{W}_{0}}{{R}^{2}}\sqrt{\mathrm{log}\left(\frac{{a}_{0}}{{a}_{0,\text{eff}}}\right)}\right)$

To save computing time, we performed 3D PIC simulations of this interaction using OSIRIS [53] by dividing the long beam into five equal beamlets each 50.9 μm long. These beamlets have different temporal synchronization (they encounter the laser peak at different distances from the focus). The 3D simulations are performed with a box size of 98 μm × 25 μm × 25 μm, resolved with 3840 × 400 × 400 cells. OSIRIS PIC results (red empty circles in figure 4(b)) are compared with the analytical predictions based on above intensity distribution functions and a numerically sampled beam. The distribution functions for the temporally non-synchronised electron beam with Δ_∥ ≠ 0 are shown in the appendix A.

For the analytical calculations and numerical sampling, we assumed the beam has a uniform density equal to the central density of the electron beam n_b = 10¹⁶ cm⁻³. This is justified by σ_x ≫ W₀ and σ_y ≫ W₀, and the highest intensity portion interacts nearly exclusively with the maximum density of the beam. The analytical calculation is therefore coherent with the simulation results, as confirmed by the comparison in figure 4(b).

Let us now consider the case of the scattering between a focused Gaussian laser and a long L ≫ z_R but thin R ≪ W₀ electron beam, where the effective laser intensity can be lower due to a longitudinal offset Δ_∥. Here, the problem becomes one-dimensional as the number of particles in the volume associated with a specific value of intensity is given by dN_b = 2n_bdV = 4n_b S_⊥dz = 4(N_b/L)dz, provided that the beam is long enough to interact with the two portions of the same laser intensity before and after the focus (which is always the case if the beam is temporally centered—otherwise if the interaction happens only on one side, the total number of particles should be divided by two). The effective laser intensity depends only on z through ${a}_{0,\text{eff}}={a}_{0}/\sqrt{1+{(z/{z}_{\mathrm{R}})}^{2}}$ and the distribution can simply be calculated as

$$\begin{equation}\frac{\mathrm{d}{N}_{\mathrm{b}}}{\mathrm{d}{a}_{0,\text{eff}}}=\frac{4{N}_{\mathrm{b}}}{L}\frac{\mathrm{d}z}{\mathrm{d}{a}_{0,\text{eff}}}.\end{equation} \tag{ 4 }$$

Just as in the previous section, the explicit distribution function in a_0,eff for a centered beam is given in table 1 (figure 5).

Zoom In Zoom Out Reset image size

Let us now consider a short L ≪ z_R beam with a Gaussian density profile ${n}_{\mathrm{b}}={n}_{0}\enspace \mathrm{exp}\left(\right.-({(x-{{\Delta}}_{\perp })}^{2}+{y}^{2})/{R}^{2}\left.\right)$ where the peak density is given by n₀ = N_b/(πR² L), x = ρ cos ϕ and y = ρ sin ϕ. We assume a longitudinally synchronized beam (Δ_∥ = 0), with an allowed transverse displacement Δ_⊥ between the beam centre and the laser propagation axis. The electrons therefore interact with the laser peak at z = 0 and the field structure reduces to ${a}_{0,\text{eff}}={a}_{0}\enspace \mathrm{exp}\left(-{\rho }^{2}/{W}_{0}^{2}\right)$. As the manifolds of constant a_0,eff are now concentric rings, the volume element associated with a specific value of a_0,eff is given by dV = Lρdρ dϕ/2, the surface element of an isosurface is dS = Lρdϕ/2, and the field gradient is given by $\nabla {a}_{0,\text{eff}}=\partial {a}_{0,\text{eff}}/\partial \rho \enspace \hat{\rho }$, with $\partial {a}_{0,\text{eff}}/\partial \rho =-2\rho \enspace {a}_{0,\text{eff}}/{W}_{0}^{2}$. We can now apply the equation (2) to obtain the particle distribution function dN_b(a_0,eff)/da_0,eff = ∫Ln_b ρdϕ/||∇a_0,eff||. This gives

$$\begin{equation}\frac{\mathrm{d}{N}_{\mathrm{b}}}{\mathrm{d}{a}_{0,\text{eff}}}=\frac{L{n}_{0}{W}_{0}^{2}}{2{a}_{0,\mathrm{e}\mathrm{ff}}}\enspace {\text{e}}^{-{\rho }^{2}/{R}^{2}}\enspace {\mathrm{e}}^{-{{\Delta}}_{\perp }^{2}/{R}^{2}}{\int }_{0}^{2\pi }{\mathrm{e}}^{2\rho {{\Delta}}_{\perp }\enspace \mathrm{cos}\enspace \phi /{R}^{2}}\mathrm{d}\phi \end{equation} \tag{ 5 }$$

where $-{\rho }^{2}={W}_{0}^{2}\enspace \mathrm{log}\enspace ({a}_{0,\text{eff}}/{a}_{0})$ and the integration result can be expressed through the modified Bessel function of the first kind ${I}_{0}(t)=(1/\pi ){\int }_{0}^{\pi }\enspace \mathrm{exp}(t\enspace \mathrm{cos}\enspace \phi )\mathrm{d}\phi$. The obtained particle distribution function is given in table 1 and can be numerically integrated applying ${N}_{+}^{\mathrm{P}\mathrm{W}}$ to every bin of the histogram to obtain the total positron yield from the laser–electron beam interaction. The calculations can be extended to a temporally non-synchronized interaction by replacing a₀ and W₀ with ${a}_{0}^{\ast }={a}_{0}(z)={a}_{0}/\sqrt{1+{(z/{z}_{\mathrm{R}})}^{2}}$ and ${W}_{0}^{\ast }={W}_{0}(z)={W}_{0}\sqrt{1+{(z/{z}_{\mathrm{R}})}^{2}}$.

In reference [52], the authors consider a spherically-symmetric Gaussian beam profile with a radius R = 6 μm, the laser spotsize W₀ = 2 μm and the Rayleigh range z_R = 15.7 μm. As the Rayleigh range is much higher than the beam length (z_R ≫ R), one can consider the beam short, just like in our calculations. They have obtained an approximate expression for the expected value of the number of positrons, which correctly predicts the order of magnitude, but has a factor of two difference compared with the simulation results in reference [52]. We have applied our distribution functions to calculate the expected number of positrons in this case, and we have obtained the correct result that is consistent with the simulation data of reference [52]. Detailed comparisons are shown in figure 6.

Zoom In Zoom Out Reset image size

This section covers the optimal focusing strategy for a wide range of laser parameters (in particular as a function of total energy content and pulse duration), as well as different electron beam energies. We assume that the electron beam is 200 μm long (flat-top longitudinal envelope) and has a Gaussian transverse shape. The electron beam is spatio-temporally synchronized with the laser (i.e., the center of the beam interacts with the laser peak at the focal plane, and they share the propagation axis). The transverse beam profile is Gaussian with σ_x = 24.4 μm and σ_y = 29.6 μm. The chosen on-axis beam density n_b = 10¹⁶ cm⁻³ corresponds to the total beam charge of Q_b = 2 nC. The results can be scaled to other values for the central beam density by introducing a factor n_b/10¹⁶ cm⁻³. A specific laser system has a fixed total energy content, which for a Gaussian transverse profile is approximately given by $\varepsilon \enspace (\mathrm{J})\sim 2.1\times 1{0}^{-5}{a}_{0}^{2}{({W}_{0}/\lambda )}^{2}\enspace \tau \enspace (\mathrm{f}\mathrm{s})$. The laser intensity (proportional to ${a}_{0}^{2}$) is therefore reversely proportional to the square of the spot size W₀. As the number of pairs produced per interacting electron ${N}_{+}^{\mathrm{P}\mathrm{W}}$ is a monotonously rising function of the effective a₀, and the number of seed electrons that would experience the high intensity is proportional to the size of the interaction volume $\sim {W}_{0}^{2}{z}_{\mathrm{R}}$, to obtain the highest possible number of positrons, one should strike the right balance between a high value of a₀ and a large W₀. In other words, there is a trade-off between using a short focal length to obtain the highest conceivable laser intensity and having a wider interaction volume where more seed electrons participate in the interaction.

What follows is a calculation of the optimal focal spot and the corresponding pair yield for lasers with energy below 1 kJ and relativistic particle beams with energies lower or equal to 20 GeV. These values include what will soon be available in several experimental facilities (e.g. SLAC [46], HiBEF [54] or ELI [43]).

For each combination of the electron beam energy and the laser total energy content, we apply the analytical expression (see table 1) to calculate the effective a₀ distribution of the interacting particles. Then, we integrate the results numerically to find the optimal spotsize and maximum positron yield for this set of parameters (as illustrated in figure 7).

Zoom In Zoom Out Reset image size

Figure 8 summarizes the optimization results covering $\sim 1000$ different parameter combinations, keeping the laser duration constant at 35 fs. For 10 GeV electrons and a 1 kJ laser, a maximum number of pairs is 10⁹, which is obtained using W₀ > 8 μm. The FACET-II 13 GeV electron beam at SLAC could generate 4 × 10⁸ pairs/shot if paired with a 300 J laser-focused to W₀ = 5.7 μm. The LUXE 17.5 GeV beam with the same laser parameters could produce 7 × 10⁸ pairs per shot, using a slightly larger W₀ = 6.8 μm.

Zoom In Zoom Out Reset image size

Similarly, figure 9 shows how to obtain optimal results as a function of the laser energy and the laser pulse duration, keeping the initial electron beam energy constant at 13 GeV. This allows estimating the positron yield at ELI Beamlines, where L4 laser specifications are at 1.5 kJ with 150 fs duration. If we assume a third of the laser energy is used to accelerate electrons, 1 kJ is available for the scattering, which can produce 2.4 × 10⁸ (n_b/10¹⁶ cm⁻³) pairs per shot using W₀ = 6.2 μm.

Zoom In Zoom Out Reset image size

The results presented in this section best correspond to the case of a wide electron beam, which corresponds to a beam from a conventional accelerator. For LWFA beams, the 'thin' beam may be a more adequate description for the laser–electron interaction in some cases (e.g., if the scattering is performed within the LWFA bubble where the electron beam transverse size is $\sim 2\enspace \mu \mathrm{m}$). We can perform a similar optimization using the thin beam effective intensity distribution functions. If the LWFA beam propagates a few cm before the collision, due to the momentum divergence, the beam size increases, and even an LWFA beam is then well described by the wide beam distribution function.

As an illustrative example, in figure 10 we show the variation in the number of positrons for different values of a₀, for synchronized and an offset of 10 μm. In this case, the offset is much smaller than the Rayleigh range, which implies a small variation in the final positron yield.

Zoom In Zoom Out Reset image size

The highlighted cases were chosen to reflect the parameters we expect to be available in near-future laser facilities. As time delay and misalignment can be considered, in cases where there are uncertainties about specific experimental parameters, our models can be used as direct support for statistical analysis on pair production. We could incorporate fluctuations in intensity, temporal and spatial synchronization, as well as the fluctuations in electron beam energy distribution, and predict their impact on the final results.

The methods outlined in this work allow to make predictions regarding the pair production in laser–electron scattering, taking into account the 3D focusing geometry, spatio-temporal synchronization, and the realistic beam shape and size. This opens a possibility for fast parameter optimization, using analytical calculations directly or combining them with simple numerical integration. The approach is faster than using full-scale Monte Carlo or PIC-QED calculations, and the results can be obtained on a single CPU. Apart from saving computing resources, the ideas from the present study can be applied for real-time optimization and data analysis in experiments.

Besides optimizing the positron yield considered in this manuscript, the equivalent intensity distributions are going to be useful to calculate the asymptotic energy spread [34–37, 55] and divergence of the interacting electron beams [39, 56–58], which are also imprinted on the emitted photon beams in the hard x-ray and gamma-ray range. The analytical description can be generalized to the tight-focusing regime beyond the paraxial approximation (W₀ < 2.5λ) considering interaction at an angle for the regions with curved wavefronts. This requires introducing corrections on the effective τ and λ, and considering that the tilted interaction angle directly affects the quantum parameter χ. This will be explored in future publications.

In summary, the findings of this manuscript are expected to be important for near-future laser–electron scattering experiments, as the calculations are fast and can provide real-time feedback during the course of an experiment.

The authors thank B Martinez for proofreading the manuscript, as well as L O Silva and T Grismayer for fruitful discussions. This work was supported by the European Research Council (ERC-2015-AdG Grant No. 695088) and Portuguese Science Foundation (FCT) Grant No. CEECIND/01906/2018. We acknowledge PRACE for awarding access to MareNostrum based in the Barcelona Supercomputing Centre.

The data that support the findings of this study are available upon reasonable request from the authors.

In this section, we derive the distribution of particles for an arbitrary temporal synchronization. Similarly, as for the synchronized beams, we map the individual beam particles according to the spatial positions where they encounter the peak of the laser pulse. The temporal synchronization is then quantified through a longitudinal offset of the electron beam center from the focal plane given by Δ_∥. Without the loss of generality, we assume that Δ_∥ > 0.

To obtain a distribution function in a_0,eff for Δ_∥ > 0, one needs to evaluate the integral on the RHS of equation (2) between z₋ and z₊, where |z₋| < |z₊| are the longitudinal coordinates of the edges of the electron beam, given by z_± = Δ_∥ ± L/4. This integration can be performed formally, but we present a more intuitive calculation relying on geometry. For Δ_∥ < L/4, a part of the beam is temporally synchronized with the laser, and some of the particles have their interaction point at the focal plane. This is illustrated in figure 11(a). The beam can be divided into two parts: the blue part on the left-hand side of the focal plane and the yellow part on the right-hand side of the focal plane. The total distribution in a_0,eff can then be written as a sum of the distributions coming from the blue and the yellow sections of the beam. Each of these beams contributes with exactly a half of the distribution function associated with a symmetrical (temporally synchronized) beam twice its size. This means we can re-use the distribution functions from table 1, modifying the beam lengths to L_b/4 = |z₋| = Δ_∥ − L/4 and L_g/4 = |z₊| = Δ_∥ + L/4. For a large temporal offset Δ_∥ > L/4, none of the beam particles interact with the laser in the focus. This is illustrated in figure 11(b). Analogously as in (a), this case can be represented by the symmetric beams distribution functions defined through z_±, but this time we subtract the contribution of the blue beam from the contribution of the yellow beam.

Zoom In Zoom Out Reset image size

This reasoning applies both to thin and wide beam geometries. In the case of a short beam, the maximum interaction value for a_0,eff will be determined by the distance from the focus through ${a}_{{\Vert}}={a}_{0}/\sqrt{1+{({{\Delta}}_{{\Vert}}/{z}_{\mathrm{R}})}^{2}}$, and the effective spotsize becomes ${W}_{0}\to W({{\Delta}}_{{\Vert}})={W}_{0}\sqrt{1+{({{\Delta}}_{{\Vert}}/{z}_{\mathrm{R}})}^{2}}$. The modified distribution functions in a_0,eff including the temporal offset were collected in table 2 for easy reference.

Table 2. Particle distributions for arbitrary temporal synchronization. The ± sign corresponds to situations where Δ_∥ < L/4 or Δ_∥ > L/4 respectively, z_± = Δ_∥ ± L/4. θ(x) is the Heaviside theta function, ${a}_{{z}_{\pm }}={a}_{0}/\sqrt{1+{({z}_{\pm }/{z}_{\mathrm{R}})}^{2}}$ and ${a}_{{\Vert}}={a}_{0}/\sqrt{1+{({{\Delta}}_{{\Vert}}/{z}_{\mathrm{R}})}^{2}}$.

Setup	Particle distribution for temporally unsynchronized beams
Wide beam	$\frac{\mathrm{d}{N}_{\mathrm{b}}}{\mathrm{d}{a}_{0,\text{eff}}}=\frac{2\pi \enspace {n}_{\mathrm{b}}\enspace {W}_{0}^{2}\enspace {z}_{\mathrm{R}}}{{a}_{0,\text{eff}}}\left(\frac{{z}_{+}}{{z}_{\mathrm{R}}}\left(1+\frac{1}{3}{\left(\frac{{z}_{+}}{{z}_{\mathrm{R}}}\right)}^{2}\right)\theta ({a}_{{z}_{+}}-{a}_{0,\text{eff}})\right.$
	$\enspace \enspace \enspace \enspace \enspace \enspace \enspace \enspace \enspace \enspace \enspace \enspace -\frac{{z}_{-}}{{z}_{\mathrm{R}}}\left(1+\frac{1}{3}{\left(\frac{{z}_{-}}{{z}_{\mathrm{R}}}\right)}^{2}\right)\theta ({a}_{{z}_{-}}-{a}_{0,\text{eff}})$
	$\left.\enspace \enspace \enspace \enspace \enspace \enspace \enspace \enspace \enspace \enspace \enspace \enspace +\frac{\sqrt{{a}_{0}^{2}-{a}_{0,\text{eff}}^{2}}}{3{a}_{0,\text{eff}}}\left(2+\frac{{a}_{0}^{2}}{{a}_{0,\text{eff}}^{2}}\right)\left(\theta ({a}_{0,\text{eff}}-{a}_{{z}_{+}})\pm \theta ({a}_{0,\text{eff}}-{a}_{{z}_{-}})\enspace \right)\right)$
Thin beam	$\frac{\mathrm{d}{N}_{\mathrm{b}}}{\mathrm{d}{a}_{0,\text{eff}}}=\frac{2{N}_{\mathrm{b}}{z}_{\mathrm{R}}}{L}\frac{{a}_{0}^{2}}{{a}_{0,\text{eff}}^{2}}\frac{1}{\sqrt{{a}_{0}^{2}-{a}_{0,\text{eff}}^{2}}}\left(\theta ({a}_{0,\text{eff}}-{a}_{{z}_{+}})\pm \theta ({a}_{0,\text{eff}}-{a}_{{z}_{-}})\enspace \right)$
Short beam	$\frac{\mathrm{d}{N}_{\mathrm{b}}}{\mathrm{d}{a}_{0,\text{eff}}}=\frac{{N}_{\mathrm{b}}\enspace {W}_{0}^{2}}{{R}^{2}\enspace {a}_{0,\mathrm{e}\mathrm{ff}}}{\left(\frac{{a}_{0,\text{eff}}}{{a}_{0}}\sqrt{1+{\left(\frac{{{\Delta}}_{{\Vert}}}{{z}_{\mathrm{R}}}\right)}^{2}}\right)}^{{\left(\frac{{W}_{0}}{R}\right)}^{2}\left(1+{\left(\frac{{{\Delta}}_{{\Vert}}}{{z}_{\mathrm{R}}}\right)}^{2}\right)},\enspace \mathrm{i}\mathrm{f}\enspace {a}_{0,\text{eff}}< {a}_{{\Vert}}$