Ultra-bright GeV photon source via controlled electromagnetic cascades in laser-dipole waves

One aim of upcoming high-intensity laser facilities is to provide new high-flux gamma-ray sources. Electromagnetic cascades may serve for this, but are known to limit both field strengths and particle energies, restricting efficient production of photons to sub-GeV energies. Here we show how to create a directed GeV photon source, enabled by a controlled interplay between the cascade and anomalous radiative trapping. Using advanced 3D QED particle-in-cell (PIC) simulations and analytic estimates, we show that the concept is feasible for planned peak powers of 10 PW level. A higher peak power of 40 PW can provide $10^9$ photons with GeV energies in a well-collimated 3 fs beam, achieving peak brilliance ${9 \times 10^{24}}$ ph s$^{-1}$mrad$^{-2}$mm$^{-2}$/0.1${\%}$BW. Such a source would be a powerful tool for studying fundamental electromagnetic and nuclear processes.

of radiation. The dipole wave can be formed by orienting, timing and focusing a number of laser pulses [14] as shown in Fig. 1. In this configuration, ART traps particles and makes them oscillate in spatiotemporal regions ideal for gaining a maximal possible energy. Furthermore, within each oscillation of the field the particles have a high probability of emitting a large portion of their gained energy in the form of a single photon. However, the dominance of this favourable single-particle dynamics requires a total power of around 100 PW. As we will see further, even for a lower power a significant number of the generated high-energy photons can decay into e + e − pairs giving rise to a cascade of pair production. As a result, before the peak intensity of the field is reached the cascade can generate an e + e − plasma dense enough to not only restrict the growth of the field intensity [7], but also to terminate the favourable ART-regime dynamics. By assessing these processes with advanced simulations, we reveal here that for powers of the 10 PW level a new physical scenario arises in which the cascade and ART in fact induce and support each other.
We demonstrate that controlling the cascade development by matching the laser pulse intensity and duration in realistic ranges makes it possible to create a unique and ultra-efficient photon source in which laser radiation is converted into a well-collimated flash of GeV photons with unprecedented brightness.
The key to our concept is matching the pulse duration, peak power and initial particle density in the target, such that the maximal field intensity is reached just before collective plasma effects start to cause a significant reduction in the energy of generated photons.This requires that particle growth in the leading half of the laser pulse be restricted, such that the number of particles remains below a certain critical value N max which would cause significant back-reaction on the field. By ensuring that this constraint is fulfilled, we can arrange for a maximal number of particles to interact with the most intense part of the laser pulses, and emit a large number of high-energy photons, as desired. To perform a theoretical analysis and thus obtain estimates for the above constraint we begin by describing how particle dynamics in the ART regime [14] governs the cascade.
The motion of particles in the ART regime, illustrated in Fig. 1, is quasi-periodic and consists of two qualitatively different phases. During the half-cycle when the electric field does not change sign, a trapped particle is accelerated and gains a kinetic energy of around amc 2 before starting to emit frequently (hereafter a ≈ 800 P/ (1PW) is the electric field amplitude in relativistic units for a dipole wave of total incoming power P); we refer to this as the acceleration phase. As the magnetic field strength rises the particle loses a majority of its energy, having possibility to emit photons with energies of up to amc 2 , and then starts gyrating; we refer to this as the turning phase, and this lasts until the electric field of the next acceleration phase starts to force the particle back in the opposite direction. Due to the probabilistic nature of emission in quantum electrodynamics, it may happen that the particles do not experience sufficient losses to be forced back, and hence leave the trapped state as shown in Fig. 1.
Using estimates for the magnetic field strength in dipole waves [22], one can obtain that for powers of order P = 10 PW the quantum efficiency parameter [16,23] χ attains values of order unity. This parameter is defined as the ratio χ = γF eff /(eE S ) of the force acting perpendicular to the direction of motion of the electron in its rest frame, to that force which it would experience in a field of Sauter-Schwinger strength E S = m 2 c 3 /eh ∼ 10 16 V/cm (e and γ are the electron charge and gamma factor respectively). The fact that χ 1 means that a significant part of an electron's energy can emitted as a single photon, which can then convert into an electron-positron (ee + ) pair. By equating the time needed for this to happen with the duration of the turning phase we can (see Methods B) estimate that this process gives rise to a pair production cascade for powers P > P min ∼ 15 PW. As we will see below, this analytically determined threshold value is in a good agreement with numerical results. However, as we demonstrate in Methods A, cascade evolution in our system is a very complex process which is not amenable to simple analytical estimates.
Thus in the interests of providing accurate predictions we turn to simulations.
We use the 3D PIC code PICADOR [24] with quantum extensions [15] (see Methods C) to simulate cascade development in a dipole wave of constant incoming power P, for several values of P. We calculate the growth rate Γ of particle number averaged over a half-period, using the relative increase of the number of trapped particles during a single turning phase. This takes into account both particle production and particles leaving the trapping state.
On time scales longer than a half-period, the particle number grows exponentially before backreaction sets in. For this exponential stage the rate Γ depends only on P and, according to the results shown in Fig. 2 (a), may be approximately fitted to the curve Γ(P) ≈ 3.
where T is the laser period and P min 7.2 PW is the numerically determined threshold value required for the number of particles to grow. We integrate the pair production rate over the leading half of a laser pulse with Gaussian profile, FWHM duration τ, and demand that by the instant of reaching the peak power the number of trapped particles rises from N 0 = 10 (taken to guarantee the cascade seeding) to N max , introduced above. By equating the Coulomb field of the trapped particles to a/2, i.e. by demanding no more than 50% collective effects (see agreement with numerical results in Fig. 2 (a), solid blue curve) we estimate and obtain the following optimal condition on the pulse duration τ, at P > P min : This curve is shown in Fig. 2 (b) and fits well with the numerically estimated (see methods) number of produced photons with energyhω > 2 GeV. As can be seen from the Figure, a large number of energetic photons are produced and escape the fields, despite a rapid cascade process. This is because their mean free path is roughly four times larger than that of electrons/positrons with the same energy, and in addition is proportional to (hω) 1/3 , which favours the escape of photons with higher energies.
Having now identified the required pulse parameters, we provide a concrete example of the possible source properties using a comprehensive PIC simulation. (For full simulation details see Methods.) The results are shown in Fig. 3 (see also Supplementary Video). We form the dipole wave using 12 laser pulses with a total peak power of 40 PW and each having 15 fs duration (FWHM for intensity), which is matched to the initial particle density of the target (a uniform sphere with diameter 3 µm and density 10 20 cm −3 ). Instead of the considered target, in reality adjusting density of stray particles in the chamber can be used to trigger the cascade starting from an arbitrarily small number of particles. (This is because the cascade can be even prevented by choosing high enough vacuum [9,25].) To account for realistic experimental conditions, and so demonstrate the robustness of our concept, the pulses are slightly de-phased (randomly within 1/30 of the wave period). The simulation shows that the concept provides 10 9 photons with energieshω > 2 GeV in a well collimated beam (< 5 • ) and of just 3 fs duration, achieving peak flux of 1.9 × 10 23 ph s −1 , brightness 7.4 × 10 26 ph s −1 mrad −2 mm −2 and brilliance 9 × 10 24 ph s −1 mrad −2 mm −2 /0.1%BW, exceeding peak values of laser Compton scattering sources in GeV range by several orders of magnitude [26]. Such a source with extraordinary high flux of GeV photons with a broad spectrum will open qualitatively new possibilities for studying photo-nuclear processes[1, 17, 18], such as triggering multiple transitions between short-living states. One could also take advantage of the even higher photon flux in the focal region, by placing nuclear targets in or as close as possible to the focal spot.

Particle dynamics and the cascade
Dynamics and radiation losses in the ART regime are highly nonlinear, as even qualitative explanations assuming continuous (classical) radiation losses show [14]. The discrete and stochastic nature of quantum emission does not destroy ART, but instead adds another layer of complexity.
To analyse typical particle motion in the ART regime we therefore use here a 'regularised' model of emission which is a course-graining of the full quantum model (as based on the locally constant field approximation and implemented in PIC codes). In the model, particles move exactly along their mean free paths. During this motion we integrate the probability of emission until it reaches unity, and then a photon is emitted with an energy equal to the average of that which it would have according to the QED rates. Local values of χ and γ are used throughout.
In Fig. 4 we demonstrate how particles are caught (green) and circulate in the trapped state (red). The regularised emission model leads to particles sitting on stable attractors, a behaviour which in reality would break off due to stochastic effects. Hence in order to demonstrate how particles leave the system (blue) we use the now common probabilistic model of emission regularly employed in extended PIC schemes [15]. In panel (a) we show the paths of particles as a function of their γ-factor and the effective magnetic field B eff , which is defined as that field which would provide the instantaneous transverse acceleration of the particle at the particle's position (in space and time). These two parameters define the mean free motion time t free between emissions (grayscale), while the average photon energy tends to roughly 1/4 of the particle energy in the limit χ 1. The dark region in the left-upper of the figure denotes the state in which emission is frequent and thus energy gains are suppressed. Hence the optimal way for a particle to penetrate into the high-gamma region (in order to emit high-energy photons) is to pass through the low-B region. As one can see, energy gain is suppressed by emission on the green trajectory, due to the strong magnetic field further from the z-axis. For the red trajectory, however, the particle gains more energy when bypassing the dark region through point (1), and then emits high-energy photons around point (2). This is because the magnetic field is weaker near the z-axis, as for the red trajectory (see panel (b)). In this way, the ART provides close-to-optimal particle dynamics as compared, for example, with the case of 'normal radiative trapping' [14]. That is why the dipole wave, apart from providing the highest possible electric field amplitude also enables the optimal strategy for particles to gain high energy while the magnetic field remains small. At the same time, ART causes particles to drift toward the optimal spatial region. Particles can occasionally emit very high energy photons, due to stochastic effects, and then leave the system: fortunately the cascade compensates this as it provides a self-sustained source of particles.
Pair production cascades are known to be complex processes [13,[27][28][29]. In our case this complexity can be seen in the black curve in the insert of Fig. 4 (a); the characteristics of the particle growth rate vary over the entire wave period. This is due to the rate being a function of where particles are located in the B-γ, and this 'position' varies greatly (note though that the particles are almost always in the region χ > 1, as denoted by the dashed curve). In particular, the cascade is dominantly of avalanche type until the point y (as the particles gain energy), while in the interval between points y and z the cascade is dominantly of shower type [30]. These two types cannot be clearly distinguished or separated, though. Moreover, the contribution of a particle to the cascade depends strongly on its distance from the z-axis (i.e. on the type of particle motion), while the spatial distribution of particles is determined by the stochasticity of emission and the rate of particle production in different states. This makes the system particularly intractable to simple analytic estimates, and hence an accurate analysis requires numerical simulations. However, particle motion, particle generation, and events of particles leaving the system are all governed by a 'self-contained' process that is, as we will see below, defined only by the parameter P. That is why the presented analysis for all values of 0 < P < 100 PW has a general relevance.
The simulations indicate that the cascade and ART appear in an interdependent way already at the power of just 7.2 PW, which is significantly lower than the 100 PW needed for the dominance of ART itself [14]. , r e f f ≈ 0.15λ (see Fig. 4 (b). As one can see from Fig. 4 of reference [14], r e f f does not depend strongly on P). Using the analytical expressions for the dipole wave [22] we can estimate that in this spatiotemporal region the effective magnetic field in relativistic units is B e f f ≈ 2.3ar e f f /λ ≈ 0.3a. For our estimate, we assume that a magnetic field of this strength acts for 1/8 of the wave period, which we associate with the typical duration of the turning phase. We also assume that particles and photons have an energy order of mc 2 a.
We can then estimate χ ≈ 0.3a 2 (mc 2 /hω). For the considered wavelength λ = 910 nm we have χ ≈ 0.5P. Thus for values of P 10 PW we can assume χ 1. In this case, we can neglect the time of photon generation, as it is roughly four times less than the one of pair production [23], and estimate the threshold P for the cascade to arise by requiring the average time of pair production to be equal to one eighth of the wave period. Using the large-χ approximation for pair production rate [23], we obtain the relation 0.37 αm 2 c 4 hmc 2 a χ 2/3 1 8 2π Using the relation a ≈ 800 P/(1PW) yields

Simulations
Numerical simulation of the above described process requires accurate accounting for the radiation losses in the entire range of regimes from weak classical (χ 1) to strong quantum (χ 1).
The former is important before the particles come close to the z-axis (point (3) of the green trajectory in fig. 4) or when the high field values have not been reached. The latter is important for simulating trapped and leaving particles. That is why for this study we developed a modified event generator [15] that reproduces exactly the particle emission without any low-energy cutoff (common for all previously known methods) for all values of χ. Thus, in particular, it reproduces the classical emission in the limit χ 1. Despite of the fact that our method provides close to the minimal possible requirements to the time step, the interval between emissions and pair production can be extremely small in the central region of interaction. To bring the computational demands into realistic range we account for an arbitrary number of QED events within a single time-step using an adaptive event generator [15]. Additionally we apply procedures for resampling the ensemble of particles to avoid the memory overload through a gradual increase of the particles/photons computational weights. The validity of our model for the considered process is accurately demonstrated in the specially dedicated paper [15], where we thoughtfully analyzed all the methodological and algorithmical aspect and benchmarked the model against analytical results and other numerical codes.
For simulations we used our numerical model integrated as a module in a supercomputer 3D particle-in-cell code PICADOR [24]. For performing simulations we used the following setup.
The simulation box of 4 µm in each direction has been used together with a thin PML layer at the boundaries to mimic a part of infinite space by generating the incoming diploe wave and absorbing the outgoing emission. The considered target was a uniform plasma region of spherical shape with the diameter of 3 µm. The initial temperature was adjusted for different initial densities to make the Debye length resolvable by the grid. In all the simulations the initial temperature was below For the numerical experiments with realistic laser pulses we used a Gaussian temporal profile for the incoming dipole wave with a duration of 15 fs (FWHM for intensity). For the considered total peak power of 40 PW the initial density of 10 13 cm −3 provided the absence of self-action within the entire duration of the numerical experiment.

Analytical estimates and calculations
The cascade development within each turning phase is not exponential. However, because the produced particles are trapped by the ART mechanism they generate very similar cascades during the next turning phase, and as a result the number of particles does grows exponentially on timescales longer than a half-period. Let Γ be the average rate over a half-period; simulation results for Γ multiplied by the laser period T are shown in Fig. 5 as a function of total power P. In order to provide some further analytic estimates we model this curve by the solid line shown in the figure, which has the form where P min 7.2 PW is the threshold value required for the cascade to be initiated, and κ = 3.21, both parameters numerically determined.
The equation describing the growth of particle number N(t) is: Assuming that the dipole wave has a Gaussian temporal shape with peak power P 0 and duration τ (FWHM for the power), so that P(t) = P 0 e −4 ln(2)t 2 /τ 2 , we can integrate (6) and a relation between the growth of the particle number from an initial value N 0 to N max at the instance of peak intensity, t = 0. Writing x ≡ P/P min we find where For P in the range P min . . . 100 the Error function may be approximated by its series expansion, but the first three nontrivial orders are needed to give good agreement over the whole range. The asymptotic expansion does not seem to work over this range. By rather crudely fitting points we can approximate, for P above threshold and less than 100 PW , f (P/P min ) ∼ 0.07(P 0.59 − 3.6) , which gives a good approximation over the whole range, the only real deviation being at the threshold, see Fig. 6.
To proceed we need estimates for N 0 and N max . As these appear inside the logarithm in (7) the estimates can be quite rough. We first choose N 0 large enough to be safely sure that some particle generates cascades: N 0 = 10, is sufficient. In order to identify the upper limit N max we equating the Coulomb field at a distance of ρ eff to 1/10 of the peak field of the dipole wave (i.e. 10% self-interaction effects): this yields N max ≈ 10 9 √ P 0 , which as we will see below fits well with simulation results.
τ T 2.23 log[10 8 P 1/2 ] P 0.59 − 3.6 (10) Numerical calculation of the photon yield for a Gaussian pulse To confirm that the proposed concept is optimal for producing a maximal number of photons we use the results of numerical experiments with constant power P to calculate the number of photons that would be generated with an arbitrary Gaussian pulse. We use the following procedure.
All intervals between instances of E = 0 are, despite the complexity of the evolution within each half-cycle, self similar. This is because the particles inherit almost no information from previous half-cycles after loosing energy due to strong magnetic field (when it peaks). Thus, assuming a slow variation of the wave amplitude, we can describe the discreet evolution of the particles number N ee + and the number of generated photons N ph (which escaped from the strong field region) at the instances of E = 0 using the following equations: N ee + (t + T /2) = N ee + (t)Y (P (t + T /4) , N ee + (t)) , N ph (t + T /2) = N ph (t) + N ee + (t) X (P (t + T /4) , N ee + (t)) , where the functions X (P, N ee + ) and Y (P, N ee + ) are defined as the relative increase in the number of particles and the number of generated photons within a single half-period between instances of E = 0, for the current value of power P and number of particles N at the beginning of the interval.
For our study we assume that N ph is the number of photons with energy above 2 GeV. We use the data obtained in the set of experiments with constant power P and interpolate the values in other points. The result is presented on the fig. 7. Using equations (11) and the numerically determined functions X (P, N ee + ) and Y (P, N ee + ) we then can calculate the expected number of photons to be generated by an arbitrary Gaussian pulse. The results are presented on the fig. 2(b).