Intrinsic energy spread and bunch length growth in plasma-based accelerators due to betatron motion

Plasma-based accelerators (PBAs), having demonstrated the production of GeV electron beams in only centimetre scales, offer a path towards a new generation of highly compact and cost-effective particle accelerators. However, achieving the required beam quality, particularly on the energy spread for applications such as free-electron lasers, remains a challenge. Here we investigate fundamental sources of energy spread and bunch length in PBAs which arise from the betatron motion of beam electrons. We present an analytical theory, validated against particle-in-cell simulations, which accurately describes these phenomena. Significant impact on the beam quality is predicted for certain configurations, explaining previously observed limitations on the achievable bunch length and energy spread. Guidelines for mitigating these contributions towards high-quality beams are deduced.


Introduction
In plasma-based accelerators (PBAs), an intense laser pulse 1 or high-energy charged particle beam 2 drives a plasma wake sustaining accelerating fields orders of magnitude higher than those achievable with conventional radiofrequency technology 3 , offering a path towards highly compact and cost-effective accelerators. Having reached GeV energies in only centimeter scales [4][5][6][7][8] , femtosecond-long electron bunches with kiloampere current 9, 10 and micron-level emittance 11,12 , PBAs are becoming increasingly attractive for applications such as free-electron lasers (FELs) 13 or future linear colliders 14 , which would strongly benefit from reduced size and cost. In particular, there is also a special interest in the production of sub-femtosecond bunches [15][16][17] to generate short X-ray pulses for ultrafast science 18 .
However, despite these major advances, achieving an energy spread below the percent level has remained an issue. This is a problem not only for applications such as FELs, which require a relative energy spread 10 −319 , but also for the beam transport after the plasma, which could cause a dramatic emittance growth 20 and further reduce the beam quality, limiting any kind of multi-stage acceleration. Although the impact on the emittance can be mitigated by tailored plasma-to-vacuum transitions 21 or the use of active plasma lenses instead of conventional quadrupoles 22 , improving the beam parameters in the plasma is essential for demanding applications. Major efforts have therefore been dedicated towards this objective, including the development of controlled injection techniques [23][24][25][26][27] as well as concepts for mitigating the correlated energy spread arising from the steep slope of the accelerating fields. These include beam loading the wake [28][29][30] , either by the accelerated bunch itself 31 or the injection of a secondary one 32 , the use of modulated 33 or tailored 34 plasma profiles, the insertion of a magnetic chicane in a multistage PBA 35 , or taking advantage of the beam-induced wakefields [36][37][38][39][40] .
Here we describe additional sources of energy spread in PBAs which should also be mitigated, as well as limitations on the minimum achievable bunch length. In particular, we show that the transverse oscillations of beam electrons (known as betatron motion), induced by the strong focusing fields in the plasma wake, can be a significant source of uncorrelated (or slice) energy spread as well as increased bunch length due to path length differences between particles with different oscillation amplitude. Although these effects were already known 41,42 , no model exists to date for their impact on the beam parameters. We present here a novel analytical theory, validated against numerical simulations with the particle-in-cell (PIC) codes OSIRIS 43 and HiPACE 44 , which accurately describes these phenomena in the assumption of relativistic electrons in a non-evolving wake. This model allows us to understand previously observed limitations, such as a finite energy spread even when the initial bunch length approaches zero 45 , and define guidelines for minimizing the impact of these phenomena for high-quality electron beams.

Results
Development of a single-particle model As a first step, a single-particle model capable of describing the coupling of the transverse oscillations into the longitudinal motion is needed. Then, from this model, statistical averages over a particle beam can be analytically obtained in order to determine the impact of this coupling on parameters such as the energy spread and bunch length.
For a relativistic particle moving close to the speed of light c, the longitudinal component of its velocity v| is the magnitude of the particle velocity v v v, is limited by c and can be decreased due to motion in the transverse planes, i.e, if v x , v y = 0. Thus, in the speed-of-light frame ξ = z − ct, where z is the longitudinal coordinate in the laboratory frame and t is the time, a particle performing betatron oscillations will experience a slippage towards the back which, at a time t, will be given by where ξ 0 is its initial position and v ξ = v z − c is its longitudinal velocity in the speed of light frame.
Determining the particle slippage thus requires obtaining an expression for v ξ (t), which, in turn, depends on the evolution of the transverse components v x and v y . In order to find these expressions, the equations of motion of a relativistic electron in a plasma wakefield have to be solved. These are given byṗ p p = −eW W W , where p p p = mγv v v is the particle momentum, m is the electron mass, γ = 1/ 1 − (v v v/c) 2 is the relativistic Lorentz factor of the particle, e is the elementary charge and W W W = (E x − cB y , E y + cB x , E z ) is the wakefield. In this expression, E i and B i (for i = x, y, z) are the different components of the electric and magnetic fields. Analytical expressions of the wakefields can be found for the linear regime 46,47 , and several models have been developed for non-linear wakes in the blowout regime 48,49 . In what follows, the wakefield gradients around the particle, K x = ∂ x W x , K y = ∂ y W y and E z = ∂ z W z , will be assumed constant and ∂ i W j = 0 for i = j. This is especially well-suited for the blowout regime, where the driver expels all plasma electrons and leaves behind an ion cavity with uniform focusing fields K x = K y = mω 2 p /2e and an approximately constant E z in a wide range of the accelerating phase. The plasma frequency is defined as ω p = n p e 2 /mε 0 , where ε 0 is the vacuum permittivity and n p is the unperturbed plasma density. For linear wakes, although these gradients are not uniform, this model can be applied for regions sufficiently close to the propagation axis, where K x and K y can be regarded as constant, and if ∆ξ c/ω p , so that the longitudinal change in K x , K y and E z can be neglected. Under these conditions, the equations of motion for a beam electron in a non-evolving wake with a phase velocity v w given by the propagation velocity of the driver, can be written aṡ where E(t) = −eE z (t)/mc, ω x (t) = K x /γ(t) is the betatron frequency and K x = eK x /m. Since E = −eE z /mc is assumed constant, then , is the energy evolution assuming constant acceleration and γ (1) dt accounts for the contribution of slippage and dephasing with respect to the wake 50 . If γ (1) /γ (0) 1, then γ γ (0) can be assumed in order to solve Eq. (2), for which analytical solutions can be found if ω x (t) is a slowly varying function 51 , i.e, ifω x /ω 2 x = E 0 /2 √ γK x 1. This yields the following expressions: where Γ(t) = γ (0) (t)/γ 0 = 1 + E 0 t/γ 0 is the relative energy evolution, A x,0 = x 2 0 + (v x,0 /ω x,0 ) 2 is the initial oscillation amplitude in the x plane and φ x,0 = − arctan (v x,0 /x 0 ω x,0 ) is the initial phase. The values x 0 , v x,0 and ω x,0 correspond to the initial transverse position, velocity and betatron frequency. The phase advance φ x (t) = t 0 ω x (t ) dt is given by Analogous expressions can be found for the y plane and are not included here for simplicity. The solutions to the transverse electron motion can now be used to obtain the longitudinal slippage, 2|v v v| and |v v v| − c −c/2γ 2 for a relativistic electron. Since the time scale of the oscillation amplitude damping, Γ −3/4 (t), in Eq. (4) is much longer than the betatron period, the time-average of this expression can be used for the electron transverse velocity. This results in the expression

2/13
where A 2 0 = A 2 x,0 + A 2 y,0 and the focusing gradient is assumed to be the same in both planes, i.e K x = K y = K. The two terms in Eq. (6) account, respectively, for the slippage due to |v v v| < c and the slippage caused by betatron motion. Most of ∆ξ occurs initially (at lower energies) and reaches a limit ∆ξ max as Γ → ∞ given by Using Eq. (6) it is now possible to calculate γ (1) and finally obtain γ(t) = γ (0) (t) + γ (1) (t), which takes into account the influence of the particle slippage on its energy. This yields where the third term is a correction to the linear energy gain, which is modified due to slippage, and the fourth one accounts for the influence of dephasing with the wake when v w < c in laser-driven cases. The following terms are higher order corrections which account, respectively, for slippage due to |v v v| < c and slippage due to betatron motion. Eqs. (3), (4), (6) and (8) allow for a complete description of the single-particle evolution within the wakefields. In order to test their validity, they have been compared with numerical solutions of Eqs. (1) and (2). Three test cases corresponding to single electrons with different initial transverse offsets are shown in Fig. 1. The electrons in these three cases (C a , C b and C c ) have, respectively, x 0 = y 0 = 1, 3, 5 µm and v x,0 = v y,0 = 0. They are injected with γ 0 = 100 and propagate 5 cm within a plasma stage with n p = 10 17 cm −3 assuming a blowout with E = −ω 2 p /2c, E 0 = ω p and K = ω 2 p /2. The wake velocity is directly determined from the group velocity of a laser driver as v w /c = ω l /(ω 2 p + ω 2 l ) 1/2 if etching effects 52 are neglected, where ω l = 2πc/λ l , assuming a wavelength λ l = 800 nm. The agreement between numerical and analytical solutions is excellent, showing that this single-particle model accurately describes the coupling between longitudinal and transverse dynamics. From these test cases it can also be seen that betatron oscillations with micron-level amplitude can lead to femtosecond-level slippage and energy differences on the percent range. When considering a particle bunch, these differences in slippage and energy gain between single particles will have an impact on parameters such as the bunch length and energy spread. One immediate consequence will be an increase of the bunch length, which will in turn also lead to increased correlated energy spread if E = 0. Another implication, more critical for FEL applications, will be an intra-bunch slice mixing which will lead to a growth of the slice energy spread.

Impact on slice energy spread
In order to obtain an expression for the induced slice energy spread, Eqs. (6) and (8) can be used to calculate the energy difference ∆γ between two beam particles, p 1 and p 2 , with the same γ 0 but different oscillation amplitude, A 0,p 1 = 0 and A 0,p 2 = 0, which at some time t are in the same infinitesimally long beam slice at ξ s due to the experienced slippage. If E and K are constant along the bunch, then E 0,p i = E 0,s − E ∆ξ p i for i = 1, 2, where E 0,s is measured at ξ s . This allows the particle energy evolution, γ p i (t), to be written in terms of the fields at the current slice. Additionally, if E ∆ξ p i (t)/E 0,s 1, then E 0,p i E 0,s can be assumed for the higher order terms in γ p i (t) while maintaining the full expression in the leading linear term. This also means that ∆ξ p 2 (t) − ∆ξ p 1 (t) A 2 0,p 2 K(Γ −1/2 s (t) − 1)/2cE 0,s , which, added to the previous considerations, leads to where Γ s (t) = 1 + E 0,s t/γ 0 . The correlation ∆γ ∝ −A 2 0 (if E < 0) within a beam slice arises because, due to slippage, particles with a higher A 0 originally come from positions ahead of their current slice. Since E < 0, the accelerating fields they have experienced as they slip towards the back are smaller than at their current position, thus leading to a lower net energy gain. Eq. (9) therefore shows that although a positive correlation ∆γ ∝ A 2 0 exists for particles starting at the same slice, as given by Eq. (8), this correlation is actually negative when looking at particles which have ended in the same slice due to slippage. The induced relative slice energy spread due to this effect can be found from the standard deviation of ∆γ/γ s if the mean slice energy is assumed to evolve asγ s (t) γ 0,s + E 0,s t, from what it is found that whereΓ s =γ s /γ 0,s , σ A 2 is the standard deviation of A 2 0 for the particles within the slice and E 0 is taken at the slice position. For an on-axis Gaussian bunch with zero average transverse momentum (x =ȳ = 0 andp x =p y = 0), it can be obtained that iγ 0 /ε n,i and σ i are, respectively, the beam's normalized emittance 53 , beta function and rms size in both transverse planes, while β m = c/ω 0 , with ω 0 = ω x,0 = ω y,0 , is the matched beta function 54 defined by the focusing fields within the plasma. For a cylindrically symmetric bunch, i.e, if ε n,i = ε n and β i = β , a "mismatch factor" F m = ((M 4 + 1)/2M 2 ) 1/2 can be defined, where M = β /β m , such that A = F m ε n , which leads to the last expression in Eq. (10). For a matched beam (M = 1) the mismatch factor is F m = 1, which also corresponds to its minimum value. The impact of betatron slippage on the slice energy spread will therefore be stronger on mismatched beams.
Analyzing Eq. (10) it can be seen that the induced slice energy spread has a fast initial growth until it reaches a maximum when the bunch energy has increased by a factor 9, i.e, whenΓ s = 9. For higher energies, the slice energy spread slowly decreases proportional toΓ −1/2 s . This behaviour arises from σ ∆ξ γ s , whose growth rate is initially faster thanγ s but gradually slows down due to the reduced slippage at higher energies.
This source of energy spread will act in addition to that due to betatron radiation (BR) 55,56 , which has also been shown to induce a correlation ∆γ ∝ −A 2 0 and generate a slice energy spread given by 57 0,s /E 0 and r e = e 2 /4πε 0 mc 2 is the classical electron radius. As both effects are induced by the same underlying cause (the betatron motion of particles), it is of interest to compare and determine their relative relevance. This can be done by obtaining the intersection points of Eq. (10) with the BR expression, for which perturbation theory is needed. This method yields that the slippage-induced energy spread will initially dominate over BR if the ratio R = C BR /C ∆ξ 0.1, where C ∆ξ = E Kσ A 2 /2cE 2 0 is the coefficient in Eq. (10). When this is the case, a transition into a BR-dominated regime occurs at an energyΓ s =Γ t which, to first order, can be obtained asΓ t R −1/2 . An illustration of both regimes can be seen in Fig. 2 for a particular set of parameters. As a general rule, although this will vary from case to case, the energy spread due to slippage will typically dominate over BR up to energies on the order of ∼ 10 GeV. This is the energy range in which FELs operate, implying that betatron slippage is the source of slice energy spread that should be minimized for these applications.
An effective way for mitigating the slippage-induced energy spread would be to beam-load the wake such that E → 0, thus suppressing the energy spread growth as given by Eq. 10. However, since this condition is difficult to achieve, especially along the whole bunch, other strategies can be derived from this equation. Regarding the accelerating structures, possible strategies include reducing plasma density, since the ratio p , or using hollow plasma channels 58 , so that K → 0. In regards to the accelerated beam, minimizing the emittance, increasing the initial energy and fulfilling matching conditions (F m = 1) would be beneficial, with these two last points being especially relevant for externally injected beams.
Minimizing this slice energy spread is especially important for PBAs based on chirp compensation schemes 32,33,35 . This is because although the average E experienced by the beam after acceleration might be 0, its value at any point of the accelerator is not, thus giving rise to the development of a slice energy spread which will remain in the de-chirped beam. It is possible, however, that the alternation between positive and negative E in these schemes might partially mitigate the correlation ∆γ ∝ A 2 0 from Eq. (9) and thus the slice energy spread, although never remove it completely. As an additional remark, it should also be noted that this correlation between oscillation amplitude and energy caused by betatron slippage has been proposed as a way of performing beam conditioning for FELs 59 by means of a magnetic chicane 60 . Thus, it might also be beneficial in some cases.

Impact on bunch length and total energy spread of ultra-short bunches
In addition to its impact on the slice energy spread, another immediate consequence of the betatron-induced slippage will be a growth of the bunch length. This is especially relevant for ultra-short bunches where the slippage experienced by the particles is comparable to or bigger than the initial bunch length σ z,0 . The contribution of the single-particle slippage to σ z can be directly obtained from the standard deviation of Eq. (6) by assuming an initially monoenergetic bunch with σ z,0 = 0. The bunch lengthening is then found to be whereΓ =γ/γ 0 1 + E 0 t/γ 0 . The last expression on the right is again obtained for the case of a cylindrically-symmetric Gaussian bunch. The bunch length evolution, when taking into account the impact of slippage, is thus given by σ z (t) = (σ 2 z,0 + σ 2 ∆ξ (t)) 1/2 for initially Gaussian bunches. When this is the case, this expression therefore establishes a limitation on the minimum bunch length achievable in a PBA. This could have a strong impact on the production of high-energy, sub-femtosecond bunches, because although certain techniques to inject such ultra-short bunches into the plasma wake have been proposed [15][16][17]31 , their sub-femtosecond duration could be lost when further accelerating them. Possible ways of mitigating this lengthening can be found by looking at Eq. (11), which yields the same conditions on the beam emittance, matching and initial energy found for Eq. (10). In this case, however, no benefits from beam-loading or reduced plasma density are expected because E does not play a role in the lengthening and the ratio K 1/2 /E 0 does not depend on the plasma density. Another remark that can be extracted from Eq. (11) is that, as expected from Eq. (7), the lengthening has a finite upper limit whenΓ → ∞ given by As a consequence of the bunch lengthening, in cases where E = 0, the correlated energy spread of the bunch will also increase. This explains previously observed limitations, where convergence to a finite energy spread was found as σ z,0 → 0 45 . This limit on the minimum energy spread as the initial bunch length is reduced can be calculated for an initially monoenergetic beam with σ z,0 = 0 from the standard deviation of Eq. (8). This leads to which, in the same way as Eq. (11), also tends to a finite value asΓ → ∞. The last expression on the right corresponds again to the case of a cylindrically-symmetric Gaussian bunch. An illustrative example of the phenomena described in this section can be seen in Fig. 3, where an externally injected witness beam is shown at the entry and after 1 cm of plasma, as obtained from an OSIRIS 2D simulation. The induced lengthening as well as the development of a quadratic correlation between energy and oscillation amplitude can be clearly seen. As predicted by Eqs. (8) and (9), this correlation is positive for the overall bunch (Fig. 3b) but negative on the slice level (Fig.  3c). The energy difference ∆γ total was calculated with respect to the average energy of the particles with A 0 0 µm. (c) Distribution of the particle energy with respect to the initial oscillation amplitude within each beam slice. The coloured lines correspond to the distribution within each 20-nm-long slice, while the thin dark lines represent the average. As expected from Eq. (9), the energy-amplitude correlation at the slice level is negative. The analytical prediction was again computed for the central beam slice using the same field values as before. The energy difference within the slice, ∆γ slice , was calculated with respect to the average energy of the particles with A 0 0 µm after removing the longitudinal energy correlation (chirp) along the beam.

Validation against 3D PIC simulations
In order to validate and illustrate the significance of the expressions for the slice energy spread and bunch lengthening, a series of 3D simulations with the PIC codes OSIRIS and HiPACE have been performed. The studied cases consisted on a beam-driven plasma stage providing a 1 GeV net energy gain to an externally injected witness beam, as it would be interesting for FEL applications or multi-staged acceleration. In the OSIRIS simulations, the impact of the slippage effects has been tested for different witness beam parameters and plasma densities, and in all cases beam loading effects could always be neglected. On the contrary, the HiPACE simulations where performed to study the possible mitigating effects of beam loading for a fixed set of plasma and beam parameters, where only the total beam charge was varied. Although only beam-driven cases have been simulated, similar results may apply for a laser driver as long as dephasing effects in the mean beam energy can be neglected, as they are not included in Eqs. (10) to (12).
In all simulated cases, a Gaussian driver with an energy of 1 GeV and a 0.1% spread has been considered. Its dimensions and charge have been defined in terms of the plasma density so that the generated blowout can be scaled with n p . For the baseline case in which n p = 7 × 10 17 cm −3 , as used in previous experiments 6 , this means a transverse size σ x = σ y = 0.4 c/ω p 2.5 µm, a length σ z = c/ω p 6.4 µm, a normalized emittance ε n = 0.4 σ x 1 µm rad and a peak electron density n b 3.78 n p for a total charge Q 265 pC and a peak current I p 5 kA. The plasma target has a flat-top longitudinal profile and is transversely uniform.
For the studies made with OSIRIS, a Gaussian witness beam is injected on-axis at ξ w ξ d − 4.5 c/ω p , where ξ d is the driver centre. Different emittance values (ε n = 0.1, 1 and 10 µm rad) as well as injection energies between 10 MeV and 1 GeV have been tested. The transverse size is matched to the focusing fields so that F m = 1. With regards to the bunch charge and duration, two different sets of parameters have been considered: σ z /c = 3 fs with 1 pC for the studies on the slice energy spread, as it could be obtained in the SINBAD facility at DESY 61, 62 , while a shorter bunch with σ z /c = 0.1 fs with 0.1 pC was used for the bunch lengthening studies. This shorter bunch duration was chosen to better highlight the impact of lengthening due to slippage on sub-femtosecond bunches. The peak current in both cases is too low to significantly modify E along the bunch, therefore allowing beam loading effects to be neglected. The slice energy spread in all tested beams is dominated by slippage, as R 0.1 andΓ R −1/2 . Also, in order to isolate the plasma stage as the only source of energy spread, a zero longitudinal momentum spread at injection has been considered.
In order to compare the simulation results with the analytical theory, the field values in Eqs. (10) and (11) need to be provided. While these could be estimated from analytical models 49,63 , they were measured directly at ξ w from a single simulation with no witness beam, obtaining E 0 0.53 ω p , E −0.37 ω 2 p /c and K 0.5 ω 2 p . This allows for a more accurate 6/13 validation of the model. A view of the plasma wakefields in one of the OSIRIS simulations is shown in Fig. 4. As assumed by the analytical model, a constant K within the blowout and a transversely homogeneous E with uniform E along the witness bunch can be seen. The final simulation results, displayed in Fig. 5, show a strong dependence of the generated slice energy spread on the plasma density and initial beam parameters, reaching values in excess of 10 −2 (or 1%) for certain configurations, way above FEL requirements. Similarly, the increase in bunch length can be on the femtosecond level for typical emittance and energy values, showcasing the difficulty of achieving GeV-class, sub-femtosecond bunches. No variation of the bunch length with the plasma density is observed, as expected from Eq. (11). The analytical model shows excellent agreement with the simulations over several orders of magnitude, with expected discrepancies in some of the higher emittance cases. The differences in Fig.  5(b) arise from the large slippage experienced by particles with higher oscillation amplitude, which causes the transverse distribution of many beam slices to be a truncated Gaussian and thus leads to a smaller energy spread than predicted. The discrepancies in Figs. 5(c) and (d) arise from the longitudinal energy correlation, neglected by Eq. (11), which induces a slight longitudinal bunch compression due to velocity differences between head and tail. In addition to the final slice energy spread and bunch length values, the analytical model also accurately predicts the evolution of these parameters during propagation within the PBA, as shown in Fig. 6. The discrepancies in the energy spread evolution in Fig. 6(a) arise from the initial dephasing between the bunch and plasma wakefield due to the low initial energy (10 MeV), which leads to an experienced accelerating field E > E 0 and, thus, a lower energy spread than predicted by Eq. (10). Furthermore, the initial oscillations on the slice energy spread at the beginning of the simulation correspond to the betatron oscillations of the witness beam, which lead to an oscillation in v z , as seen previously in Fig. 1(b). When the v z of the beam particles is minimum (maximum), more (less) slippage occurs and the energy spread growth is faster (slower). This also leads to small oscillations on the bunch length, as seen in Fig. 6(b). The impact of the v z oscillations on the particle slippage had to be neglected in the theoretical model in order to obtain the analytical expression in Eq. (6). Thus, this effect is not reproduced by the analytical curves shown in Fig. 6.
In order to explore the potential of beam loading for mitigating the generated slice energy spread, an additional set of simulations was performed with HiPACE. The same driver parameters as in the OSIRIS simulations and a plasma density of 7 × 10 17 cm −3 were used. The witness beam parameters were fixed to an initial energy of 100 MeV, ε n = 1 µm rad and σ z /c = 3 fs, while its total charge was varied from 1 to 100 pC. The final energy spread of the central bunch slice after a 1 GeV energy gain is shown in Fig. 7 for the different bunch charges. Only the central slice energy spread is displayed in this case because the effect of beam loading on E varies along the bunch. The analytical points were calculated from the average value during acceleration of the fields at the center of the bunch.
The results show again good agreement between theory and simulations, and that beam loading can significantly reduce the development of slice energy spread (up to a factor of ∼ 20 in this case). For the particular set of parameters under consideration, the energy spread minimum is achieved for a total bunch charge of 40 pC (corresponding to a peak current of ∼ 5.3 kA) and grows again for higher charges because beam loading overcompensates E , which becomes positive. Peak currents in this range of parameters have been experimentally demonstrated 9, 10 and, if optimized, could therefore provide an effective way of mitigating the slice energy spread growth. Discrepancies between the analytical curve and the simulations appear, however, around the energy spread minimum. This is because, in theory, the slice energy spread due to slippage could be completely suppressed for an ideally beam-loaded case in which E = 0 along the whole bunch. This, however, requires an optimized trapezoidal current profile 30 , while the bunch in the simulated case is Gaussian and therefore does not experience a uniform E . Thus, even if E = 0 at the center, a certain energy spread will develop in the central slice due to the slippage of particles originally coming from slices ahead, where E = 0. In addition, other effects apart from particle slippage can lead to increased slice energy spread in beam-loaded cases, such as transverse variations of the accelerating fields in the linear and weakly non-linear regimes 66 .

Discussion
The presented analytical model and simulation results show that the impact of the betatron-induced slippage could be very significant on beam parameters such as the slice energy spread (up to ∼ 1%) and bunch length (growth on the femtosecond range). These effects should therefore be minimized, and guidelines for their mitigation can be derived from the analytical expressions together with the simulated cases. These show that an effective method of suppressing the slippage-induced slice energy spread growth is the use of beam loading (so that E → 0), while another possible solution would be performing the acceleration in hollow plasma channels (K → 0). Otherwise, in order not to rely on these two options, the simulations show that for achieving a slice energy spread < 10 −3 , as typically required by FELs, a plasma density 10 17 cm −3 , a normalized emittance 1 µm rad as well as an initial energy 100 MeV become necessary. The same considerations apply regarding the production of sub-femtosecond bunches, although beam loading or lower plasma densities do not provide any benefit in this case. For applications requiring multiple accelerating stages, such as plasma-based particle colliders, these results imply that the slice energy spread due to slippage would mostly be generated in the first accelerating stage, decreasing at higher energies and eventually being dominated by the emission of betatron radiation.

3D Particle-in-Cell simulations with OSIRIS
The 3D PIC in simulations corresponding to the results shown in Figs. 4, 5 and 6 where carried out with the OSIRIS code 43 using a standard finite differences field solver 67 . Due to the long propagation distance within the plasma, a co-moving simulation box was considered. The simulation grid on the cases shown in Figs. 5 and 6 consisted on 146 × 106 × 106 cells and had dimensions of 8.66 × 6.30 × 6.30 in units of ω p /c. In the graphs shown in Fig. 4 used for visualizing the fields, the simulation grid was slightly longer, with 170 × 106 × 106 cells and dimensions of 10 × 6.30 × 6.30 in units of ω p /c. This was done in order to show the tail of the wakefields, which was otherwise cut off. The time step in all cases was chosen to fulfil the stability criterion 67 . The number of particles per cell was 16 for the beam driver, 1 for the plasma and between 64 and 4096 for the witness beam, depending on the case, in order to have a meaningful total number of particles. The initialization of the electron beams was done using the n_accelerate method in OSIRIS, which progressively increases the longitudinal particle momentum over the specified number of time steps in order to obtain almost self-consistent electromagnetic fields before the beam enters the plasma.

3D Particle-in-Cell simulations with HiPACE
The 3D PIC in simulations corresponding to the results shown in Fig. 7 where carried out with the HiPACE code 44 , which is based on the quasi-static approximation. The simulation grid consisted on 320 × 106 × 106 cells and had dimensions of 8.66 × 6.30 × 6.30 in units of ω p /c. A constant time step size of 5 was used. The number of particles per cell was 16 for the beam driver, 3 for the plasma and 512 for the witness beam in order to have a meaningful total number of particles. The initialization of the particle species in HiPACE is performed by computing the self-consistent fields in each time step by means 9/13 of fast Poisson solvers.

2D Particle-in-Cell simulation with OSIRIS
The data displayed in Fig. 3 showcasing the bunch lengthening and energy spread development comes from a single 2D simulation with OSIRIS of a laser-driven PBA. This simulation considered a plasma density n p = 10 17 cm −3 , a Gaussian laser pulse with a peak power of 1.5 PW and a total energy of 25 J, a peak normalized vector potential a 0 = 4, a wavelength λ = 800 nm, a spot size w 0 = 54 µm and a FWHM duration τ = 15 fs. The witness bunch has a Gaussian profile based on achievable parameters by the ARES linac at SINBAD 61, 62 , featuring a total charge of 0.7 pC, an initial energy of 100 MeV, a duration σ z /c = 0.78 fs, a transverse size with σ x = 5.0 µm and σ y = 5.3 µm, an emittance ε x = 0.70 µm rad and ε y = 0.74 µm rad and an energy spread σ γ /γ = 0.37%. This simulation was also performed with a standard finite differences field solver 67 considering a co-moving simulation box with 5000 × 500 cells and dimensions of 5.95 × 29.76 in units of ω p /c. The time step was chosen to fulfil the stability criterion 67 . The number of particles per cell was 8 for the plasma and 64 for the witness beam. The initialization method of the witness beam is the same as in the 3D case described above.
Numerical solution of the single-particle equations of motion As shown in Fig. 1, the equations of motion for a single particle (Eqs. (1) and (2)) were solved numerically to compare against the analytical expressions. These equations were solved with a Runge-Kutta method of order 4 with a time step of 0.01 fs. The fields experienced by the particles take into account dephasing effects assuming a laser driver with λ = 800 nm.

Calculation of the slice energy spread
Due to the limited amount of particles in a PIC simulation, calculating σ s γ /γ for an infinitesimally short slice as in Eq. (10) is not possible. However, this quantity can be estimated for a finite slice in which the longitudinal energy correlation is removed. The uncorrelated energy distribution can then be obtained as γ unc = γ − δ ξ −ξ , where δ is the slope of the correlation. Additionally, since δ might change along the bunch, the bunches are divided into n slices in which σ s γ /γ is calculated. From there, for the cases shown in Figs. 5 and 6, the slice energy spread of the whole bunch is obtained as the weighted average of the value on each slice. This process is repeated by increasing the number of slices until convergence to a value is reached while still having a significant amount of particles per slice. For the cases to study the effect of beam loading shown in Fig. 7, the value of the slice energy spread only given for the central slice and not as an average.