We rely on high-fidelity simulations using OSIRIS to illustrate this high-laser-coupling-efficiency, low-energy-spread concept. In Fig. 1A, a 125-pC witness electron beam is accelerated through a tri-plateau density profile to ~600 MeV in a wake driven by a laser pulse with only 340-mJ energy. The resulting energy transfer efficiency from the laser to the electron beam is >21%. High efficiency is achieved for relatively low laser normalized vector potential a_0\equiv \frac{e{A}_0}{m_e{c}^2} in the blowout regime [38,39]. According to the scaling law [40], the beam energy gain is E_b\propto a_0^{5/2} and the laser pulse energy is E_l\propto a_0^{7/2} (assuming a round laser pulse cτ ∼ w₀), and the energy transfer efficiency is thus \mathit{\Gamma }\sim \frac{E_b}{E_l}\sim a_0^{-1}. Here, we chose a₀ = 1.67 to ensure the excitation of plasma wake in the blowout regime and a relatively high energy transfer efficiency. The resulting laser power is not sufficient for self-guiding [41], so a tri-plateau plasma structure with a parabolic transverse density of the form n_p\left(r\right)=n_{p0}\left(1+\frac{\mathrm{\Delta }n{r}^2}{w_0^2}\right) is used, where Δn = 0.5 such that the laser with a spot size w₀ is free of diffraction throughout the 6-mm-long plasma. Because the density of each plateau increases in discrete steps (Fig. 1B), the size of the wake cavity shrinks as it transits from one stage to the next stage (Fig. 1A). The on-axis densities, n_p0, of the 3 plateaus are 2 × 10¹⁸ cm⁻³, 3.13 × 10¹⁸ cm⁻³, and 8.3 × 10¹⁸ cm⁻³, and the corresponding lengths of plateau are 1.73, 2.67, and 0.15 mm, respectively. Adjacent plateaus are connected by a linear up-ramp of 0.5 mm, and the slopes at both the entrance and the exit for the first and last plateau are 0.25 mm long. Most of the energy gain occurs in the uniform density plateau regions—in this sense, the exact lengths of the density up-ramps between the stages are not critical. The laser pulse has a Gaussian transverse profile, focused to a beam waist w₀ = 10.7 μm at the midpoint of the density up-ramp of the first plateau, and a sin²-shaped temporal envelope of the field with a 21.5-fs full width at half maximum (FWHM) pulse length. The laser pulse is nonlinearly chirped to partially compensate the nonlinear dispersion (see Materials and Methods). For ease in interpretation, a bi-flattop electron beam with 3.76-μm beam length is initialized behind the drive laser with an initial energy of 51.1 MeV (γ = 100) and 0.1% energy spread, corresponding to the situation of external injection or staging. The resulting beam energy is insensitive to the specific longitudinal profile for a Gaussian-like laser pulse, and >21% energy efficiency with <1% energy spread can be achieved. The initial emittance is 0.1 mm·mrad, and the beam is focused to the midpoint of the first up-ramp with a spot size of 0.5 μm. Compared with the uniform (the yellow line in Fig. 1B) and bi-plateau (the red lines in Fig. 1B) channel structure, the final energy gain of the tri-plateau structure is almost tripled, significantly increasing the acceleration efficiency. The length of the uniform plasma is 8 mm, and the lengths of the 2 plateaus for the bi-plateau are 1.73 and 2.92 mm, respectively. Numerical parameters are given in Materials and Methods.

Because ultra-relativistic electrons essentially move at the speed of light, the accelerated e-beam will gradually outrun the accelerating phase of the plasma wake excited by lasers (dephasing). This causes the acceleration gradient felt by the electron beam to decrease or even become negative as the laser propagates in a uniform plasma. Since the dephasing length is typically much shorter than the pump depletion length, the acceleration gradient will be significantly reduced far before the laser energy is depleted, thus greatly limiting the acceleration efficiency. Due to this dephasing process, the spatial gradient of the energy gain in the laboratory frame gradually reduces and appears parabolic as depicted by the yellow dashed lines in Fig. 1C. On the other hand, during the acceleration process, the laser loses energy in a nearly linear manner as shown by the yellow dashed lines in Fig. 1E, and eventually, only one-third of the initial energy is transferred to the plasma wake. In our tri-plateau scheme, a rephasing process occurs in both the density up-ramp and the higher density plateau that follows it. In the rephasing process, the wavelength of the plasma wake shrinks due to the increase of the plasma density and e-beam finds itself once again close to the maximum acceleration phase, significantly improving the acceleration efficiency. Since the acceleration gradient also increases with approximately square root of the density, the slope of the gain curve (blue line) increases sharply from one plateau to the next. As shown in Fig. 1C, the energy gain of each stage is insensitive to the plasma density, albeit the acceleration distance (being proportional to the dephasing length) becomes shorter as the density increases. Using the scaling laws in the blowout regime [40], for the accelerating field, E_z\propto n_p^{1/2}{a}_0^{1/2}, and the dephasing length, L_d\propto n_p^{-1}, the energy gain scales as \mathrm{\Delta }E\propto n_p^{-1/2}{a}_0^{1/2}. Due to the focusing and steepening of the laser pulse, a₀ increases so that ΔE does not change much as n_p increases. As a result, the energy gain is nearly tripled compared with that of the uniform plasma case.

As shown by the blue solid line in Fig. 1C, the energy gain curve now increases almost linearly in each stage, and hence, a high beam loading efficiency is maintained continuously for the 125-pC beam charge. The beam loading efficiency (wake-to-beam efficiency), defined as the ratio of energy extracted by the particles from the wake to the energy cost of the laser, is calculated by N〈∆W〉/(W_i − W_o) in each plateau, where N is the electron number of the beam, 〈∆W〉 is the averaged energy gain of the electron beam in this plateau, and W_i/W_o is the laser energy entering/exiting this plateau. In the first 2 plateaus, the wake-to-beam efficiency is 49.6% and then drops to 22.3% in the third plateau. More laser energy is utilized in the tri-plateau plasma; 67.5% of the laser energy is consumed in the proposed scheme, while only 34.0% is consumed in a single uniform plasma channel. The total efficiency is the product of how much of the laser energy is transferred to the wake (roughly the laser energy utilization percentage) times the beam loading efficiency, which is >30% in this case. As a result, a total efficiency of >21% is obtained, which is the highest value achieved to date for a fully self-consistent LWFA concept.

We note that, in this proposed scheme, the laser energy is efficiently utilized because the pump laser pulse has a nonlinear chirp with frequency increasing from the front to the middle and then decreasing from the middle to the back, allowing the pump waveform to self-compensate for photon deceleration due to wake formation and maintain its shape. The nonlinear frequency chirp is of the form \frac{k\left(\xi \right)}{k_0}=c_0+c_1\xi +c_2{\xi }^2, where k(ξ) and k₀ are the local and central wave number and ξ ≡ ct − z. The chirp coefficients c₀, c₁, and c₂ are 1.2, 1.79 × 10⁻³k₀, and -9.7\times {10}^{-5}{k}_0^2, respectively. Such a nonlinear chirp can be realized by either nonlinear cross-phase modulation [42] or nonlinear pulse compression by customized gratings [43,44]. In the blowout regime [38,39], the driving laser pulse with the optimal pulse length resides at the first half of the ion cavity where the photon deceleration rate r_ph < 0 (see Materials and Methods). The occurrence and effects on pulse distortion of photon deceleration/pump depletion are illustrated in Fig. 2. The left column frames are taken from a simulation with a transform-limited pulse, while the right column is for a laser with a nonlinear frequency chirp. For these examples, a laser with identical parameters as above is sent through a single plateau with n_p0 = 3 × 10¹⁸cm⁻³ and Δn = 0.5. For each case, the axial lineout and spectrogram of the laser are shown at propagation distances, 0, 3.81, and 5.08 mm, respectively. Initially, the middle of the pulse experiences a more intense photon deceleration than the head and tail of the pulse (see the inset in Fig. 2A). This can be seen from Fig. 2A, C, and E, where the on-axis spectrogram for the transform-limited laser evolves into a V shape, indicating as expected that the middle part suffers from a faster photon deceleration. The pulse is significantly lengthened as a result of the resulting envelope distortion (Fig. 2E). On the other hand, as shown in Fig 2B, D, and F, the initial nonlinear chirp nearly perfectly balances the redshifting, and the waveform is well preserved throughout the simulation. Off-axis evolution of the spectrogram has a similar pattern shown by Fig. 2, which is supplied in Materials and Methods.

Such a pulse can propagate much further than the pump depletion limit of a transform-limited pulse of similar pulse width and peak intensity. The importance of this chirp is seen in the gray line of Fig 1C, where it can be seen that no energy gain arises in the third stage without the chirp (transform-limited pulse). The reason for saturation of the energy gain can be seen in Fig. 1D to be that the pulse length increases rapidly for the laser pulse that has no chirp in the final stage. In addition to the high efficiency, the beam energy spread of <1% can be obtained owing to the DBL. The mechanism of how the DBL effect maintains a low energy spread will be elaborated on in the later section.

The use of stepwise plasma density plateaus connected with density up-ramps not only significantly increases the beam energy gain but also is capable of accelerating a large amount of beam charge with an extremely small energy spread through a DBL effect. The beam loading effect can be understood by tracking the accelerating field felt by the beam [13]. In Fig. 3A to D, the laser field, lineout of the accelerating field (green line), plasma, and beam density are plotted at the start to end points of each transition section (point a to point d) in Fig. 1B. In Fig. 3E to I, the beam energy versus axial position is shown at the same propagation distances. It can be seen that the plasma wake is overloaded () during the first and at the entrance of the second plateau (Fig. 3A and B), leading to a negative energy chirp in the electron beam. On the other hand, in the third plateau, and during the exit of the second and third plateaus, the wake is underloaded (\frac{d{E}_z}{\mathit{d\xi }}>0) (Fig. 3C and D) such that the energy chirp caused by the first plateau is almost perfectly removed. The transition from an over- to underloaded wake happens in the second plateau (between Fig. 3B and C).

This mechanism can be qualitatively described by the nonlinear theory of the blowout regime [38]. We find that the gradient of E_z (see Materials and Methods) can be expressed as

where \lambda =k_p^2{\int }_0^{+\infty }\frac{n_b}{n_p}\mathit{rdr} is the normalized linear charge density of the beam, E_p\equiv \frac{m_{e}c{\omega }_p}{e} is the characteristic electric field in the plasma, and r_b is the blowout radius. The linear charge density λ, which is independent of plasma density, can be assumed to not depend on the propagation distance. However, the value of r_b at the location of the beam changes (increases) from dephasing within an individual plateau, and k_p increases when crossing into the next plateau. For the parameters simulated, during the first plateau, the dominant term is -\frac{2\lambda }{k_p{r}_b^2}, while it becomes \frac{k_p}{2} in the third plateau. In the second plateau, the dominant term changes from -\frac{2\lambda }{k_p{r}_b^2} to \frac{k_p}{2} as r_b increases because of the increased plasma density and dephasing. The final energy chirp is obtained by integrating Eq. 3 with respect to the propagation distance, and eventually vanishes at the optimal distance (Fig. 3I). As seen in Fig. 3E to I, the beam energy chirp increases first and then decreases, and the residual energy chirp is almost fully compensated, which agrees with the analysis above. The final relative energy spread reaches as low as 0.83% (FWHM), while a considerable beam charge of 125 pC is accelerated.

Assuming a fixed a₀ and λ₀, and that the laser is guided by a channel, the key physics can be viewed approximately self-similar, and scaling laws can be straightforwardly obtained (see Materials and Methods). When the plasma density changes from n_p to n^′_p, the characteristic lengths of the key physics can be scaled by a factor F=\sqrt{n_p/n_p^{\prime }}. If we proportionally scale the focal waist and pulse length of the laser through w_0^{\prime }=F·w_0 and τ′  = F · τ, the acceleration gradient, maximum beam energy gain, and the number of electrons that can be accelerated scale as E_z′  = F⁻¹E_z, 〈∆W〉′  = F² · 〈∆W〉, and N′ = F · N, respectively, according to the scaling law of LWFAs [40] (see Materials and Methods for the derivation). The energy transfer efficiency η′ is an invariant under this scaling since.

We carried out simulations with F = 2 and F = 4 to validate the scaling law, and the results are summarized in Table 1. In these simulations, the initial parameters of laser, plasma, and electron beam are scaled from case 1 and the results (energy gain, efficiency) show excellent agreement with what the scaling law predicts as shown in Fig. 4B. In each case, the final energy spread remains below 1%. Case 3 corresponds to a 500-pC electron beam being accelerated to ~10 GeV by a 22-J laser pulse, which is within the capability of modern petawatt lasers. The calculated energy transfer efficiency greatly exceeds the results in previous LWFA experiments and simulations [45,46]. Experimental verification of case 3 can be attempted in the near term. Furthermore, we can extrapolate the proposed scheme to F = 16 in column 4. In this case, an LWFA driven by a 1.4-kJ, 340-fs laser pulse would accelerate a 2-nC electron beam with ~150-GeV energy with an average acceleration gradient of 6 GeV/m. Laser with such energy and beams with such charge are already available. The proposed scheme could thus provide the electron arm of a Higgs factory, with ~20% efficiency and sub-1% energy spread. For comparison, a scaled case (case 5) using a uniform plasma channel is also listed in Table 2. In this case, only about one-third of the energy gain in case 4 is achieved, and the energy spread grows up to 15%. To scale the scheme to even higher energy approaching 10 TeV, in-depth studies on the limitations caused by physics such as ion motion should be conducted in the future.

The propagation and the evolution of the laser pulse are of crucial importance in LWFA research. Thus, high-fidelity modeling of the laser–plasma interaction including nonlinear dispersion effects is critical for the calculation of energy gain and energy transfer efficiency. Standard PIC simulations utilize a second-order Yee solver to solve the electro-magnetic fields, which could suffer from significant numerical dispersion, numerical Cerenkov radiation (NCR), and self-fields from the beam for the propagation distances considered here. These numerical errors may cause the PIC simulations to overestimate the dephasing effect that is closely related to the energy efficiency from laser pulse to the electron beam, and the growth of the emittance and energy spread of the beams. In this article, a customized finite difference EM solver, which can eliminate the numerical dispersion, NCR, and self-forces, is adopted [47–50].

All the simulations are conducted using the quasi-three-dimensional (3D) version of the PIC code OSIRIS. For the F = 1 case, the simulation window has a dimension of 150k_0^{-1}\times 400k_0^{-1} with 500 × 2,000 cells in the r and z directions, respectively. This corresponds to cell sizes of \mathrm{\Delta }r=0.30k_0^{-1} and \mathrm{\Delta }z=0.20k_0^{-1} (where the wave number k₀ = 2π/λ₀ and the wavelength λ₀ = 800 nm). The time step is chosen as \text{∆}t=0.1c^{-1}{k}_0^{-1}. In the scaled cases of F = 2 and F = 4, the sizes of the simulation window and transverse cell size are scaled accordingly, while the longitudinal cell size is kept at \mathrm{\Delta }z=0.20k_0^{-1}. There are 10 macroparticles in each cell and 8 duplications in the azimuthal direction. Only the m = 0 and m = 1 azimuthal modes are kept in the simulations.

Due to the self-phase modulation of the laser pulse in the plasma, the laser pulse loses its energy by the frequency downshifting (photon deceleration) with the rate expressed as [51–53]:

where k is the local wave number of the laser pulse. This photon deceleration rate has a distribution along ξ (the inset of Fig. 2B) and is almost nonevolving during most of its propagation length. In the proposed high-efficiency scheme, the nonlinear chirp is added to the laser pulse to offset this photon deceleration. Another option is to shape the pulse to provide a wake with a constant photon deceleration rate throughout the entire pulse [54].

In addition to the on-axis waveform evolution of the nonlinearly chirped pulse (Fig. 2), off-axis waveform evolutions with different transverse position up to 11.1 μm are provided in Fig. 5 for integrality, noting that the matched w₀ = 10.7 μm. Each row shows the laser time–frequency distribution at a different time when the laser pulse propagates 0, 7.39z_R and 11.36z_R, with z_R = 447 μm is the Rayleigh length, while each column shows the spectrograms with different transverse positions at x₂ = 3.7 μm, 7.4 μm, and 11.1 μm, noting that the beam waist of the initial laser pulse w₀ = 10.7 μm. The off-axis waveforms show similar evolution with the on-axis waveform (Fig. 2B, D, and F).

To derive the expression of E_z, the nonlinear theory of the bubble regime [38] is followed here. In the ultra-relativistic limit where the maximum blowout radius k_pR_b ≫ 1, the trajectory of the innermost particle is given bywhere r_b(ξ) is the radial position of the innermost particle in the bubble sheath and \lambda =k_p^2{\int }_0^{+\infty }{n}_b/n_p\mathit{rdr} is the normalized linear charge density of the electron bunch. The longitudinal electric field E_z can be expressed as a function of r_b[40],where ψ ≡ ϕ − A_z/c is the pseudo-potential of the plasma wake. The gradient of the accelerating field dE_z/dξ is obtained via differential of Eq. 5,

In Eq. 6, the term \frac{k_p}{2}{\left(\frac{d{r}_b}{d\mathit{\xi }}\right)}^2 is omitted since it is negligible for most part of the bubble.

The initial parameters for the scaled simulations are guided by LWFA scaling laws [40], starting from the scaling of the accelerating wakefield structure. For convenience, real and normalized units are used in this section. Through simulations, it has been found that a relatively stable propagation of the ion cavity is realized when k_p{w}_0\simeq k_p{R}_b\simeq 2\sqrt{a_0}, where w₀ is the initial laser beam waist and R_b is the blowout radius of the ion cavity [38,39]. For the parameters being considered here, these relationships hold in a parabolic density channel if the density at the bottom of the channel is used. The connection between the plasma density and the size of the ion cavity is established as n_p ∝ F⁻² because k_p\propto \sqrt{n_p}, where F is the scaling factor of the radius of the ion cavity. To efficiently drive the wake, the laser pulse length should fill the first half of the ion cavity; thus, cτ_FWHM ∝ R_b, where τ_FWHM is the pulse duration. Hence, the laser pulse energy scales with F³. The length of the plasma structure L_acc is characterized by the dephasing length L_{\varphi }=\frac{4}{3}\frac{k_0^2}{k_p^2}\frac{\sqrt{a_0}}{k_p}, providing L_acc ∝ F³. Using the scaling law of the accelerating field, E_z\propto \frac{m_{e}c{\omega }_p}{e}\propto F, the expected beam energy gain is obtained as 〈∆W〉 ∝ F². Another initial parameter to be determined is the loaded electron number N, and based on [13,38,55], it should scale with the electron number that is expelled from the ion cavity, leading to N\propto n_p{R}_b^3\propto F. The relevant formulas and scaling factors are summarized in Table 2.

As for the PIC simulations, the simulation window should scale with the ion cavity, e.g., r_window ∝ F and z_window ∝ F. The transverse cell size ∆r scales with R_b to resolve the plasma wavelength, while the longitudinal cell size ∆z is fixed to resolve the laser wavelength.