Scaling laws for direct laser acceleration in a radiation-reaction dominated regime

We study electron acceleration within a sub-critical plasma channel irradiated by an ultra-intense laser pulse ($a_0>100$ or $I>10^{22}~\mathrm{W/cm^2}$). In this regime, radiation reaction significantly alters the electron dynamics. This has an effect not only on the maximum attainable electron energy but also on the phase-matching process between betatron motion and electron oscillations in the laser field. Our study encompasses analytical description, test-particle calculations and 2-dimensional particle-in-cell simulations. We show single-stage electron acceleration to multi-GeV energies within a 0.5 mm-long channel and provide guidelines how to obtain energies beyond 10 GeV using optimal initial configurations. We present the required conditions in a form of explicit analytical scaling laws that can be applied to plan the future electron acceleration experiments.

In the interaction of under-dense plasma with a strong laser pulse, plasma electrons can be accelerated to relativistic velocities [1]. Depending on the laser pulse duration, different mechanisms are responsible for electron acceleration. Using laser-wakefield acceleration (LWFA) tecnique, one can obtain quasi-monenergetic electron energy distribution [2][3][4][5]. The highest electron energy achieved in experiments with LWFA to-date is 7.8 GeV [6]. For an efficient LWFA, the laser pulse duration τ 0 should be on the order of 1/ω p , where ω p is the electron plasma frequency. Short lasers (τ 0 ∼ 30 fs) are required to achieve this.
For a long laser pulse τ 0 1/ω p , where LWFA is not efficient, an alternative acceleration mechanism can be explored. Upon laser propagation through the underdense plasma, the ponderomotive force expels electrons sideways resulting in the formation of a positively charged plasma channel. The formation of longitudinal plasma wave used for particle acceleration in LWFA is prevented by the length of the laser pulse. However, the formation of long-range transverse fields is still possible, due to the charge displacement and the electron current. Electrons propagating at relativistic velocities aligned with the channel axis undergo betatron oscillations, similarly as in LWFA bubble regime. At the same time, they are also experiencing the transverse laser fields. If the frequency of betatron oscillations is close to the Doppler-shifted laser frequency, then the laser energy can be efficiently coupled to the electron (so-called betatron resonance). This regime of electron acceleration is referred to as the direct laser acceleration (DLA) [7]. The rich physics of this setup has raised a keen interest among researchers in the past few years [8][9][10][11][12][13][14][15][16][17][18][19][20][21][22][23][24][25]. It was reported that the efficiency of DLA may be enhanced using various strategies: parametric amplification of betatron oscillations [26][27][28][29], applying additional longitudinal electric field [30,31] or by "breaking of the adiabaticity that precludes electron energy retention" [32]. Even a few DLA experiments have been performed using weakly relativistic laser intensities [33,34].
Ultra-intense laser beams will become available through the next generation of 10 PW-class laser systems [35][36][37]. Interaction with a very intense pulse can provide a strong acceleration through DLA, but the underlying physics becomes even more complex than before. Relativistic electrons in a strong electromagnetic field lose energy by radiation, which in turn alters their trajectories. [38]. This phenomenon is called radiation reaction (RR) and is expected to affect electron motions at laser intensities exceeding 10 22 W/cm 2 . A small local change in electron trajectory can make a difference as to whether the betatron resonance is achieved by that electron or not [29]. It is therefore necessary to assess the effect of radiation reaction on DLA for the lasers of ultra-high intensities, i.e. when a 0 100, where a 0 is the dimensionless normalized vector potential defined by a 0 86 I 0 [10 22 W/cm 2 ] λ 0 [µm], with I 0 and λ 0 representing the laser intensity and wavelength respectively.
The question regarding the role of radiation reaction in DLA was addressed in several recent works. The studies encompass the effect of radiation reaction on the formation of an ultradense helical electron bunch [39,40] and emission of γ-rays [41,42]. Somewhat counter-intuitively, radiation damping can be favourable for particle acceleration. One way to contribute is through allowing particle trapping near the channel axis in situations where all electrons would be expelled if no radiation was emitted. The trapped electrons interact with the peak laser intensity and can potentially achieve multi-GeV energies [43][44][45]. Another way to contribute is through the change of resonant conditions. Even though radiation reaction may reduce the global maximum achievable energy of accelerated electrons, it can strongly enhance acceleration of some electrons [46].
Obtaining a precise description of DLA is a challenge even without considering radiation reaction. By extending the mathematical tools previously developed for classical relativistic DLA at moderate laser intensities (a 0 10) [18], this paper aims to pave the path for a quantitative description of DLA including radiation reaction. In particular, we explore how RR changes the onset of betatron resonance and provide analytical predictions for the final electron energy. We discuss different aspects that affect or limit the acceleration, and find the promising parameter ranges and scaling laws to guide future experiments. The laws are verified with a comprehensive numerical study, involving test particles in ideal conditions, as well as full-scale self-consistent 2D particle-in-cell simulations. The setup considered in this work is a pre-formed plasma channel irradiated by an intense laser pulse that accelerates plasma electrons which are initially at rest (somewhere within the channel or at the channel walls). Our findings show that a fraction of particles that can achieve the betatron resonance is greatly increased by RR, which allows for high charge content in the accelerated electron beam. We quantify which are the necessary initial conditions and expected asymptotic energies for those particles.
This manuscript is organized as follows. In section I, we provide the analytical description of electron dynamics within a plasma channel in the presence of an intense laser field. We consider cases with and without radiation reaction in a simplified field configuration. In section II, test particle simulations are presented, the obtained data are compared with the analytical predictions, and the validity of used approximations is verified. Section III features 2D Particle-In-Cell (PIC) simulations with a pre-formed plasma channel, while the summary of the work is given in section IV.

I. ELECTRON DYNAMICS IN A PLASMA CHANNEL WITH A LASER FIELD
In this section, we analyze the acceleration of a single electron in a symmetric plasma channel irradiated by an ultra-intense linearly polarized laser. Using DLA scheme, the electron can be efficiently accelerated only if it achieves the betatron resonance. Whether this condition is fullfilled depends on initial conditions of the electron (e.g. distance from the axis, initial momentum, etc.) and the background plasma density. In addition, as the electron performs betatron oscillations in a strong electromagnetic background field, it loses its energy by emitting photons which, in turn, alters its dynamics [38]. All these effects have an impact on achieving the betatron resonance, and thus on the overall dynamics of the electrons within the plasma channel.
We first introduce our configuration and an analytical description of a particle performing the coupled oscilla- tions. This motion has an invariant I that can be used to obtain the resonance condition. The next step is to estimate how much energy particles can gain over a certain acceleration distance, still without radiation reaction. We show later on that the principal effect of radiation reaction is to limit to the maximum energy achievable on a resonant trajectory. More specifically, radiation reaction breaks the invariant I and it is possible to quantify this change over a resonant cycle. A decrease in I results in a gradual change of the resonant condition over time, which allows the electron to become resonant even if its initial conditions were far from optimal.

Integral of motion in simplified electromagnetic configuration
A rigorous analytical description of our setup requires a few simplifications. Similarly as in previous works, in this section we consider the laser to be a plane wave, and the channel fields to be a linear function of the distance from the propagation axis x [18,19,27]. Consequently, the channel fields are radially symmetric with respect to the x-axis. The electron can be placed initially at different radial positions inside the plasma channel or on channel walls (see figure 1). We consider electrons starting at rest, and in the (x, y) plane (the plane defined by the laser polarisation and propagation direction). The laser propagates in positive x direction and the fields are given by E L = E 0 sin φŷ, B L = B 0 sin φẑ, where E 0 and B 0 are the amplitudes of the electric and magnetic field and φ is the phase of the wave. The phase velocity of the laser is assumed to be equal to the speed of light c, which is justified by the low plasma density and high laser intensity (this is verified in section II).
The electromagnetic field experienced by the electron is the combination of the laser field and the fields emerging due to the displacement of plasma electrons in the channel (channel fields) [7]. These self-generated quasistatic channel fields are the radial electric field and the azimuthal magnetic field where r = yŷ + zẑ is perpendicular to the channel axis and v = vx is the velocity of the flow. The numerical factor f depends on the fraction of electrons within the plasma channel and takes values between 0 ≤ f ≤ 1 [17-19, 21, 26, 29]. The transversely expelled electrons generate the radial electric field, while electrons accelerated forward within the channel form a current that generates the azimuthal magnetic field. Usually, the higher the background plasma density, the lower the value of f [44]. In other words, the channel fields are linearly dependent on a radial distance from the channel axis, but the electric field E C and the magnetic field B C do not necessarily have the same magnitude. The total electromagnetic field experienced by the electron is then given by The field structure defined in equation (2) induces electron oscillations due to the laser field as well as betatron oscillations at the same time. The background plasma density n p affects the electron motion since the magnitude of the self-generated channel fields is proportional to the plasma frequency: ω 2 p = 4πe 2 n p /m e . We, therefore, expect DLA to be sensitive to the initial conditions of the electron, the intensity of the laser pulse and the density within the plasma channel.
Without radiation reaction, the electron motion in the channel is only governed by the Lorentz force: where γ = 1 + (|p|/m e c) 2 is the relativistic Lorentz factor and p is the electron momentum. For particles propagating in the positive x-direction, both E C and B C contribute in a similar manner: they provide a restoring force that pushes electrons towards the channel axis. In fact, the value of the numerical factor f from equation (1) is not important because the restoring force is actually proportional to |E C | + |B C |. From the Hamiltonian of the electron, one obtains an integral of motion I. For a particle that is initially in the (x, y) plane with z = 0, the integral of motion can be written as [18] where p x is the component of the electron momentum in the direction of wave propagation. For a better intuitive understanding, it is useful to note that the first two terms of the integral of motion are the same as for the particle interacting with a plane wave in vacuum, while the third term accounts for the transverse harmonic oscillations within the channel. The value of I defines a limit on the energy attainable in the system if all dissipation mechanisms are neglected. For example, in the case of the electron initially at rest placed on the channel axis (y = 0), the value of the integral of motion is I = 1 [18].

Betatron resonance
Due to the interaction with the channel fields, the electrons perform betatron oscillations with the frequency given by ω β = ω p / √ 2γ [47]. We define the betatron period to be T β = 2π/ω β . For electrons propagating in the same direction as the laser, the frequency of the quiver motion due to the laser field depends on how fast the electron moves in the x-direction. If the velocity of the electron is v x while the laser propagates at approximately the speed of light c (underdense plasma) the Doppler-shifted laser frequency can be expressed as The electron in the channel can achieve an efficient acceleration if the conditions for the betatron resonance are met, i.e. when the Doppler-shifted laser frequency ω L is close to the betatron frequency [7] ω L ω β .
If we assume that the momentum of an ultra-relativistic electron accelerated in the channel is predominantly pointing in the forward direction (p x p y m e c), then the forward momentum can be expressed [18] as p x γm e c[1 − (1/2) (p y /γm e c) 2 ]. After replacing p x into equation (4), we obtain [25] I 1 2γ Equation (6) allows to obtain an expression for p y m e c √ 2γI cos ψ, where the amplitude of the oscillations is written as y = R = R 0 sin ψ, and ψ is the betatron phase. The Doppler-shifted laser frequency is now The resonant condition ω L ω β applies to an average over the entire cycle, and not for the instantaneous values of ψ. For the condition to be satisfied, the required particle energy is I being an invariant, we can also write I = I 0 and I 0 = 1+[ω p R 0 /(2c)] 2 for a particle initially at rest. This means that the maximum γ * max and average γ * in a resonant cycle can be expressed as using cycle-averaged value cos 4 ψ = 3/8. According to equation (9), the resonant electron energy does not depend on the laser intensity but solely on the initial position and the background plasma density. However, the intensity will determine how long it may take an electron to get accelerated towards the resonant energy. In the next subsection, we will discuss the work of the electromagnetic field for optimal particle acceleration, with a special consideration of resonant electrons.

Optimal acceleration for quasi-resonant electrons
The work performed by the electric field of the laser on one particle is given by dW = −e E · p/(γm e ) dt. In our configuration, the fields are purely transverse and only the p y component is relevant This expression is composed of two oscillatory functions. It is therefore crucial to determine under which conditions the averaged work does not vanish. Intuitively, the initial phase of the particle should be favourable and the oscillation frequencies should be of the same order. This can be understood by looking at figure 2. The field is doing constructive work when p y is anti-parallel with E L . If T L T β , then during one betatron oscillation, there are many oscillations of the laser field whose work in the positive and negative half-cycles would cancel making the total work vanish. By approximating the slowly varying function sin ψ as a step function, then at best we can obtain one entire laser cycle T L of the constructive work within one T β (sections of the cycle with constructive contribution are shaded in figure 2). This would give sin φ cos ψ max | sin φ| T L /T β (2/π)T L /T β . If we apply the same idea, but account on average for the fact that betatron oscillations are not a step function, we get | sin φ cos ψ| max (2/π) 2 T L /T β . The work of the field is therefore expected to be less efficient when there is a large discrepancy between the oscillation periods T L and T β . Inserting our estimates into equation (10) and expressing the values of T L and T β as a function of I, ω 0 , and ω p , we can estimate the maximum attainable constructive power dW dt This represents the power for the best case scenario -there is no guarantee that any particles would meet this condition. Figure 3 illustrates a resonant vs. nonresonant particle trajectory. It is clear that some particles never get accelerated, despite being in the strong field. A strict calculation would be possible for resonant particles for which sin φ cos ψ max 1/2. This yields a solution similar to what equation (11) predicts for T L T β .
Equation (11) considers the local amplitude of the electric field E 0 . As we are usually dealing with laser pulses, this local field will vary in space and time. We could assume, for example, a sin 2 shape function for the laser temporal envelope. This would lower the total work done by the laser from the beginning to the time when the particle gets to the highest field region approximately by a factor of two. From that we can derive the expected cutoff energy versus time which can also be presented as a function of the acceleration distance l acc From (12) and (13), we observe that particles which start with a low I can experience a more efficient acceleration. Moreover, their resonant γ is lower according to equation (9), so the resonance may be easier to achieve. Consequently, for any setup with a predetermined acceleration distance, we may have particles at a preferential initial I to achieve the betatron resonance at a highest value of γ. We can then extract the value for the optimum initial distance from the channel axis The particles that were initially at the distance R opt as defined by equation (14) have the best conditions for achieving the betatron resonance and getting accelerated. This remains true as long as the interaction is without radiation losses, i.e. the integral of motion is unchanged I I 0 . We draw attention to the reader that the above calculations did not account for a transverse shape of the laser pulse, as the derivation considers a plane wave laser with a temporal envelope. However, one could extend this model by using a lower effective a 0 to account for the transverse envelope function.

Effects of Radiation reaction
Radiation reaction must be present in addition to the Lorentz force to accurately resolve the motion of a particle when ultra-intense laser beams are considered. We model the radiation reaction term using the Landau-Lifshitz (LL) equation [38] dp where the leading order term of the radiation reaction force is given by (16) We mentioned before that the electron dynamics without radiation reaction does not depend on the value of f in equation (1) [21]. We will now show that even the radiation reaction force at the first order in I/γ does not depend on the value of f . For simplicity, we assume again the laser to be a plane wave that is linearly polarized along the y-axis, and propagates along x. We consider an electron that starts in the plane z = 0 (in fact in the sub-plane y > 0 without loss of generality).
The electromagnetic fields experienced by the electron are given by equation (2), where E C and B C are positive by construction. The laser fields E L and B L can take any sign, but they are equal to one another in magnitude: E L = B L = E 0 sin φ. This field configuration cannot change p z if it was initially zero, so the particle is confined in the (x, y) plane. We can simplify the term in square brackets of equation (16) considering the directions of E and B from equation (2) and replace p x and p y according to the previous values where We can draw an important conclusion from equation (17). In an extremely intense laser, it is likely to have E * C E 0 . However, even though E * C has a significantly smaller magnitude, the channel field dominates the radiation reaction as long as E 0 E * C γ/I. This inequality can be violated for a large transverse momentum component p y p x . It is not excluded that some particles could be in this situation before getting accelerated in the positive x-direction.
The energy loss can be approximated by is the damping constant [38]. Radiation reaction sets a limit on the maximum electron energy attainable in the interaction with an ultra-intense laser pulse. The limit is reached when all the energy an electron can gain during one betatron cycle is radiated out during the same time interval. As a result, the particle energy is equal before and after such a cycle. In other words, no further acceleration is possible. For a relativistic particle that satisfies the condition E 0 E * C γ/I, we can estimate the maximum energy by equating dγ/dt from equation (10) and equation (19) which leads to Equation (20) represents the acceleration limit when the channel field is the dominant cause of radiation reaction and can be applied to quasi-resonant particles. The result can be simplified to Radiation reaction-induced reduction of I When electrons lose energy due to radiation emission, the integral of motion I defined in equation (4) is not conserved anymore. We can make a quantitative estimate of how I changes during one resonant betatron cycle. Radiation losses change the instantaneous energy and momentum of the particles, which also reduces the integral of motion I. For relativistic particles, we can assume that most of the radiation is emitted in the direction of motion (this emission is contained within a cone that has an opening angle θ ∼ 1/γ). According to this assumption, the momentum and energy are reduced proportionally: dγ/γ dp/p dp x /p x . The integral of motion therefore decreases through the reduction of γ − p x /(m e c), in the following way: where we used the proportionality mentioned above to express dp x as a function of dγ. Using equation (22) we can derive that the integral of motion is reducing at a rate where dγ/dt is given by equation (19). We can estimate how much I decreases during one resonant betatron cycle. We first perform a change of variables dt = γ * /(I cos 2 ψ) dψ, and use equation (8) to remove direct γ dependence. If we assume a limit where the channel field dominates, using 2π 0 cos 8 ψ sin 2 ψ dψ = 7π/128 we get I decreases to In the other limit, where the laser field dominates, we have 2π 0 cos 6 ψ dψ = 5π/8, which gives Both equations (24) and (25) predict an asymptotic value that particles with a large initial value of I 0 would tend to after a resonant cycle. If the predicted asymptotic I is much smaller for one of the limits (RR due to the laser or due to the channel fields), it means that this limit represents the dominant contribution to the radiation reaction. As it turns out, the I limit predicted by equation (24) is frequently lower than the one predicted by equation (25). An example for a 0 = 500 is given in figure 4. The figure illustrates that even one resonant betatron cycle can be enough to reduce the integral of motion very close to the value where the radiation reaction is not very strong anymore and cannot make a significant change during the subsequent cycles. Inserting the asymptotic value of I into equation (21) will give us the maximum electron energy allowed due to the radiation reaction.
For predictions of the maximum energy in the system, one should compare equations (9), (13) and (21), using the lowest energy predicted among the three as a final result. One should consider the I 0 for equation (13), because it applies to acceleration from the beginning of the interaction, while equations (9) and (21) should be considered with the value of I reduced by radiation reaction.

II. TEST PARTICLE SIMULATIONS, COMPARISONS WITH THE ANALYTICAL MODEL
To validate our model and probe parameters optimal for particle acceleration, we now use a test-particle approach in a simplified electromagnetic configuration, similarly as in section I. We consider a wide range of initial conditions varying the particle distance from the channel axis, background plasma density and the laser intensity.  (25) for a laser with a0 = 500, while the coloured lines represent the predictions of equation (24) for different values of ωp/ω0. For reference, I = I0 is given as a solid black line.
This data can illustrate how limits presented in the previous section affect particle acceleration.
The laser is represented as a plane wave with a temporal envelope [21,29], which propagates along the positive x direction at the speed of light and is polarised along the y-axis. We consider λ 0 = 1 µm and a Gaussian temporal profile with duration τ = 150 fs defined as full-width-at-half-maximum (FWHM) in the electromagnetic field amplitude. Due to the slowly-varying laser envelope function, the initial carrier-envelope phase does not impact the laser-electron interaction. The test particle equations of motion are integrated using the 4th order Runge-Kutta method with a time step ∆t = λ 0 /(1000c) satisfying the condition for accurate calculation of electron dynamics in a strong laser field [48]. The numerical factor f introduced in equation (1) is here set to 0.5, such that the electric and magnetic field of the channel have an equal magnitude. The plasma channel is assumed to be axially symmetric, with a typical field structure given by equation (1). In our calculations, the length of the plasma channel is 528 µm. In the following analysis, the distances are often expressed in units normalized to the laser frequency where c/ω 0 ≈ 0.16 µm. The electrons always start at rest at an initial distance R 0 from the channel axis. This corresponds to a case when an intense laser propagates through a pre-formed plasma channel and accelerates electrons originating from the channel walls. We vary R 0 from 0 to 100 c/ω 0 for different values of ω p /ω 0 ranging from 0.01 to 1. For an optical laser, this corresponds to a plasma density range between 10 17 and 10 21 parts./cm 3 , and distance from the channel axis up to ∼ 15 µm. Radius and ω p /ω 0 are sampled with 100 values each, which makes for 10 4 test cases just by varying these two parameters. We also consider several values of a 0 , with examples of strong and weak radiation reaction.
Without radiation reaction Figure 5 summarizes the results obtained from testparticle simulations without radiation reaction for three different values of a 0 . Each data point represents the maximum achieved value of relativistic Lorentz factor γ of an electron as a function of the initial value for I 0 ω 0 /ω p . Some particles achieve energies over 10 GeV, and many of them are resonant according to equation (8). This is verified in figure 5 where γ * max is shown as a black solid line, while cycle-averaged value γ * is represented with a dashed line. Particles are resonant at different energies for different values of I 0 ω 0 /ω p , and, as shown in section I, the resonance condition does not depend on the laser intensity. The asymptotic energy, therefore, has to be estimated by considering how much work towards particle acceleration the laser can invest during the available acceleration distance. According to equation (13), this limit is a function of a 0 , I 0 and ω p . The color-coded solid lines show the maximum possible electron energy for quasi-resonant particles given by equation (13) corresponding to the acceleration distance of 528 µm. For a specific set of parameters a 0 and ω p , the most favourable initial R 0 for acceleration is the one that allows to achieve the betatron resonance at a highest electron energy within the available acceleration distance. This optimal R opt corresponds to the intersection between the dashed and the coloured solid curves in figure 5, and is given by equation (14). The acceleration is in principle possible even far from betatron resonance, with a slightly lower maximum allowed energy. However, as we have mentioned before, the probability for constructive work performed by the laser towards particle acceleration is low for high values of I 0 , where R 0 R opt . The existence of an optimal initial radial position is better illustrated in figure 6, where panels (a) and (b) display the test electron maximum energy as a function of R 0 and ω p /ω 0 for the same dataset summarized in figure 5. We can also see sudden jumps in the maximum achieved energy between points with similar parameters, which suggest a resonant nature of the acceleration mechanism, as highlighted in the previous works [7,8].

Effects of radiation reaction
When considering extreme laser intensities, we expect radiation reaction to be important for the particle dynamics, and consequently for acceleration. This section addresses the radiation reaction effects. We present test particle simulations for the same parameters as in the previous section, but now including RR in the equations of motion. In the code, RR is implemented as a continuous process of losing electron energy described by the Landau-Lifshitz equation [38,49]. We do not include the electron Gaunt factor correction in these simulationswe verify later in this section that this does not change our findings. The results of test-particle simulations for a 0 = 100 and a 0 = 500 are shown in figure 6 (c) and (d).
The differences between cases with and without RR for a 0 = 100 are small, both in the absolute values of the energy obtained as well as in the parameters favourable to achieve betatron resonance. This suggests that RR is not strong for a 0 = 100. In contrast, the maximum energy achieved with RR is about a factor of two lower than without RR for a 0 = 500. But the most striking difference is in the parameter range where the maximum energy is achieved. It appears that all particles with R 0 > R opt can get accelerated. Since radiation reaction affects the electron dynamics and thus the evolution of the phase-matching process between the particle oscillations in the channel and the laser field, an efficient electron acceleration can be achieved for a wider range of initial conditions. This is accomplished through the reduction of the integral of motion, which allows particles with an initially large I 0 to eventually converge to a lower quasi-resonant value of I that enables energy retention.
We have compared the theoretical estimate for reduction of I during one resonant cycle with the data ob-tained from test-particle simulations. The predictions of equation (24) against the data for a 0 = 500 are shown in figure 7. The values of I for our parameter range are restricted by three boundaries, as shown in figure 7 (a). The solid line represents a situation of weak radiation reaction where the integral of motion does not change. The dotted line corresponds to the highest plasma density we considered (ω p /ω 0 = 1), while the dashed line is for electrons starting at the furthest initial position from the channel axes. Figure 7 (b) shows the integral of motion calculated from the data at the moment when the electron reaches its maximum energy. Equation (24) is a good predictor for the order of the final value of I, even if we consider only one resonant betatron cycle. Some particles will have enough time to perform multiple betatron oscillations, and achieve the asymptotic value of I where the radiation reaction becomes negligible. However, since this value is of the same order, for simplicity we will continue considering a single betatron cycle for computing the final I. After inserting this corrected value of I into equations (9) and (21), the equation that predicts a lower energy is a good estimate of the asymptotic energy of the electron, provided that the acceleration distance was long enough to achieve this energy according to equation (13).
We can verify this by presenting the test particle simulations with radiation reaction in a similar fashion as in figure 5. Figure 8 shows the maximum particle energy as a function of I 0 ω 0 /ω p for two different laser intensities, a 0 = 400 and a 0 = 600. If we assume that the electron is initially at the maximum distance R max = 100 c/ω 0 from the channel axis, then equation (21) provides the limit on the attainable electron energy which is illustrated by the dashed line in figure 8. The dotted line represents the expected value for ω p = ω 0 according to equation (21). For ω p /ω 0 < 0.3, there is a window of parameters that allows energies above 10 GeV. However, we do not see them in our results, because of the finite interaction time. As has been stated in the previous subsection, longer acceleration distance allows to accumulate more energy. This is illustrated in figure 9, where all parameters are kept equal, except the acceleration distance. The upper panels show the same plots as figure 7 for a 0 = 500 and l acc = 528 µm or l acc = 1056 µm. The lower panels show the obtained energy vs. expected energy according to the minimum predicted by equations (9) and (21). There is a clear convergence of the data towards the maximum allowed energy we predict, but there are also outliers, especially for ω p /ω 0 < 0.3. These particles do not attain their allowed maximum, either because there is not enough time, or because their initial integral of motion is already so low that the radiation reaction does not significantly affect their motion and they can easily remain far from the betatron resonance. A fraction of these particles eventually does converge towards the predicted values in the example with a longer propagation distance.  (b) Comparison with data points from test-particle simulations at the moment when the electron reaches the maximum energy. The considered laser intensity for both panels corresponds to a0 = 500.

Approximations: limits and validity
We have introduced several approximations in our analytical model, that we hereby justify.
The first question that might arise is that of the choice of numerical factor f in equation (1)  eter study for a 0 = 500 in figure 10 (a) and (b). There are differences for individual test particles (represented by individual pixels), but the expected maximum energy for a given region of parameter space remains unchanged. However, if the electric field turned out to be very low (i.e. f → 0), particles would need an initial transverse momentum to be efficiently accelerated, because there would be no electrostatic potential in the y-direction.
For simplicity, we have neglected the semi-classical correction to the radiation reaction, modelled by the Gaunt factor g(χ e ) that reduces the amount of emitted power as a function of parameter χ e . The parameter χ e characterizes the interaction of an electron with a strong electromagnetic field [50,51], and it is approximately equal to the ratio of the field as seen by the electron in its rest frame to the Schwinger limit defined by E S = m 2 e c 3 /( e). Our analysis considered only the leading order term of the radiation reaction, while the correction function is the next order in χ e [52]. The final electron energies within the channel depend on where the radiation reaction shuts down, and therefore, the most important section of the trajectory should be in the regime where χ e 1. In that case, the g-factor is a very small correction. We have verified this by introducing the g-factor and showing the expected electron energies on a parameter study for a 0 = 500 in figure 10 (c). The differences from results in panel (a) are very small, which confirms this approximation was adequate.
The phase velocity of the laser propagating through plasma is considered luminal in our model. The actual phase velocity can be expressed as v ph c [1 + n p /(n c a 0 )], where n c = ω 2 0 m e /(4πe 2 ) represents the non-relativistic critical plasma density. This evidently becomes important for dense plasmas and moderate laser intensities. However, as phase difference accumulation between the particle and the wave is an important factor in the acceleration process, there is a possibility that even a small change in the phase velocity of the laser may be relevant for particle dynamics. The luminal approximation is justified as long as the phase accumulation due to the transverse motion of the particle is much larger than the effect of super-luminosity. This is equivalent to stating that the difference between the speed of light and particle average velocity in the direction of laser propagation c − v x is much larger than v ph − c. We can estimate using p x from section I, that v x c 1 − (1/2)( p y /(γm e c)) 2 which can be approximated as v x c [1 − (1/2)(I/γ)]. To neglect superluminosity, the following condition must be met I/(2γ) (ω p /ω 0 ) 2 /a 0 . As I reduces, and γ increases over time, the left-hand side of the inequality has a lowest value in the asymptotic conditions, where the particles have reached their maximum allowed energy. We can then use the expression for γ * from equation (9), as well as I ∼ 50 ω p /ω 0 according to equation (24). The condition that must be satisfied then becomes a 2 0 5×10 3 (ω p /ω 0 ). The condition is verified for our parameters, which can also be confirmed from figure 10 (d). Therefore, similarly as in reference [46] we conclude superluminosity is a minor effect in our conditions. According to this estimate, for optical lasers a 0 300, the superluminal phase velocity would become important for targets of solid density n p 100 n c , and negligible for gas jets where n p ∼ n c .

III. PARTICLE-IN-CELL SIMULATIONS
We also compare the maximum electron energy predicted by theory against 2D PIC simulations performed with OSIRIS [53,54]. We consider several examples of initial conditions for R 0 and ω p /ω 0 , and two values for the laser intensity a 0 = 400 and a 0 = 600. The simu-  figure 11. Such a channel can be generated by a preceding laser pulse, as presented in reference [44]. The existence of a plasma channel is beneficial for laser guiding and can allow for longer propagation distances than if we rely on self-focusing.
Plasma density at the channel axis is n p . At the channel radius R 0 , the plasma density reaches its peak 4n c , while outside the channel it is set to 2n c . The protons are mobile and have the same initial density profile as the electrons. The total length of the plasma slab is l acc = 528 µm, and there is a 5.5 µm ramp at the entrance and the exit of the plasma. We use the moving window, and periodic transverse boundary conditions. The laser duration is 150 fs, defined at FWHM in the electromagentic field amplitude. The transverse profile is Gaussian, with a focal spot radius W 0 = 3.2 µm, and the laser wavelength is λ 0 = 1 µm. Any effects potentially associated with super-luminal phase velocity of the laser are naturally included in the PIC simulations. We will show this later. The radiation reaction is incorporated as follows: for χ e ≤ 0.2 we use the Landau-Lifshitz equation [38] with implementation details given in reference [49]. For χ e > 0.2 it is modeled with a QED Monte-Carlo based algorithm [55,56]. The correction to the classical radi-ation reaction due to the electron Gaunt factor was not considered in the PIC simulations, as it is not expected to affect our conclusions according to Section II.

Energy cutoff
In the following, we compare the highest energy achieved in PIC simulations with the scalings presented in section I. In one PIC simulation, parameters a 0 and ω p /ω 0 are fixed, but particles can have different R 0 , and the maximum available integral of motion I 0 will depend on the channel radius and laser spotsize. According to sections I and II, the particles with the maximum initial I 0 will have the best conditions for acceleration. We therefore use the channel radius to calculate the characteristic I 0 which gives us the estimate for the expected cutoff energy in a PIC simulation. Figure 12 shows how maximum energy in the simulation box grows as a function of the propagation distance, both with and without radiation reaction. If radiation reaction is neglected, the electron energy should scale according to equation (13). This trend is confirmed by the results from 2D PIC simulations. However, radiation reaction limits the energy gain. Predictions of equations (9) and (21) are shown together with the data from the 2D PIC simulation with radiation reaction that lies below At the distance R0, the plasma density reaches its peak 4nc, while outside the channel it is set to 2nc, where nc is the critical plasma density. The incoming laser pulse is propagating along the positive x-direction. all limits. As expected, the energy cutoff is lower with radiation reaction but the number of accelerated particles is higher. For example, at t = 210 T 0 , shortly after the radiation reaction effects become visible in global quantities, the simulation with radiation reaction has 2.5 times more particles with energy above 2 GeV than the simulation without radiation reaction, even though the latter has a cutoff energy nearly twice as high.
Among the analytical limits, equation (9) predicts the lowest expected value for the considered parameters. According to our previous analysis, we would expect in this case to have an energy cutoff given by γ * max . Instead, the data with radiation reaction converges towards the cycle-averaged value γ * , which is lower. This is not surprising, because our previous analysis considered ideal conditions: the laser transverse intensity distribution was neglected, as well as the energy depletion during propagation. In a realistic setup, particles do not feel the maximum laser intensity everywhere in the channel, and they do not necessarily originate exactly from the edges of the channel. The channel fields are a linear function of R just in a region close to the channel axis, and we have no guarantee that particles initially further than a certain threshold radius would be accelerated because of the finite transverse laser profile. The asymptotic I can also be lower than the prediction of equation (24) applied to the I 0 that corresponds to the channel radius. All these effects cannot be formally taken into account for analytical predictions. But if we probe several different combinations of parameters, taking γ * instead of γ * max provides a very good agreement with the energy cutoff measured from 2D PIC simulations. This is shown in figure 13 (a). The energy distribution of the electrons is presented in figure 13 (b). The first two examples illustrate how channel radius changes the electron distribution function. They have the same laser intensity FIG. 12. The maximum electron energy as a function of the propagation distance for the case without (blue) and with (orange) radiation reaction for ωp = 0.5 ω0, R0 = 47 c/ω0 and a0 = 600. Solid blue line represents the analytical estimate given by equation (13) without radiation reaction (considering I0). Orange lines are obtained using equations (9) and (21) with radiation reaction (considering I given by equation (24)). Crosses show the corresponding results from 2D PIC simulations. a 0 = 600 and the same background plasma density for which ω p /ω 0 = 0.5, but the channel radius changes by a factor of two. The energy cutoff does not change significantly, which suggests that the narrower channel is wide enough for electrons to achieve the maximum allowed energy. However, the number of particles that approach this value is significantly higher in the case with the wider channel, which is reflected on the slope of the spectrum. Other examples show that the cutoff can be increased by reducing the plasma density. The optimal initial distance from the axis depends on the plasma density according to equation (14), so reducing the density to ω p /ω 0 = 0.32 does not increase much the energy cutoff if the channel radius remains small R 0 = 47 c/ω 0 . However, if one combines the wider channel with the lower plasma density, the electron energy cutoff evidently increases. By applying this strategy, it is possible to obtain more energetic electrons even using lower laser intensities (e.g. a 0 = 400 in figure 13).

The electromagnetic field structure
We proceed to a more detailed analysis of PIC results, to illustrate the connection with the analytical model presented in Section I. In figure 14 we present density and field profiles for a simulation with parameters ω p = 0.27 ω 0 , R 0 = 81 c/ω 0 and a 0 = 600 at time . The detailed structure of the electromagnetic field within the smaller box marked in (e) is presented in panels (b), (d) and (f). Inside the channel, the laser is partially absorbed, but the wavefronts are nearly straight, as can be verified by observing the transverse electromagnetic field components E y and B z . The longitudinal electric field E x generated by the density fluctuations within the channel is noisy, and can potentially work towards or against acceleration of individual particles. Panels (g) and (h) show the total focusing field associated with the plasma channel, given by equation (1). As these fields are weaker than the laser field, one can obtain only the total value |E C | + |B C | by subtracting E y − B z . Doing so reveals a long-range field structure, linearly dependent on the distance from the axis. The amplitude of the focusing field is of the same order as E x , with a difference that the E x average is close to zero on scales larger than a laser wavelength. The color-coded lineouts in panel (h) confirm that the channel field can be regarded as a linear function of the distance from the axis up to approximately R 0 /2, where the slope is defined by the effective plasma density. For a flat density plasma slab, the effective density is the background plasma density n p . For a channel with a parabolic density gradient like ours, the average plasma density for the region around the channel axis between −R 0 /2 and R 0 /2 can be calculated analytically and directly applied to equations presented in Section I. The channel field predicted by equation (1)  Laser-electron dephasing depends on the velocity of the particle along the laser propagation direction β x and phase velocity of the laser. One can measure the phase velocity in PIC simulations from the temporal evolution of the electric field along the channel axis, shown in panel (g). The phase velocity is not constant, owing to the 2D laser dynamics. The measured average phase velocity is 1.00075c, on the same order as the analytical prediction in Section II that gives the value 1.00066c for our n eff p 0.4 n c and a 0 = 600. The β x depends on the ratio p y /p x ∼ 0.1 (c. f. figure 15), which confirms that for resonant particles the transverse electron motion has about an order of magnitude stronger effect on dephasing than superluminosity (and more than one order of magnitude for non-resonant particles). We can therefore apply our analytical model for studying the particle motion. Figure 15 shows the particles perform a similar asymptotic motion at late times, despite the different initial dynamics. Some particles achieve resonance faster than others. The interaction with the laser can start from rest, or after pre-acceleration at the laser-plasma interface. The integral of motion associated with the initial distance R 0 /2 from the axis is I 0 ∼ 160. According to equation (24), the expected integral of motion I after one full resonant cycle would be I ∼ 60, while after two cycles I ∼ 50. Among our selected particles, there are many that completed two full resonant cycles. The temporal evolution of trajectories in p y -energy phase-space is also shown, with a highlighted single-particle example in panel (d). This illustrates how the most energetic particles gain energy over time: they first oscillate in the region of low energy, unable to retain the energy gained over the oscillation cycle; as the oscillations become quasi-resonant, the particles can retain a fraction of the energy obtained, and they accelerate sloshing back and forth from the points on the parabola defined by their current integral of motion through p y = m e c √ 2γI. The definition of the parabola follows from equation (6) when cos ψ 1. The phasespace associated with simulation particles initially at rest within the channel between 680 < x * ω 0 /c < 725 is shown in panels (f) and (h) for times t = 2501 ω −1 0 and t = 4000 ω −1 0 , respectively. The colorscale is logarithmic, to highlight all possible combinations of p y and energy. The widest parabola is associated with the maximum initial integral of motion I 0 = 160 available among the considered particles. As the time passes, the invariant I of some particles can be reduced. The two inner parabolas show the asymptotic value after one (dashed line) and two (solid line) resonant cycles. The energetic particles are eventually expected inside the parabola associated with the asymptotic I, as shown in panel (h). According to equation (20), the expected asymptotic value of the Lorentz factor for I = 50 is γ ∼ 1.1 × 10 4 . We can also predict the amplitude of oscillations around the channel axis for these particles is about 22 c/ω 0 by inserting p y = 0 in equation (6). All these estimates are in agreement with our simulation results.

IV. CONCLUSIONS
The betatron resonance in a plasma channel can result in efficient electron acceleration, but strongly depends on the initial particle position, channel density, laser intensity and acceleration distance. Particles with an initial radius R 0 R opt , where R opt is defined by equation (14), have the best chance to attain the maximum energy in the system. At very high laser intensities, radiation reaction causes a decrease in the maximum electron energy, however, it allows accelerating a higher number of particles as it affects the phase-matching process between betatron and laser oscillations. In fact, nearly all particles with R 0 > R opt have a very good chance to be accelerated due to radiation reaction, provided that R opt is still within the range where channel fields are nearly a linear function of the distance from the laser propagation axis. To that end, it may be optimal to use plasma channel radius ∼ 2R opt .
Analytical predictions of the maximum electron energy as a function of channel density and electron initial radial position are provided by equations (9) and (21) using radiation-reaction reduced value of I given by the lowest estimate from equations (24) and (25). The required acceleration distance to achieve this energy can be estimated through equation (13). Our findings are confirmed by numerical integration of test particle trajectories, as well as 2D full-scale PIC simulations. According to the presented calculations, using near-future lasers (10 PW-class, 150 fs) we can obtain multi-GeV electrons in a single-stage acceleration within a 0.5 mmlong plasma channel. We envisage a roadmap towards electron energies in excess of 10 GeV if the channel and laser parameters are matched for an optimal outcome.
Our findings can be extended to acceleration of externally injected beams, with an arbitrary beam loading (i.e. high current within the channel, similar as in reference [46]). In that case, the background magnetic field due to the current within the channel may be stronger than expected from the background plasma density. For such a setup, one should estimate an effective (higher) value for the plasma density n eff p that would correspond to a comparable current density and the same order of magnitude for the self-consistent background magnetic field within the plasma. Our scaling laws can then be applied using this effective n eff p instead of n p to calculate the asymptotic cutoff energy of the electron beam. An example for a 0 = 200 and a M B = 10 −2 a 0 from reference [46], can be mapped using I ∼ 40 and ω p /ω 0 ∼ 0.3. Our model predicsts the asymptotic electron energy is ∼ 7 GeV, which is in agreement with their simulations that show 7.5 GeV.    This trajectory is overlaid with the py-energy phase-space (blue scale) and analytical estimates for the connection between py and energy when betatron phase is cos ψ 1 (the point of local maximum for py) in panels (f) t = 2501 ω −1 0 and (h) t = 4000 ω −1 0 . The outer solid line is associated with the highest initial integral of motion in the system, while the inner lines show the expected values after one (dotted line) and two (solid line) resonant cycles according to equation (24). Panel (g) shows time evolution of the transverse electric field along the x-axis inside a simulation window moving at the speed of light. The figure represents a magnified portion of space slightly over three wavelengths long, such that phase velocity can be measured directly.