J. Beckers, F. Kaspar and J. Knollmüller on behalf of the COMPASS
Collaboration Technische Universität München, Physik Department, James-Franck-
Straße 1, 85748 Garching bei München Excellence Cluster Origins,
Boltzmannstraße 2, 85748 Garching bei München

# Progress in the partial-wave analysis methods at COMPASS

J. Beckersins:tum<EMAIL_ADDRESS>F. Kasparins:tumins:cl J.
Knollmüllerins:tumins:cl for the COMPASS Collaboration
florian.kaspar@tum.dejakob.knollmueller@tum.deins:tumins:tumins:clins:cl

###### Abstract

We study the excitation spectrum of light and strange mesons in diffractive
scattering. We identify different hadron resonances through partial wave
analysis, which inherently relies on analysis models. Besides statistical
uncertainties, the model dependence of the analysis introduces dominant
systematic uncertainties. We discuss several of their sources for the
$\pi^{-}\pi^{-}\pi^{+}$ and $K^{0}_{S}K^{-}$ final states and present methods
to reduce them. We have developed a new approach exploiting a-priori knowledge
of signal continuity over adjacent final-state-mass bins to stably fit a large
pool of partial-waves to our data, allowing a clean identification of very
small signals in our large data sets. For two-body final states of scalar
particles, such as $K^{0}_{S}K^{-}$, mathematical ambiguities in the partial-
wave decomposition lead to the same intensity distribution for different
combinations of amplitude values. We will discuss these ambiguities and
present solutions to resolve or at least reduce the number of possible
solutions. Resolving these issues will allow for a complementary analysis of
the $a_{J}$-like resonance sector in these two final states.

##  1 Introduction

At COMPASS, we study the spectrum of light mesons. They are produced as short-
lived intermediate states in so-called diffractive-dissociation reactions, in
which a high-energy hadron beam interacts strongly with a proton target. The
decay daughters are measured in the spectrometer. In order to disentangle the
contributions from the interfering states and extract their parameters, we
perform a partial-wave analysis (PWA). In the following, we will present the
advances in the PWA methods at COMPASS, in the $K_{S}^{0}K^{-}$ and
$\pi^{-}\pi^{-}\pi^{+}$ final states. In this paper, we focus on the studies
in the former channel, while a detailed report on the methods developed in the
latter is given in [1].

##  2 Partial-wave analysis method

We will first summarise the method of partial-wave analysis, detailed in [2].
In order to gain information about the states produced in the diffractive-
dissociation process, we model its intensity distribution in the invariant
mass $m_{X}$ and in the $n$-body-phase-space variables $\tau_{n}$. We separate
the amplitude describing the process into contributions with specific quantum
numbers, so-called partial waves denoted by the index $a$, so that

$\mathcal{I}(m_{X};\tau_{n})=\left|\sum_{a}\,\mathcal{T}_{a}(m_{X})\,\psi_{a}(m_{X};\tau_{n})\,\right|^{2}\
.$ (1)

In eq. 1, each partial-wave amplitude is further subdivided into two parts.
The decay amplitude $\psi_{a}$ describes the phase-space distribution of the
daughter particles of a state with quantum numbers $a$. The production
amplitude $\mathcal{T}_{a}$ is what we want to extract from the data, as it
contains information about the produced mesonic states. To do this in a quasi-
model-independent manner, the data is usually divided in kinematic bins of
$m_{X}$, and then fitted with our model to extract the values of the
production amplitudes in each bin. These may then serve as input to a second
fit, in which their mass dependence is modelled in terms of resonances and the
latter’s parameters extracted.

##  3 Ambiguities in the partial-wave analysis of the $K_{S}^{0}K^{-}$ final
state

Figure 1: Partial-wave intensities of the ambiguous amplitudes. The model and
calculated solutions are drawn as orange and blue curves respectively. The
intensities obtained in the PWA fit are shown as black dots with grey error
bars. See also [1].

###  3.1 Non-uniqueness of the amplitude values

For any final state with two spinless particles, such as $K_{S}^{0}K^{-}$, the
decomposition of the intensity in eq. 1 is not unique in each $m_{X}$ bin, and
several sets of amplitudes $\\{\mathcal{T}_{a}\\}$ result in the same
intensity distribution in phase space. To show this, and following [3], we
start by simplifying eq. 1 to111We use the reflectivity basis for the
amplitudes [5]. We neglect contributions with spin-projection $M\neq 1$
because of the dominance of Pomeron exchange at COMPASS energies [2] and the
suppression of higher $M$ contributions. We drop indices
$M^{\varepsilon}=1^{+}$.

$\mathcal{I}(\theta,\phi)=\big{|}\sqrt{2}\,\underbrace{\sum_{J=1}^{J_{\mathrm{max}}}\,\mathcal{T}_{J}\,Y_{J}^{1}(\theta,0)}_{a(\theta)}\,\big{|}^{2}\
|{\sin{\phi}}|^{2}\ .$ (2)

By expressing the spherical harmonics $Y_{J}^{1}(\theta,0)$ as functions of
the sole variable $u\equiv\tan{(\theta/2)}$, $a(\theta)$ in eq. 2 can be
written as a polynomial in $u^{2}$, which complex coefficients are linear
combinations of $\\{\mathcal{T}_{J}\\}$. Finding the polynomial’s
$J_{\mathrm{max}}-1$ roots $r_{k}$, also called Barrelet Zeros [4], enables us
to write $a(\theta)$ as

$(1+u^{2})^{J_{\mathrm{max}}}\,a(u)=u\
c\prod_{k=1}^{J_{\mathrm{max}}-1}(u^{2}-r_{k})$ (3)

by root decomposition. The product above enters only as absolute value in the
intensity model from eq. 2. Therefore, the latter is invariant under complex
conjugation of a single or several roots. Since the values of the roots
$r_{k}$ depend non-linearly on $\\{\mathcal{T}_{J}\\}$, this operation
evidently leads to different amplitude values. We can compute them from a
starting set of amplitudes, by obtaining the roots of
$a(\theta\,;\\{\mathcal{T}_{J}\\})$ numerically. The new polynomials and sets
of amplitude values can then be calculated from the $2^{J_{\mathrm{max}}-1}$
combinations of complex-conjugated roots.

###  3.2 Pseudodata study of the ambiguities

We have studied these ambiguities in the PWA of the $K_{S}^{0}K^{-}$ final
state. To this end, we use a model containing four partial-wave amplitudes,
parameterised close to physical amplitudes, i.e. continuous in $m_{X}$ and
containing resonant-like structures. The model amplitude of the $J^{P}=3^{-}$
wave is zero. First, we compute the exact distributions of the ambiguous
amplitudes, by sampling in $m_{X}$ and, at each point, using the values of the
model amplitudes as input for the computation of the ambiguous solutions. Fig.
1 shows the resulting intensity distributions. The model is drawn as an orange
curve, and the ambiguous solutions in blue. Observing fig. 1, we see that the
ambiguous amplitudes are also continuous in $m_{X}$, but that the ambiguities
lead to distributions that are significantly different from the starting
model, especially in the smaller waves $J^{P}=1^{-}$ and $3^{-}$. We also note
that the highest-spin wave (here $J^{P}=4^{+}$) is unaffected by the
ambiguities. In a second step, we generate pseudodata according to the
amplitude model and perform a PWA.222To ensure that we find all solutions, we
perform a large number of fitting attempts with random starting values. The
resulting partial-wave intensities are shown as black dots in fig. 1. Overall,
the amplitudes found by the fit agree with the calculated distributions. Due
to finite data however, the intensity distribution is distorted in some mass
bins and the number of solutions is reduced. In addition, we see that the
amplitude of the $J^{P}=4^{+}$ wave is still invariant.

###  3.3 Resolving the ambiguities

The ambiguities presented above are inherent to the $K_{S}^{0}K^{-}$ final
state, but their influence on the partial-wave analysis can be alleviated in
the COMPASS data by the choice of waveset. Looking at the lower left plot of
fig. 1, we see that intensity appears in the ambiguous solutions of the
$J^{P}=3^{-}$ wave although the model amplitude is zero. Leaving this wave out
of the sum, effectively setting its amplitude to zero, hence does not prevent
to correctly find the model solution. This, however, changes the intensity
distribution of other solutions with a non-zero $3^{-}$ amplitude such that
they are no longer indistinguishable. This can be used in the COMPASS
analysis: in our data, contributions with odd spins $J$ are suppressed.
Excluding them from the waveset should matter only little in the PWA but
effectively resolve the ambiguities in most mass bins.333We could also achieve
this by manually sorting out ”unphysical” solutions, i.e. with significant
intensity in these waves, but this way resolves the ambiguities directly in
the fit procedure. We have confirmed this by performing a PWA fit of the
pseudodata without the $3^{-}$ wave: there, only the model solution is found.

##  4 Continuity for partial-wave analysis

We have seen, for the $K_{S}^{0}K^{-}$ final state, how ambiguities and noise
can impact the quality of the fit results. This affects all such analyses,
since the number of parameters in the fit models are large compared to the
available data and the likelihood often shows multiple optima.

We have developed a new approach that tries to reduce the uncertainties on the
obtained solutions by making use of additional prior information. We know that
the underlying signals are continuous and follow the reaction kinematics,
meaning they should be suppressed both at threshold and for large final-state
masses. In the binned approach, we extract the amplitudes in a quasi-model-
independent way. We want our new model to keep this flexible behaviour. It is
based on Information Field Theory (IFT) [6, 7, 8] and allows us to fit smooth
curves to all bins simultaneously, thus enforcing our prior requirements. We
are also able to extend the model description with resonance
parameterisations, allowing us to directly perform resonance fits in a single
step.

We studied the new method on simulated data of the $\pi^{-}\pi^{-}\pi^{+}$
final state. We generated data according to the continuous model and tried to
recover the input using both the binned approach and the new method. Fig. 2
demonstrates how combining information across multiple bins can drastically
reduce the overall uncertainty of the fit as compared to the usual method.

Figure 2: _Left:_ Intensity of input model (orange) split in resonant (blue)
and non-resonant (green) components. _Right:_ Comparison of the two different
fitting methods. The binned fit (black points) shows large statistical
uncertainties. The new method (red) recovers the input model (orange) with
high accuracy. It also correctly separates the components (blue and green).
More examples can be found in [1].

##  5 Outlook

As discussed, the ambiguities in the partial-wave analysis of the
$K_{S}^{0}K^{-}$ final state can be reduced. We are proceeding with the
partial-wave analysis of the COMPASS $K_{S}^{0}K^{-}$ data and are extracting
parameters of the $a_{J}$ states appearing therein. We are applying the
continuous PWA method on the large COMPASS $\pi^{-}\pi^{-}\pi^{+}$ dataset to
measure the contributing $a_{J}$ and $\pi_{J}$ resonances with increasingly
higher precision and to improve our sensitivity to small signals. In addition,
we are working on applying this novel method on other channels measured at
COMPASS. Initial results on simulated data show that this approach may help to
separate the different ambiguous solutions in $K_{S}^{0}K^{-}$.

###### Acknowledgements.

We would like to thank Philipp Frank and Torsten Enßlin (both Max-Planck
Institute for Astrophysics), as well as Stefan Wallner (Max-Planck Institute
for Physics), who, together with the authors, made a first attempt to adapt
NIFTy for partial-wave analyses. We also thank the COMPASS Hadron Subgroup.
Some of the results in this publication have been derived using the NIFTy
package [8].
Funded by the DFG under Germany’s Excellence Strategy - EXC2094 - 390783311
and BMBF Verbundforschung 05P21WOCC1 COMPASS.

## References

  * [1] F. Kaspar, J. Beckers J. Knollmüller in Proceedings of the MESON 2023 Conference, EPJ Web Conf.  (forthcoming)
  * [2] B. Ketzer, B. Grube D. Ryabchikov Prog. Part. Nucl. Phys.1132020 103755
  * [3] S. U. Chung Phys. Rev. D561997 7299–7316
  * [4] Barrelet Nuov. Cim. A81972 331-371
  * [5] S. U. Chung and T. L. Trueman Phys. Rev. D111975 633–646
  * [6] T. A. Enßlin, M. Frommert F. S. Kitaura Phys. Rev. D802009 105005
  * [7] P. Arras, P. Frank, P. Haim others Nat Astron62022 259–269
  * [8] P. Arras and M. Baltac and T. A. Ensslin and P. Frank and S. Hutschenreuter and J. Knollmüller and R. Leike and M. Newrzella and L. Platz and M. Reinecke and others Astrophysics Source Code Library2019 https://gitlab.mpcdf.mpg.de/ift/NIFTy


# Multi-GeV Wakefield Acceleration in a Plasma-Modulated Plasma Accelerator

J. J. van de Wetering<EMAIL_ADDRESS>John Adams
Institute for Accelerator Science and Department of Physics, University of
Oxford, Denys Wilkinson Building, Keble Road, Oxford OX1 3RH, United Kingdom
S. M. Hooker  John Adams Institute for Accelerator Science and Department of
Physics, University of Oxford, Denys Wilkinson Building, Keble Road, Oxford
OX1 3RH, United Kingdom  R. Walczak  John Adams Institute for Accelerator
Science and Department of Physics, University of Oxford, Denys Wilkinson
Building, Keble Road, Oxford OX1 3RH, United Kingdom  Somerville College,
Woodstock Road, Oxford OX2 6HD, United Kingdom

###### Abstract

We investigate the accelerator stage of a Plasma-Modulated Plasma Accelerator
(P-MoPA) [Phys. Rev. Lett. 127, 184801 (2021)] using both the paraxial wave
equation and particle-in-cell (PIC) simulations. We show that adjusting the
laser and plasma parameters of the modulator stage of a P-MoPA allows the
temporal profile of pulses within the pulse train to be controlled, which in
turn allows the wake amplitude in the accelerator stage to be as much as
$72\%$ larger than that generated by a plasma beat-wave accelerator with the
same total drive laser energy. Our analysis shows that Rosenbluth-Liu detuning
is unimportant in a P-MoPA if the number of pulses in the train is less than
$\sim 30$, and that this detuning is also partially counteracted by increased
red-shifting, and hence increased pulse spacing, towards the back of the
train. An analysis of transverse mode oscillations of the driving pulse train
is found to be in good agreement with 2D PIC simulations. PIC simulations
demonstrating energy gains of
$\sim$1.5\text{\,}\mathrm{G}\mathrm{e}\mathrm{V}$$
($\sim$2.5\text{\,}\mathrm{G}\mathrm{e}\mathrm{V}$$) for drive pulse energies
of $2.4\text{\,}\mathrm{J}$ ($5.0\text{\,}\mathrm{J}$) are presented. Our
results suggest that P-MoPAs driven by few-joule, picosecond pulses, such as
those provided by high-repetition-rate thin-disk lasers, could accelerate
electron bunches to multi-GeV energies at pulse repetition rates in the
kilohertz range.

##  I Introduction

First proposed in 1979 by Tajima and Dawson, laser wakefield accelerators
(LWFA) [1] can outperform the accelerating gradients of conventional radio-
frequency cavities by up to three orders of magnitude. In a LWFA, an
ultrashort laser pulse excites a plasma wave “wake” by pushing free electrons
via the ponderomotive force whilst the heavier ions remain approximately
stationary. The separation of electrons and ions generates large electric
fields, which propagate with a phase velocity set by the group velocity of the
laser pulse, and are therefore suitable for accelerating relativistic charged
particles.

For a laser pulse to efficiently excite a plasma wave, its duration must be
less than half the plasma period $T_{p}=2\pi/\omega_{p}$, which is in the 100
fs range for plasma densities of interest. Because sufficiently short joule-
scale pulses were not yet feasible, Tajima and Dawson proposed the plasma
beatwave accelerator (PBWA) [1]. In this configuration, two long co-
propagating pumps with a carrier frequency mismatch approximately equal to the
plasma frequency $\Delta\omega=\omega_{1}-\omega_{2}\approx\omega_{p}$
interfere to form a train of short pulses via beatwave modulation. The
resulting cosine-squared intensity modulation corresponds to a train of pulses
spaced by $T_{p}$ that can resonantly excite a plasma wave [1, 2, 3, 4]. More
generally, we will refer to resonant wakefield excitation with trains of
uniformly or non-uniformly spaced pulses as multipulse-LWFA (MP-LWFA) [5].

The great majority of recent LWFA experiments have instead used single high
intensity ultrashort $(<$100\text{\,}\mathrm{f}\mathrm{s}$)$ laser pulses from
Ti:sapphire laser systems enabled by the development of chirped pulse
amplification (CPA) [6]. At high pulse energies ($>$1\text{\,}\mathrm{J}$$),
these laser systems are restricted to low ($\sim$ 0.1-10 Hz) repetition rates
[7] and have poor ($<0.1\%$) electrical-to-optical energy efficiencies [8].
These limitations reduce the number of applications for which these LWFAs
offer an advantage over conventional, radio-frequency particle accelerators.

To tackle these issues, a new scheme was recently proposed [9], which we have
dubbed the Plasma-Modulated Plasma Accelerator (P-MoPA). This approach aims to
drive high repetition rate, multi-GeV plasma accelerators with thin-disk
lasers which are both efficient and can provide multi-joule laser pulses at
kHz repetition rates [10, 11, 12]. These laser systems cannot be used to drive
LWFAs directly since the duration of the pulses they provide is too long:
$\tau\gtrsim$1\text{\,}\mathrm{p}\mathrm{s}$$ for joule-scale pulses [13, 14,
15]. To overcome this limitation, in a P-MoPA the drive pulse is spectrally
broadened by adding multiple sidebands spaced by the plasma frequency using a
‘modulator’ stage. This is followed by a dispersive optical system that
removes the spectral phase exhibited by the sidebands, to generate a train of
short pulses that can be used to resonantly excite a wakefield in an
accelerator stage.

In this paper we study the performance of multi-GeV MP-LWFAs driven by P-MoPA
pulse trains; it builds on our earlier work that established the range of
parameters for which operation of the modulator stage was stable [16]. We
derive a full 3D analytic theory of the pulse train and wakefield evolution in
long plasma channels, and compare the results with PIC simulations using WarpX
[17]. We use these results to optimize the performance of P-MoPAs, which are
set by both the dynamics of the accelerator stage, as well as the previously
derived constraints on the modulator stage [16]. We explore the key
differences between the pulse trains formed by P-MoPA and PBWA. We find that a
beneficial feature of the P-MoPA scheme is the ability to control the
longitudinal profile of the pulses in the train independent of their pulse
spacing; this allows for more resilient pulse trains and more efficient wake
excitation than PBWA, as well as offering a degree of control over the onset
of depletion effects.

Figure 1: [Color online]. Snapshot at $z=$2.1\text{\,}\mathrm{m}\mathrm{m}$$
of the relative plasma density perturbation $\delta n/n_{0}$ (blue-green)
resonantly excited by a P-MoPA pulse train with normalized intensity $|a|^{2}$
(red-yellow). The transverse electron density profile of the channel,
corresponding to Eq. (III) (green), and the longitudinally-averaged transverse
intensity profile (purple) of the drive pulses are shown on the end panel of
the plot. Results are from a 2D simulation using the PIC code WarpX [17] with
$W_{\text{drive}}=$1.2\text{\,}\mathrm{J}$$,
$\tau_{\text{drive}}=$1\text{\,}\mathrm{p}\mathrm{s}$$,
$\lambda_{L}=$1030\text{\,}\mathrm{n}\mathrm{m}$$,
$R=$30\text{\,}\mathrm{\SIUnitSymbolMicro}\mathrm{m}$$,
$w_{0}=$30\text{\,}\mathrm{\SIUnitSymbolMicro}\mathrm{m}$$,
$n_{00}=$2.5\text{\times}{10}^{17}\text{\,}\mathrm{c}\mathrm{m}^{-3}$$ and
modulator parameter $\beta=1.2$ as defined in Eq. (IV).

##  II Scaling Laws

Since the P-MoPA accelerator stage necessarily operates in the (quasi-)linear
regime, we expect the energy gain scaling laws for a P-MoPA to be similar to
those of the single pulse linear regime [18]. In this section we outline these
scaling laws, taking into account the effect of finite spot size effects on
the laser group velocity.

To extend the acceleration length beyond the Rayleigh range, the pulse train
is guided in a pre-formed plasma channel. A parabolic plasma channel matched
to a Gaussian pulse of spot size $w_{0}$ has a transverse electron density
profile given by [19, 20, 21]

$\displaystyle n_{0}(r)=n_{00}+\Delta n(r/w_{0})^{2}$ (1)

where $n_{00}$ is the on-axis plasma density and $\Delta n=(\pi
r_{e}w_{0}^{2})^{-1}$ is the channel depth parameter with $r_{e}$ being the
classical electron radius. The corresponding group velocity of a Gaussian
pulse propagating in such a matched channel, including finite spot size
effects, is given by [21]

$\displaystyle\frac{v_{g}}{c}=1-\frac{\omega_{p0}^{2}}{2\omega_{L}^{2}}-\frac{2c^{2}}{\omega_{L}^{2}w_{0}^{2}}\,,$
(2)

where $\omega_{p0}=\omega_{p}(r=0)$ indicates the on-axis plasma frequency. We
note that the contributions to the reduction in $v_{g}$ of the finite spot
size and plasma density are equal when $w_{0}=2c/\omega_{p0}$, which for an
on-axis plasma density of
$n_{00}=$1.0\text{\times}{10}^{17}\text{\,}\mathrm{c}\mathrm{m}^{-3}$$
corresponds to $w_{0}=$34\text{\,}\mathrm{\SIUnitSymbolMicro}\mathrm{m}$$.

From Eq. (2), the dephasing length, defined as the propagation distance over
which an electron propagating at $c$ slips relative to the wakefield by one-
quarter plasma wavelength $\sim\lambda_{p0}/4$ (assuming no channel effects on
wake structure), is

$\displaystyle
L_{\text{d}}=\frac{\lambda_{p0}}{4}\left(\frac{\omega_{p0}^{2}}{2\omega_{L}^{2}}+\frac{2c^{2}}{\omega_{L}^{2}w_{0}^{2}}\right)^{-1}.$
(3)

Note that a method has been proposed [22] which allows MP-LWFAs to accelerate
beyond the dephasing length by introducing a short region of higher plasma
density, which lets the accelerated electron bunch to rephase into the next
accelerating bucket back at the resonant density without deceleration.

Integrating the equations of motion for a relativistic particle moving in a
sinusoidal wakefield gives the energy gain $\Delta W_{e}$ before the onset of
dephasing as

$\displaystyle\frac{\Delta
W_{e}}{m_{e}c^{2}}\approx\frac{2}{\pi}\frac{eE_{\text{max}}L_{\text{d}}}{m_{e}c^{2}}=2\,\frac{E_{\text{max}}}{E_{\text{wb}}}\frac{\omega_{L}^{2}}{\omega_{p0}^{2}+4c^{2}/w_{0}^{2}}$
(4)

where $E_{\text{max}}$ is the peak accelerating field and
$E_{\text{wb}}=m_{e}\omega_{p0}c/e$ is the cold wavebreaking field. This shows
that for a given spot size, the maximum energy gain cannot be increased
indefinitely by reducing the plasma density as the dephasing length becomes
dominated by finite spot size effects. This is a well known effect in LWFA
[23], where the prescription is to scale the spot size with the plasma density
such that

$\displaystyle w_{0}\gtrsim\frac{2c}{\omega_{p0}}=\frac{\lambda_{p0}}{\pi}\,.$
(5)

On top of preventing an excessive reduction in group velocity, scaling the
spot size this way also prevents the wakefield from becoming largely radial
[24, 23]. As the required spot size increases with the plasma wavelength, the
total laser energy must necessarily increase. The total energy of the pulse
train $W_{\text{drive}}$ required to excite a given accelerating field
relative to the wavebreaking field $E_{\text{max}}/E_{\text{wb}}$ scales with
the volume enclosed by the spot size and the plasma length, i.e.
$W_{\text{drive}}\sim w_{0}^{2}\lambda_{p0}\sim\lambda_{p0}^{3}$. This grows
faster with the plasma wavelength than the electron energy gain $\Delta
W_{e}\sim\lambda_{p0}^{2}$ and the maximum bunch charge $Q_{e}\sim
n_{00}w_{0}^{2}\lambda_{p0}\sim\lambda_{p0}$, but the overall laser-to-bunch
energy efficiency $Q_{e}\Delta W_{e}/eW_{\text{drive}}$ has no plasma density
scaling. These scaling laws, which we have derived for pulse trains in the
linear regime, are identical to that of a linear single pulse LWFA [18].

Another limit to acceleration is pump depletion. Similar to the dephasing
length, the depletion length will also be reduced by finite spot size effects.
This is due to the extra pulse energy spent exciting the aforementioned
transverse wakefield for pulses with spot sizes comparable to the plasma
wavelength. The depletion length can be estimated by comparing the laser
energy with the energy contained within the linear wakefield and its
transverse counterpart spanning the depletion length $L_{\text{pd}}$. Assuming
that the pulse train is resonantly spaced and is comprised of individual
linearly-polarized Gaussian pulses of the form
$a_{m}^{2}=a_{0,m}^{2}\exp(-(\xi-\xi_{0,m})^{2}/L^{2})$, where $m=1,2,\ldots
N$, and that each pulse is short (i.e. $L\ll\lambda_{p0}$) so that each pulse
increases the normalized wake amplitude [25] by
$\sqrt{\pi/2}\,a_{0,m}^{2}k_{p0}L/4$, the pump depletion length is given by

$\displaystyle
L_{\text{pd}}=\left(\frac{\pi}{8}\frac{E_{\text{max}}}{E_{\text{wb}}}\right)^{-1}L_{\text{d}}\,.$
(6)

Note that this derivation is identical to that for single pulse LWFA. This
expression shows that the depletion length is always larger than the dephasing
length in the linear wakefield regime, which is also true for single pulse
LWFA [18]. To make efficient use of the available laser energy, we would
ideally want the pump depletion length to be larger but similar to the
dephasing length. As this is not feasible, the P-MoPA accelerator stage should
at least be operated in the quasi-linear regime, whilst avoiding nonlinear
detuning effects associated with the large amplitude wake which will be
discussed later. We also note that this treatment of the depletion length is
too simplistic for pulse trains since, as we will show later, it does not
capture the different rates at which each of the pulses deplete.

##  III Channel Effects on Energy Gain

Driving the accelerator within a pre-formed plasma channel causes the laser to
remain focused well beyond the Rayleigh range. However, the transverse
variation of plasma density within the channel leads to a concomitant
variation of the plasma frequency that can also influence the structure of the
excited wakefield and hence the energy gain of an injected electron bunch.

It is well known for single-pulse LWFA that plasma channels can strongly
influence wakefield structure via the curvature of wake phase-fronts [26]. For
example, shallow parabolic channels matched to the spot size $w_{0}$, i.e.
$\frac{1}{2}k_{p0}^{2}w_{0}^{2}\gg 1$, increase the overlap between the
accelerating and focusing regions of the wakefield by up to a factor of two at
sufficient distances behind the laser pulse. This can be beneficial for
acceleration as it can increase the dephasing length by a factor of up to two
without sacrificing much accelerating gradient. However, deep matched
parabolic channels, i.e. those for which $\frac{1}{2}k_{p0}^{2}w_{0}^{2}<1$,
yield a rapid conversion of axial field into radial field with distance behind
the drive pulse, meaning that LWFA is only possible in the first few buckets
behind the drive pulse. This quick conversion of axial field into radial field
readily outpaces other known wakefield decay mechanisms such as the ion motion
modulational instability [27] which only becomes relevant on the ion plasma
frequency timescale, corresponding to about forty plasma wavelengths for
hydrogen. Therefore the transverse profile of the plasma channel is the
dominant effect determining damping of the accelerating wakefield, and hence
the useful number of pulses in the train.

When using a single pulse driver, the channel effect on wake structure only
needs to be considered for one plasma wavelength behind the driver. However,
for a pulse train comprised of $N$ pulses, the channel influence on wake
structure is relevant over $N$ plasma wavelengths. The depth of a parabolic
channel will therefore determine the maximum useful number of pulses used in a
MP-LWFA, as the axial wakefield contribution of the first pulse in a train
will disappear by the time the $N^{\text{th}}$ pulse passes. There are two
ways to allow for longer trains. The first way is to operate well in the
shallow channel limit $\frac{1}{2}k_{p0}^{2}w_{0}^{2}\gg N$. However, this
forces the required matched spot size to become large at low plasma densities
and large pulse train lengths (e.g.
$w_{0}\gg$75\text{\,}\mathrm{\SIUnitSymbolMicro}\mathrm{m}$$ at
$n_{00}=$1.0\text{\times}{10}^{17}\text{\,}\mathrm{c}\mathrm{m}^{-3}$,N=10$),
bringing the required laser energy to drive the wakefield in the quasi-linear
regime to the 10-joule scale for a $1030\text{\,}\mathrm{n}\mathrm{m}$ laser
wavelength. The other solution is to use square-profile channels, which have a
much weaker effect on the wakefield structure. Fortunately, the transverse
density profiles of CHOFI plasma waveguides are closer to square than
parabolic in shape [28, 29, 30], which makes them well suited to driving
P-MoPAs.

An example PIC simulation of the wake structure driven by a P-MoPA pulse train
is presented in Figure 1, where the pulse train is guided in a quasi-square
channel of the form

$\displaystyle\frac{n_{0}(r)-n_{00}}{\Delta n}=$
$\displaystyle\begin{cases}(r/R)^{10}&r<1.2R\\\ (1.2)^{10}&1.2R\leq
r<1.2R+d\\\ (1.2)^{10}\left(1-\frac{r-1.2R-d}{d}\right)&1.2R+d\leq
r<1.2R+2d\\\ 0&\text{otherwise}\end{cases}$ (7)

where $R$, $\Delta n=(\pi r_{e}R^{2})^{-1}$ and
$d=$10\text{\,}\mathrm{\SIUnitSymbolMicro}\mathrm{m}$$ are the channel wall
radius, depth and thickness parameters respectively. Resonant excitation of a
large amplitude $\sim 50\%$ wakefield on axis is evident. Also apparent is a
large amplitude wakefield on the inner wall of the channel. This is due to the
steep plasma density gradient

$\displaystyle\frac{\mathrm{d}\ln n_{0}}{\mathrm{d}r}>k_{p}\,,$ (8)

which is characteristic of plasma channels at these spot sizes and low $\sim
10^{17}$ cm-3 on-axis densities. The wake amplitude in this region continues
to grow even after the laser driver has already passed due to the
aforementioned conversion of the wakefield from the longitudinal to transverse
component [26]. This accumulated wall-boundary wake can eventually exceed the
on-axis wake amplitude, suggesting that large-amplitude wakefields driven by
pulse trains in plasma channels would undergo transverse wavebreaking effects
before longitudinal wavebreaking.

##  IV P-MoPA Pulse Shaping

The plasma modulator offers a degree of control of the temporal profile of the
pulses in the pulse train, allowing for further optimization of the
accelerator stage. The pulse train formed after compression will be comprised
of pulses evenly spaced by the modulator plasma period with pulse heights and
durations determined by the amount of spectral modulation introduced by the
modulator stage. The additional bandwidth of the spectrally-modulated drive
pulse enables a P-MoPA to form pulses substantially shorter than their
separation, unlike the cosine-squared pulses created with the beatwave method.
This allows larger wake amplitudes to be driven in a P-MoPA with the same
drive pulse energy.

According to 3D spectral modulation theory [16], assuming that the spectral
modulation remains primarily in the fundamental channel mode, the pulse trains
formed using the “ideal compressor” described by Jakobsson et al [9], with a
seed wake of the form $\delta n(r,\xi)=\delta
n_{s}\cos(k_{p0}\xi)\exp(-2r^{2}/w_{0}^{2})$, will have the following temporal
amplitude envelope

$\displaystyle
f_{\text{acc}}(\xi;\beta)=f_{\text{mod}}(\xi)\left(J_{0}(\beta)+2\sum_{n=1}^{\infty}J_{n}(\beta)\cos(nk_{p0}\xi)\right)\,,$
$\displaystyle\beta=\frac{2\Omega_{s}L_{\text{mod}}}{v_{g,\text{mod}}}$ (9)

where $\xi=z-ct$ is the longitudinal co-moving coordinate, $0\leq
f_{\text{mod}}(\xi)\leq 1$ is the input temporal amplitude envelope,
$L_{\text{mod}}$ is the modulator length, $v_{g,\text{mod}}$ is laser group
velocity, and $\Omega_{s}=(\omega_{p0}^{2}/8\omega_{L})(\delta n_{s}/n_{00})$
is the rate of spectral modulation of a Gaussian pulse co-propagating with a
seed wake of on-axis amplitude $\delta n_{s}$ in a parabolic channel [16]. The
modulator parameter $\beta$ determines the effective number of sidebands
generated in the drive pulse spectrum and controls the temporal profile of the
pulse train produced by the plasma modulator (whilst obeying the modulator
stability condition outlined in [16]). Figure 2 shows how the temporal
intensity profiles of the pulses within the train, their full-length at half
maximum duration, and their contrast, varies with $\beta$. It can be seen that
adjusting the parameter $\beta$ allows the duration and contrast of the pulses
within the train to be controlled.

The peak intensity occurs at $\beta=j_{0,1}\approx 2.405$, where $j_{m,n}$
indicates the $n^{\text{th}}$ nonzero root of the Bessel function $J_{m}(x)$.
This can be shown by evaluating the temporal envelope at the pulse center
using well-known identities for Neumann series of Bessel functions [31, 32]

$\displaystyle\frac{f_{\text{acc}}(0;\beta)}{f_{\text{mod}}(0)}$
$\displaystyle=J_{0}(\beta)+2\sum_{n=1}^{\infty}J_{n}(\beta)$
$\displaystyle=1+\int_{0}^{\beta}J_{0}(s)ds$ (10)

which shows that its global maximum occurs at $\beta=j_{0,1}$, where it takes
the value $f_{\text{acc}}\approx 2.470\,f_{\text{mod}}$. As evident in Figure
2, increasing $\beta$ above this value leads to the formation of shorter
pulses, but these have temporal wings, which wastes energy since they do not
contribute to wake excitation. Thus when optimizing the P-MoPA scheme we only
need to consider $\beta\leq j_{0,1}$.

The additional control that the parameter $\beta$ provides enables more
efficient wake excitation than is possible with PBWA. This is shown in Figure
2(b), which shows that $\beta$ in the range of 1.1–1.6 optimizes the wake
amplitude according to linear wakefield theory [18]. We can also show in the
linear regime that the ratio between the amplitude of the wake potential
driven by a P-MoPA pulse train and an ideal beatwave pulse train with a
temporal intensity modulation of the form $1+\cos(k_{p0}\xi)$ containing the
same laser energy is given by

$\displaystyle\frac{\phi_{\text{P-MoPA}}}{\phi_{\text{PBWA}}}$
$\displaystyle=\frac{\int_{-\pi}^{\pi}\left[J_{0}(\beta)+2\sum_{n=1}^{\infty}J_{n}(\beta)\cos(ns)\right]^{2}\cos(s)ds}{\int_{-\pi}^{\pi}\left[1+\cos(s)\right]\cos(s)ds}$
$\displaystyle=4\sum_{n=0}^{\infty}J_{n}(\beta)J_{n+1}(\beta)$
$\displaystyle=2\beta\left[J_{0}^{2}(\beta)+J_{1}^{2}(\beta)\right]\,.$ (11)

This ratio is maximized for $\beta\approx 1.43$, at which value the wake
amplitude in a P-MoPA is $72\%$ larger than that driven by an ideal PBWA with
the same total drive energy.

(a)  (b)

Figure 2: [Color online]. (a) Pulse train temporal profiles for various values
of the modulator parameter $\beta$ and their respective FWHM durations
(dashed). (b) FWHM duration (blue, solid), peak prominence (orange, dashed)
and the wakefield amplitude according to 1D linear theory (green, dot-dashed,
arb. units) as a function of $\beta$.

##  V The Rosenbluth-Liu Limit

Resonant wakefield acceleration requires the driving pulses to be spaced by
the plasma period. However, as the plasma wave amplitude increases beyond the
linear regime, its wavelength increases due to relativistic effects. The
consequent build-up of plasma wave phase difference relative to the laser
pulse centroids sets a limit on the maximum wake amplitude achievable.

The saturation effect arising from this detuning was first derived by
Rosenbluth and Liu [2] in the context of plasma beatwave excitation. In their
derivation, two infinite-duration pumps with normalized vector potential
amplitudes $a_{1}$ and $a_{2}$ interfere to form a cosine beatwave pulse train
resonant with the plasma, which results in the saturated wakefield amplitude
scaling with the third root of the intensity of the beatwave modulation

$\displaystyle\frac{E_{\text{sat}}}{E_{\text{wb}}}=\left(\frac{16a_{1}a_{2}}{3}\right)^{1/3}.$
(12)

P-MoPA pulse trains, which are also spaced uniformly by design, are thus
expected to exhibit a similar limit. However, for a P-MoPA, an important
difference is that the train will comprise a relatively small number of
pulses.

According to 1D nonlinear plasma wave theory [18], the nonlinear correction to
the plasma wavelength for a plasma wave of amplitude
$\phi_{a}=E_{a}/E_{\text{wb}}$ at non-relativistic laser intensities in the
$\phi_{a}^{2}\ll 1$ limit is given by

$\displaystyle\lambda_{p}^{\text{NL}}\approx\lambda_{p}(1+3\phi_{a}^{2}/16)\,,$
$\displaystyle k_{p}^{\text{NL}}\approx k_{p}(1-3\phi_{a}^{2}/16)\,.$ (13)

Resonant excitation of the wakefield eventually saturates due to detuning
caused by the progressive growth of the plasma wavelength. The co-moving
coordinate $\xi_{\text{sat}}$ at which this detuning occurs is approximately
given by

$\displaystyle\int_{\xi_{\text{sat}}}^{\infty}(k_{p}-k_{p}^{\text{NL}})d\xi\approx\frac{\pi}{2}\,.$
(14)

Since the pulse train in a P-MoPA pulse train corresponds to a sequence of
short pulses spaced by $T_{p}$, within an envelope determined by the duration
of the ps-duration drive pulse, most of its energy will be contained within
only a few pulses close to its center. To account for this we model the
original pulse envelope as a temporally Gaussian pulse of the form

$\displaystyle
f_{\text{mod}}^{2}(\xi)=\exp\left(-\frac{\xi^{2}}{2\sigma^{2}}\right)\,,$ (15)

where $\sigma=c\tau_{\text{drive}}/2\sqrt{2\ln 2}$ with $\tau_{\text{drive}}$
being the FWHM duration of the drive pulse. From linear theory, ignoring
detuning, the wakefield potential amplitude $\phi_{a}(\xi)$ excited by the
pulse train is given by

$\displaystyle\frac{\phi_{a}(\xi)}{\phi_{\text{max}}}=\frac{\int_{\xi}^{\infty}d\xi^{\prime}f_{\text{acc}}^{2}(\xi^{\prime})}{\int_{-\infty}^{\infty}d\xi^{\prime}f_{\text{acc}}^{2}(\xi^{\prime})}\approx\frac{1}{2}\,\text{erfc}\left(\frac{\xi}{\sqrt{2}\,\sigma}\right)\,,$
(16)

where $\phi_{\text{max}}=E_{\text{max}}/E_{\text{wb}}$ is the maximum
wakefield amplitude resonantly excited by the pulse train. For the pulse train
to efficiently drive a wakefield, resonant detuning must not occur before the
majority of the pulses have already passed. If we demand that at least $\sim
90\%$ of the pulse train energy passes before the full $\pi/2$ detuning
occurs, then combining Eqs. V, 14 and 16 results in the following condition

$\displaystyle-\frac{\xi_{\text{sat}}^{90\%}}{\lambda_{p}}>0.27\,N_{\text{eff}}\,,$
(17)

where $N_{\text{eff}}=2\tau_{\text{FWHM}}/T_{p}$ is defined as the effective
number of pulses in the pulse train. Substituting this expression along with
Eq. (16) into Eq. (14) yields the resonant detuning limit for wakefields
driven by P-MoPA pulse trains:

$\displaystyle\phi_{\text{max}}^{90\%}=\frac{E^{90\%}_{\text{max}}}{E_{\text{wb}}}<\frac{2.8}{\sqrt{N_{\text{eff}}}}\,.$
(18)

For example, for $\tau_{\text{FWHM}}=$2\text{\,}\mathrm{p}\mathrm{s}$$ at
$n_{00}=$2.5\text{\times}{10}^{17}\text{\,}\mathrm{c}\mathrm{m}^{-3}$$, a
P-MoPA would form a pulse train comprised of $N_{\text{eff}}=18$ pulses.
Substituting this into Eq. (18) then limits the maximum accelerating gradient
before detuning to
$E_{\text{max}}<$32\text{\,}\mathrm{G}\mathrm{V}\mathrm{/}\mathrm{m}$$, or
$<66\%$ of the cold wavebreaking field. If we operate the P-MoPA in the linear
and quasi-linear regime, for which $\phi^{90\%}_{\text{max}}\lesssim 0.5$,
then Eq. (18) becomes $N_{\text{eff}}<30$. We therefore conclude that
wakefield saturation caused by this detuning mechanism will not be an issue
for a drive train comprising $\sim 30$ or fewer pulses.

A further consideration for P-MoPA pulse trains is that the individual pulses
can be much shorter than in PBWA, such that $a\ll 1$ can no longer be assumed.
Thus, to maximize the wake amplitude, the relativistic effects associated with
the quiver motion in regions of high intensity, rather than the relativistic
effects associated with the fluid motion, should also be corrected for by
using a slightly higher plasma density in the accelerator than the modulator
stage to counteract this detuning mechanism.

Aside from limiting the wakefield amplitude, resonant detuning can also lead
to focusing/defocusing effects caused by pulses co-propagating with high
amplitude regions of the wakefield. These effects are discussed in the next
section.

##  VI Relativistic Self-Phase Modulation and Transverse Mode Excitation

Up until now we have been treating the pulse train envelope as unchanging as
it propagates along the accelerator. While this assumption is useful for
deriving scaling laws and isolating certain physics, in reality the pulse
train is a dynamic 3D structure which will evolve under the influence of both
relativistic nonlinearities and its own excited wake. Neglecting depletion for
now, the two main effects to consider for MP-LWFA in long channels are
relativistic self-phase modulation (SPM) and tranvserse mode dynamics caused
by the excitation of higher order transverse channel modes by the
focusing/defocusing effects of the large amplitude wake.

Relativistic SPM causes further red/blue-shifting near the head/tail of each
pulse, especially for the shortest duration, highest $a_{0}$ pulses.
Transverse mode excitation results in spot size oscillations and hence
oscillations in the wakefield structure, which become more pronounced towards
the back of the pulse train where the wake amplitude is the greatest. These
oscillations could negatively affect beam quality, since a propagation-varying
focusing field could increase the energy spread and emittance of the
accelerated bunch.

Before significant pump depletion has occurred, i.e. while the frequencies
within the laser spectrum satisfy $|\omega-\omega_{L}|/\omega_{L}\ll 1$, the
evolution of the pulse train propagating in an axisymmetric channel of the
form $n_{0}(r)=n_{00}+\delta n_{0}(r)$ can be described using the paraxial
wave equation in the weakly relativistic limit [33, 34, 35, 36]

$\displaystyle\left[\frac{i}{\omega_{L}}\frac{\partial}{\partial\tau}+\frac{c^{2}}{2\omega_{L}^{2}}\Delta_{\perp}\right]a=$
$\displaystyle\frac{\omega_{p}^{2}}{2\omega_{L}^{2}n_{0}}\left[\delta
n_{0}(r)+\delta n(r,\xi;|a|^{2})-n_{0}(r)|a|^{2}/4\right]a$ (19)

where $a(r,\theta,\xi,\tau)$ is the envelope of the normalized vector
potential of the pulse, $\omega_{L}$ is the laser frequency and the
propagation is described in co-moving coordinates $\xi=z-v_{g0}t$, $\tau=t$,
with $v_{g0}/c=(1-\omega_{p0}^{2}/\omega_{L}^{2})^{1/2}$ defined as the group
velocity of electromagnetic plane waves in uniform plasma of density $n_{00}$,
corresponding to the on-axis plasma channel frequency $\omega_{p0}$, and
$\Delta_{\perp}=\partial^{2}/\partial r^{2}+(1/r)\partial/\partial
r+(1/r^{2})\partial^{2}/\partial\theta^{2}$ is the tranvserse Laplacian.

Figure 3: [Color online]. Results of a 2D PIC simulation of the propagation of
a pulse train in a quasi-square plasma channel. For this simulation, $\xi$ as
defined after Eq. (VI) at varying propagation distances $z$ in a quasi-square
channel of the form given in Eq. (III) with
$n_{00}=$2.5\text{\times}{10}^{17}\text{\,}\mathrm{c}\mathrm{m}^{-3}$$,
$\Delta n=$1.26\text{\times}{10}^{17}\text{\,}\mathrm{c}\mathrm{m}^{-3}$$,
$R=$30\text{\,}\mathrm{\SIUnitSymbolMicro}\mathrm{m}$$,
$W_{\text{drive}}=$2.4\text{\,}\mathrm{J}$$,
$\tau_{\text{drive}}=$1\text{\,}\mathrm{p}\mathrm{s}$$ and modulator parameter
$\beta=1.2$. The pulse train at the entrance of the channel initially has a
Gaussian transverse profile with spot size
$w_{0}=$30\text{\,}\mathrm{\SIUnitSymbolMicro}\mathrm{m}$$. For each plot, the
top panel displays the laser intensity profile $|a|^{2}$ along with the local
longitudinal variation of the laser centroid $x_{c}(\xi)$ (dashed) and local
effective spot size $w_{\text{eff}}(\xi)$ (solid). The bottom panel displays
the wake density perturbation $\delta n/n_{00}$.

Pre-formed plasma channels have often been modelled using parabolic transverse
profiles, since these guide Laguerre-Gaussian modes. The long
($\gtrsim$100\text{\,}\mathrm{m}\mathrm{m}$$), low-attenuation and kHz
production rate plasma channels required for this scheme have been
experimentally demonstrated using the CHOFI technique [28, 29, 30], which
forms deep, low-loss channels with on-axis densities as low as
$1\text{\times}{10}^{17}\text{\,}\mathrm{c}\mathrm{m}^{-3}$. These CHOFI
channels are often better modelled with transverse profiles closer to square
rather than parabolic [29], meaning that wake phase-front curvature effects
can be neglected. In the following analysis we will use a simplified profile
and treat the plasma channel as an infinite square well of the form

$\displaystyle n_{0}(r)=\begin{cases}n_{00}&r\leq R_{\infty}\\\
\infty&r>R_{\infty}\end{cases}\,,$ (20)

where $R_{\infty}$ indicates the wall radius parameter specifically for this
infinite square well. Note that the analysis that follows can be generalized
to any axisymmetric, monotonically increasing channel profile. Assuming that
the channel is unperturbed and that relativistic effects are negligible,
substituting Eq. (20) into Eq. (VI) yields the following orthogonal set of
Bessel mode solutions for $r\leq R_{\infty}$

$\displaystyle
a_{pm}(r,\theta,\xi,\tau)=\alpha_{pm}(\xi)\mathcal{J}_{pm}(r,\theta)e^{-i\Omega_{pm}\tau}\,,$
$\displaystyle\mathcal{J}_{pm}(r,\theta)=\frac{J_{1}(j_{0,1})}{\left|J_{|m|+1}(j_{|m|,p+1})\right|}J_{|m|}\left(j_{|m|,p+1}\frac{r}{R_{\infty}}\right)e^{im\theta}\,,$
$\displaystyle\Omega_{pm}=\frac{j_{|m|,p+1}^{2}c^{2}}{2\omega_{L}R_{\infty}^{2}}\,,$
(21)

where the integers $p\geq 0$ and $m$ are the radial and azimuthal indexes
respectively. Note that the effective spot size $w_{\text{0,eff}}$ of the
fundamental mode is related to the wall radius $R_{\infty}$ by

$\displaystyle
w_{\text{0,eff}}=\sqrt{\frac{\int_{0}^{R}2r^{2}\mathcal{J}_{00}^{2}(r)rdr}{\int_{0}^{R}\mathcal{J}_{00}^{2}(r)rdr}}\approx
0.66R_{\infty}\,.$ (22)

which is analogous to the result for the standard step-index waveguide in the
limit of infinite $V$-number [37]. Similar to a previous analysis of the
P-MoPA modulator stage [16], long-propagation of pulse trains in plasma
channels are expected to excite higher order transverse channel modes, with an
amplitude that scales with the wake amplitude. To model this effect, we first
assume a fixed multipulse-driven wake of the form

$\displaystyle\delta n(r,\xi)=\delta n_{a}(\xi)\mathcal{J}_{00}^{2}(r)\,,$
(23)

which is consistent with a pulse train primarily in the fundamental mode
exciting a wake in the (quasi-)linear regime. Applying time-dependent
perturbation theory (TDPT) to Eq. (VI) for pulses propagating in a square
channel with only the $m=0$ radial transverse modes present, we find the
following set of PDEs for the coefficients $\alpha_{p}(\xi,\tau)$ of each of
the square channel’s radial modes $\mathcal{J}_{p0}(r)$

$\displaystyle i\frac{\partial\alpha_{p}(\xi,\tau)}{\partial\tau}=$
$\displaystyle\frac{\omega_{p0}^{2}}{2\omega_{L}}\sum_{n}\alpha_{n}(\xi,\tau)\Big{\langle}\mathcal{J}_{p0}\Big{|}\frac{\delta
n}{n_{00}}-\frac{|a|^{2}}{4}\Big{|}\mathcal{J}_{n0}\Big{\rangle}e^{i\left(\Omega_{p0}-\Omega_{n0}\right)\tau}\,,$
$\displaystyle
a(r,\xi,\tau)=\sum_{p}\alpha_{p}(\xi,\tau)\mathcal{J}_{p0}(r)e^{-i\Omega_{p0}\tau}\,,$
$\displaystyle\langle\mathcal{J}_{p^{\prime}m^{\prime}}|(...)|\mathcal{J}_{pm}\rangle\equiv\frac{\int_{0}^{2\pi}d\theta\int_{0}^{R}rdr\mathcal{J}_{p^{\prime}m^{\prime}}^{\ast}(...)\mathcal{J}_{pm}}{J_{1}^{2}(j_{0,1})\pi
R_{\infty}^{2}}\,.$ (24)

We note that TDPT requires that $\pi r_{e}(2R_{\infty}/j_{0,1})^{2}\delta n\ll
1$, which is _not_ true everywhere in the accelerator stage due to the large
amplitude wake. However, for a resonant pulse train the zeros of $\delta
n_{a}(\xi)$ are located near the laser centroids, so TDPT _is_ still valid
provided that the duration of each pulse is short, so that $\delta n$ is
sufficiently small in regions of high intensity. Hence we can still use TDPT
to describe the evolution of the pulse train despite the large amplitude wake.

Assuming that most of the light is initially in the fundamental mode,
transverse mode transitions will be dominated by excitations of only the first
radial mode, with higher order mode excitations being much weaker. Hence
substituting our expression for the wake and approximating Eq. (VI) as a two-
level system description of the fundamental and first radial transverse modes
yields

$\displaystyle\frac{\partial}{\partial\tau}\begin{pmatrix}\alpha_{0}\\\
\alpha_{1}\end{pmatrix}=\bm{D}(\xi,\tau)\begin{pmatrix}\alpha_{0}\\\
\alpha_{1}\end{pmatrix}\,,$
$\displaystyle\bm{D}(\xi,\tau)=-i\gamma(\xi)\begin{pmatrix}d_{00}&d_{01}e^{-i\omega_{w}\tau}\\\
c.c.&d_{11}\end{pmatrix}\,,$
$\displaystyle\gamma(\xi)=\frac{\omega_{p0}^{2}}{2\omega_{L}}\left(\frac{\delta
n_{a}(\xi)}{n_{00}}-\frac{|a_{0}|^{2}f_{\text{acc}}^{2}(\xi)}{4}\right)\,,$
$\displaystyle\omega_{w}=\Omega_{10}-\Omega_{00}\approx\frac{12.3\,c^{2}}{\omega_{L}R_{\infty}^{2}}\,$
$\displaystyle
d_{00}=\langle\mathcal{J}_{00}|\mathcal{J}_{00}^{2}|\mathcal{J}_{00}\rangle\approx
0.57\,,\,\,\,d_{01}=\langle\mathcal{J}_{00}|\mathcal{J}_{00}^{2}|\mathcal{J}_{10}\rangle\approx
0.27\,,$ $\displaystyle
d_{11}=\langle\mathcal{J}_{10}|\mathcal{J}_{00}^{2}|\mathcal{J}_{10}\rangle\approx
0.48$ (25)

where $\omega_{w}$ is the spot size oscillation frequency, $c.c.$ indicates
the complex conjugate of the upper anti-diagonal entry, and we have assumed
the fundamental intensity profile
$|a|^{2}\approx|a_{0}|^{2}f_{\text{acc}}^{2}(\xi)\mathcal{J}^{2}_{00}(r)$ for
the weakly relativistic term. This linear system of PDEs can then be solved
analytically

$\displaystyle\begin{pmatrix}\alpha_{0}(\xi,\tau)\\\
\alpha_{1}(\xi,\tau)\end{pmatrix}=\bm{A}(\xi,\tau)\begin{pmatrix}\alpha_{0}(\xi,0)\\\
\alpha_{1}(\xi,0)\end{pmatrix}\,,$
$\displaystyle\bm{A}(\xi,\tau)=\exp\left[-i\gamma(\xi)\tau\begin{pmatrix}d_{00}&d_{01}\text{sinc}\left(\frac{\omega_{w}\tau}{2}\right)e^{-i\omega_{w}\tau/2}\\\
c.c.&d_{11}\end{pmatrix}\right]\,,$ $\displaystyle
w_{\text{eff}}(\xi,\tau)\approx
w_{0,\text{eff}}\left(1-0.40\,\frac{\alpha_{0}^{\ast}\alpha_{1}e^{-i\omega_{w}\tau}+c.c.}{|\alpha_{0}|^{2}}\right)\,.$
(26)

This solution describes the excitation of the first radial transverse mode due
to transverse (de-)focusing and self-modulation, both of which can arise
through the effects of the wake and self-focusing. The interference between
the fundamental and first radial mode then leads to a modulation of the local
effective spot size $w_{\text{eff}}(\xi,\tau)$. Evaluating this matrix
exponential yields a unitary matrix with diagonal entries of $\mathcal{O}(1)$
and anti-diagonal entries of $\mathcal{O}(\gamma/\omega_{w})\ll 1$. This
ensures that, as long as the original TDPT condition is satisfied in regions
of high intensity, laser light initially in the fundamental transverse channel
mode will also remain primarily in the fundamental mode, with small
$\mathcal{O}(\gamma^{2}/\omega_{w}^{2})$ and $\mathcal{O}(|\delta
w|^{2}/w_{0}^{2})$ oscillating exchanges in energy between the fundamental and
first radial modes, where $\delta w/w_{0}\ll 1$ is some small initial spot
size mismatch.

Since the bulk of the wake is already excited by the zeroth order intensity
$\sim\mathcal{J}_{00}^{2}(r)$, we can neglect the wake contribution of the
$\mathcal{O}(\gamma/\omega_{w})\ll 1$,
$\sim\mathcal{J}_{00}(r)\mathcal{J}_{10}(r)$ first order correction to the
laser transverse mode dynamics. However, the near-axis wakefield itself can be
significantly affected as laser energy continuously oscillates towards and
away from the axis, resulting in propagation-dependent accelerating and
focusing fields which could, in principle, impact both the emittance and
energy spread of an injected electron bunch.

An example of this phenomenon in 2D PIC simulations is evident in Figure 3,
which shows snapshots of the laser intensity envelope and density perturbation
as it propagates in a quasi-square channel. We can see that transverse mode
excitation is the strongest for the trailing pulses which co-propagate with
the largest amplitude wake, as described by Eq. (VI). This results in an
oscillation of the wake structure and hence the maximum on-axis accelerating
field $E_{\text{max}}$. (This oscillation of on-axis field can also be
observed directly in Figure 5, as we discuss below.) Despite the transverse
mode excitation causing a break-up of the trailing pulses, Figures 3 and 4
show that the pulse train is still able to resonantly excite a wakefield even
in the quasi-linear regime with $E_{\text{max}}/E_{\text{wb}}\sim 0.5$.

To compare the spot size oscillation frequency $\omega_{w}$ for a square
channel given in Eq. (VI) with the PIC simulations presented in Figure 5, a
correction has to be made for the fact that the plasma channel used in all of
the presented PIC simulations is quasi-square with walls of finite height and
width of the form given in Eq. (III). For this quasi-square channel we define
an effective channel radius, $R_{\infty}^{\text{eff}}$, as that radius for
which the fundamental mode, $\mathcal{J}_{00}(r)$, of an infinite square
channel best matches the lowest-order mode of the quasi-square channel. For
the $R=$30\text{\,}\mathrm{\SIUnitSymbolMicro}\mathrm{m}$$ case, this results
in an effective channel radius of approximately
$R_{\infty}^{\text{eff}}\approx$46\text{\,}\mathrm{\SIUnitSymbolMicro}\mathrm{m}$$.
Substituting this into Eq. (VI) then gives a spot size oscillation wavelength
of $2\pi c/\omega_{w}\approx$6.6\text{\,}\mathrm{m}\mathrm{m}$$, which is in
good agreement with the observed oscillation period of the on-axis
accelerating field in Figures 5 (a) and (b).

Various strategies can be employed to minimize this effect. First, and most
obviously, the mode of the input pulse train should be matched to the lowest-
order mode of the channel as closely as possible, which would minimize
oscillations associated with $|\delta w|/w_{0}$. Next the following quantity,
which is a measure of the oscillation of the wakefield structure behind the
pulse train, should be minimized:
$\left|\left(\int_{-\infty}^{\infty}d\xi(\gamma(\xi)/\omega_{w})f^{2}_{\text{acc}}(\xi)\right)\big{/}\left(\int_{-\infty}^{\infty}d\xi
f^{2}_{\text{acc}}(\xi)\right)\right|$. For a given wake amplitude $\delta
n_{a}(\xi)/n_{00}$, the two ways to minimize this quantity are: (i) to reduce
detuning effects to ensure that the plasma wave zeros are as close to the
pulse centroids as possible; and (ii) to increase the modulator parameter
$\beta$, which shortens the duration of each pulse in the train, so that they
only sample the plasma wave close to its zero-crossings; although we note that
this latter strategy also increases SPM, which hastens the onset of depletion
effects towards the heads of the pulses.

##  VII Depletion Effects

We also have to consider that the paraxial description which led to the result
in Eq. (VI) breaks down when depletion becomes non-negligible, which becomes
important when driving quasi-linear wakes where the dephasing length is a
large fraction of the depletion length.

One caveat to the earlier analysis of the depletion length in Eq. (6) is that
pulse trains do not deplete uniformly. Instead, trailing pulses, which co-
propagate with steeper density gradients, will deplete faster than those near
the beginning of the pulse train. We thus expect the pulses towards the rear
of the train to have individual depletion lengths up to half as short as in
Eq. (6), whilst the pulses near the beginning of the train will experience
minimal depletion.

The non-uniform depletion effects include group velocity dispersion (GVD)
between the pulses, caused by the trailing depleted pulses slowing down more,
and detuning spot size oscillations as depleted pulses have shorter spot size
oscillation periods due to their shorter Rayleigh lengths.

As demonstrated in Figure 4, the pulse train is initially uniformly spaced and
can therefore suffer from Rosenbluth-Liu detuning, though it remains a small
effect due to the low number of pulses. However, as the pulse train depletes,
the strongly red-shifting trailing pulses fall behind and become resonant with
the wake. We also observe a small portion of the heads of each of the trailing
pulses strongly blue-shifting and catching up with the train.

The onset of depletion effects is also expedited by relativistic SPM, as this
increases the rate at which the heads of each pulse in the pulse train red-
shifts and causes pulse steepening [38], which is observed in Figure 4(a). On
one hand, depletion causes bunch dephasing to occur sooner than predicted by
Eq. (3) by a reduction in group velocity. On the other hand, the resulting
reduction in instantaneous laser frequency increases the ponderomotive force
(which scales with $F_{p}\propto\omega^{-1}$ when holding photon number
density constant), causing the wakefield amplitude to gradually increase with
propagation distance, similar to the effect observed experimentally by M. J.
V. Streeter et al. [39].

(a)  (b)

Figure 4: [Color online]. Results from the same PIC simulation presented in
Figure 3. (a) The on-axis normalized axial electric field
$E_{z}/E_{\text{wb}}$ (solid, black) and the normalized laser intensity
envelope $|E/E_{L}|^{2}$ with $E_{L}=m_{e}\omega_{L}c/e$ and its shift in
instantaneous frequency $\Delta\omega/\omega_{L}$ (color scale) after 2 mm
(top) and 100 mm (bottom) of propagation. (b) The evolution of the transverse-
averaged intensity envelope $\langle|E/E_{L}|^{2}\rangle_{\perp}$. Note that
$|E/E_{L}|\neq|a|$ at large propagation distances due to the significant red-
shifting of the trailing pulses.

##  VIII Multi-GeV P-MoPA electron accelerator

Figure 5 summarizes the performance of the accelerator stage of a P-MoPA for
several different drive laser, and modulator stage parameters. The parameters
selected take into account the constraints set by the modulator stage (which
sets a limit on $W_{\text{drive}}$ for a given modulator spot size and plasma
density [16]), relativistic detuning effects, optimizing the wake amplitude
(by using a suitable modulator parameter $\beta$) and the evolution with $z$
of the wakefield (caused by SPM and transverse mode excitation). The figures
show the evolution of the laser spectra, maximum on-axis accelerating
wakefield, and the energy of a test electron bunch. For these 2D PIC
simulations, a 1 pC, 35 MeV Gaussian electron bunch was manually placed in
phase with the peak accelerating wakefield with dimensions
$x_{\text{rms}}=$4\text{\,}\mathrm{\SIUnitSymbolMicro}\mathrm{m}$,\,z_{\text{rms}}=$1.5\text{\,}\mathrm{\SIUnitSymbolMicro}\mathrm{m}$$
for $n_{00}=$2.5\text{\times}{10}^{17}\text{\,}\mathrm{c}\mathrm{m}^{-3}$$.
The bunch dimensions were scaled up by a factor of $\sqrt{2.5}$ for the
$n_{00}=$1.0\text{\times}{10}^{17}\text{\,}\mathrm{c}\mathrm{m}^{-3}$$
simulation.

(a)  $W_{\text{drive}}=$2.4\text{\,}\mathrm{J}$$,
$R=$30\text{\,}\mathrm{\SIUnitSymbolMicro}\mathrm{m}$$,
$n_{00}=$2.5\text{\times}{10}^{17}\text{\,}\mathrm{c}\mathrm{m}^{-3}$$ (b)
$W_{\text{drive}}=$2.4\text{\,}\mathrm{J}$$,
$R=$30\text{\,}\mathrm{\SIUnitSymbolMicro}\mathrm{m}$$,
$n_{00}=$2.5\text{\times}{10}^{17}\text{\,}\mathrm{c}\mathrm{m}^{-3}$$ (c)
$W_{\text{drive}}=$5.0\text{\,}\mathrm{J}$$,
$R=$36\text{\,}\mathrm{\SIUnitSymbolMicro}\mathrm{m}$$,
$n_{00}=$1.0\text{\times}{10}^{17}\text{\,}\mathrm{c}\mathrm{m}^{-3}$$

Figure 5: [Color online]. 2D PIC simulation results of accelerator performance
at different modulation parameters $\beta$ in a quasi-square channel of the
same form used in Fig. 3. The top panels display the normalized laser spectrum
$|\hat{E}(\omega)|$ evolving with the propagation distance, whilst the bottom
panels show the evolution of the normalized injected bunch electron energy
spectrum $\hat{Q}_{e}(W_{e})$ and maximum on-axis accelerating field
$E_{\text{acc}}^{\text{max}}$. Note that the oscillation of the on-axis
accelerating field due to spot size oscillations is a continuous sinusoidal
motion which disappears beyond $50\text{\,}\mathrm{m}\mathrm{m}$ in (a,b) and
$130\text{\,}\mathrm{m}\mathrm{m}$ in (c), and only appears to be jagged in
these plots due to sub-sampling.

For $n_{00}=$2.5\text{\times}{10}^{17}\text{\,}\mathrm{c}\mathrm{m}^{-3}$$,
the observed propagation distance at which the maximum energy gain is achieved
for the $\beta=1.2$ case corresponds well with
$L_{\text{d}}=$93\text{\,}\mathrm{m}\mathrm{m}$$ predicted by Eq. (3) using
$R_{\infty}^{\text{eff}}=$46\text{\,}\mathrm{\SIUnitSymbolMicro}\mathrm{m}$\Rightarrow
w_{0,\text{eff}}=$30\text{\,}\mathrm{\SIUnitSymbolMicro}\mathrm{m}$$. The
laser-plasma energy transfer is found from the PIC simulations to be $16\%$.
According to Eq. (3), $L_{\text{d}}$ is independent of $\beta$, but Figure
5(b) shows that increasing $\beta$ to 2.4 causes $L_{\text{d}}$ to decrease to
$80\text{\,}\mathrm{m}\mathrm{m}$. The higher value of beta is also found to
increase the laser-plasma energy transfer to $19\%$, primarily due to the
earlier onset of pump depletion via SPM, as is evident by comparing the laser
spectra in Figs. 5(a) and (b). Despite the earlier onset of pump depletion for
the higher $\beta$ case, the mean energy gain for both values of $\beta$ is
approximately $\sim$1.5\text{\,}\mathrm{G}\mathrm{e}\mathrm{V}$$. This value
is in good agreement with Eq. (4), which predicts an energy gain of $\Delta
W_{e}=$1.4\text{\,}\mathrm{G}\mathrm{e}\mathrm{V}$$ for a relative wake
amplitude of $E_{\text{max}}/E_{\text{wb}}=0.5$, consistent with the amplitude
observed in PIC simulation shown in Figure 4.

Figure 5(c) shows the result for a lower density accelerator stage with
$n_{00}=$1.0\text{\times}{10}^{17}\text{\,}\mathrm{c}\mathrm{m}^{-3}$$ and
$w_{0}=$36\text{\,}\mathrm{\SIUnitSymbolMicro}\mathrm{m}$$ at $\beta=1.2$. The
evolution of the pulse train spectrum, accelerating field, and bunch energy is
shown for lengths up to $400\text{\,}\mathrm{m}\mathrm{m}$. A bunch of energy
around 2.5 GeV is obtained, which is in good agreement with the value of
$\Delta W_{e}=$2.6\text{\,}\mathrm{G}\mathrm{e}\mathrm{V}$$ given by Eq. (4)
for $E_{\text{max}}/E_{\text{wb}}=0.46$ (which is slightly lower than the
previous simulations due to the scaled laser energy density
$W_{\text{drive}}/w_{0}^{2}\lambda_{p0}$ being $8.5\%$ smaller). Figure 5(c)
shows that dephasing occurs at $z=L_{d}=$300\text{\,}\mathrm{m}\mathrm{m}$$,
which is close to the value of $296\text{\,}\mathrm{m}\mathrm{m}$ predicted
from Eq. (3). It is clear from these simulations that P-MoPAs operating at
multi-GeV energy gains should be possible.

All three simulations shown in Figure 5 exhibit the expected oscillation in
peak on-axis accelerating field due to transverse oscillations of the pulses
in the train, as discussed in the previous section and observed in Figure. 3.
In addition, Figure 5 shows that the accelerating field initially grows in
magnitude with propagation. This increase arises from: (i) SPM-induced pulse
steepening; (ii) the effect shown in Figure 4, whereby greater red-shifting of
the trailing pulses in the train causes the pulse separation to increase
towards the back of the train, partially circumventing Rosenbluth-Liu
detuning; and (iii) an increase in the ponderomotive force exerted by the
trailing pulses, which arises from their greater red-shift.

##  IX Conclusion

We have studied the operation of the accelerator stage in a P-MoPA using both
the paraxial wave equation in the weakly relativistic limit and 2D PIC
simulations.

The ability to adjust the modulator parameter $\beta$ was found to bring
several advantages. First, unlike in the plasma beat-wave accelerator, it
allows the temporal envelope of pulses within the pulse train to be
controlled. For the optimal value of $\beta=1.43$, the wake amplitude was
found to be $72\%$ larger than that generated in a PBWA with the same total
drive energy. Second, adjusting $\beta$ can minimize transverse mode
excitation, and hence oscillation of the wakefield structure, which is likely
in turn to minimize increases in the energy spread and emittance of an
injected electron bunch.

Rosenbluth-Liu detuning was found to be much less important in a P-MoPA than
it is in a PBWA, for two reasons. First, our analysis shows that in the linear
and quasi-linear regimes Rosenbluth-Liu detuning would not occur appreciably
if the number of pulses in the P-MoPA train was less than $\sim 30$. Second,
in a P-MoPA operating in the quasi-linear regime, this detuning is at least
partially counteracted by the increase in pulse spacing towards the back of
the train caused by increased red-shifting.

We presented an analysis of transverse mode oscillations within the paraxial
approximation, which was found to be in good agreement with 2D PIC
simulations. PIC simulations for several combinations of laser and plasma
parameters showed good agreement with expressions derived for the spot size
oscillation frequency, pump depletion, dephasing length and energy gain. These
simulations demonstrate that particle energy gains of several GeV should be
possible in a P-MoPA.

We note that in this study we ignored the effects of the accelerated bunch on
the wakefield. In addition to decreasing the acceleration gradient — which has
the compensating advantage of providing a means to reduce the bunch energy
spread — beamloading effects could introduce a coupling between the pulse
train and the injected bunch dynamics. This will be a subject for future
research.

The work presented in this paper confirms that P-MoPAs driven by few-joule,
picosecond pulses, such as those provided by high-repetition-rate thin-disk
lasers, have the potential to accelerate electron bunches to multi-GeV
energies at pulse repetition rates in the kilohertz range.

###### Acknowledgements.

This work was supported by the UK Engineering and Physical Sciences Research
Council (EPSRC) (Grant No. EP/V006797/1), the UK Science and Technologies
Facilities Council (Grant No. ST/V001655/1]), InnovateUK (Grant No. 10059294),
United Kingdom Research and Innovation (UKRI) ARCHER2 Pioneer Projects
(ARCHER2 PR17125) [40] and the Ken and Veronica Tregidgo Scholarship in Atomic
and Laser Physics. This publication arises from research funded by the John
Fell Oxford University Press Research Fund. This research used the open-source
particle-in-cell code WarpX [17] https://github.com/ECP-WarpX/WarpX, primarily
funded by the US DOE Exascale Computing Project. Primary WarpX contributors
are with LBNL, LLNL, CEA-LIDYL, SLAC, DESY, CERN, and Modern Electron. We
acknowledge all WarpX contributors. This research was funded in whole, or in
part, by EPSRC and STFC, which are Plan S funders. For the purpose of Open
Access, the author has applied a CC BY public copyright licence to any Author
Accepted Manuscript version arising from this submission. The input decks used
for the PIC simulations presented in this paper are available at
https://doi.org/10.5281/zenodo.10061009

## References

  * Tajima and Dawson [1979] T. Tajima and J. M. Dawson, Laser electron accelerator, Phys. Rev. Lett. 43, 267 (1979).
  * Rosenbluth and Liu [1972] M. N. Rosenbluth and C. S. Liu, Excitation of plasma waves by two laser beams, Phys. Rev. Lett. 29, 701 (1972).
  * Esarey  _et al._ [1988] E. Esarey, A. Ting, and P. Sprangle, Relativistic focusing and beat wave phase velocity control in the plasma beat wave accelerator, Applied Physics Letters 53, 1266 (1988), https://pubs.aip.org/aip/apl/article-pdf/53/14/1266/7767409/1266_1_online.pdf .
  * Gibbon [1990] P. Gibbon, The self‐trapping of light waves by beat‐wave excitation., Physics of Fluids B: Plasma Physics 2, 2196 (1990), https://doi.org/10.1063/1.859401 .
  * Hooker  _et al._ [2014] S. M. Hooker, R. Bartolini, S. P. D. Mangles, A. Tünnermann, L. Corner, J. Limpert, A. Seryi, and R. Walczak, Multi-pulse laser wakefield acceleration: a new route to efficient, high-repetition-rate plasma accelerators and high flux radiation sources, Journal of Physics B: Atomic, Molecular and Optical Physics 47, 234003 (2014).
  * Strickland and Mourou [1985] D. Strickland and G. Mourou, Compression of amplified chirped optical pulses, Optics Communications 56, 219 (1985).
  * Dawson  _et al._ [2012] J. W. Dawson, J. K. Crane, M. J. Messerly, M. A. Prantil, P. H. Pax, A. K. Sridharan, G. S. Allen, D. R. Drachenberg, H. H. Phan, J. E. Heebner, C. A. Ebbers, R. J. Beach, E. P. Hartouni, C. W. Siders, T. M. Spinka, C. P. J. Barty, A. J. Bayramian, L. C. Haefner, F. Albert, W. H. Lowdermilk, A. M. Rubenchik, and R. E. Bonanno, High average power lasers for future particle accelerators, AIP Conference Proceedings 1507, 147 (2012), https://aip.scitation.org/doi/pdf/10.1063/1.4773687 .
  * Hidding  _et al._ [2019] B. Hidding, S. Hooker, S. Jamison, B. Muratori, C. Murphy, Z. Najmudin, R. Pattathil, G. Sarri, M. Streeter, C. Welsch, M. Wing, and G. Xia, Plasma wakefield accelerator research 2019–2040: A community-driven uk roadmap compiled by the plasma wakefield accelerator steering committee (pwasc), PWASC (2019).
  * Jakobsson  _et al._ [2021] O. Jakobsson, S. M. Hooker, and R. Walczak, Gev-scale accelerators driven by plasma-modulated pulses from kilohertz lasers, Phys. Rev. Lett. 127, 184801 (2021).
  * Herkommer  _et al._ [2020] C. Herkommer, P. Krötz, R. Jung, S. Klingebiel, C. Wandt, R. Bessing, P. Walch, T. Produit, K. Michel, D. Bauer, R. Kienberger, and T. Metzger, Ultrafast thin-disk multipass amplifier with 720 mj operating at kilohertz repetition rate for applications in atmospheric research, Opt. Express 28, 30164 (2020).
  * Nagel  _et al._ [2021] S. Nagel, B. Metzger, D. Bauer, J. Dominik, T. Gottwald, V. Kuhn, A. Killi, T. Dekorsy, and S.-S. Schad, Thin-disk laser system operating above 10 kw at near fundamental mode beam quality, Opt. Lett. 46, 965 (2021).
  * Produit, Thomas  _et al._ [2021] Produit, Thomas, Walch, Pierre, Herkommer, Clemens, Mostajabi, Amirhossein, Moret, Michel, Andral, Ugo, Sunjerga, Antonio, Azadifar, Mohammad, André, Yves-Bernard, Mahieu, Benoît, Haas, Walter, Esmiller, Bruno, Fournier, Gilles, Krötz, Peter, Metzger, Thomas, Michel, Knut, Mysyrowicz, André, Rubinstein, Marcos, Rachidi, Farhad, Kasparian, Jérôme, Wolf, Jean-Pierre, and Houard, Aurélien, The laser lightning rod project, Eur. Phys. J. Appl. Phys. 93, 10504 (2021).
  * Paschotta  _et al._ [2001] R. Paschotta, J. Aus der Au, G. Spühler, S. Erhard, A. Giesen, and U. Keller, Passive mode locking of thin-disk lasers: effects of spatial hole burning, Appl. Phys. B 72, 267 (2001).
  * Südmeyer  _et al._ [2009] T. Südmeyer, C. Kränkel, C. Baer, O. Heckl, C. Saraceno, M. Golling, R. Peters, K. Petermann, G. Huber, and U. Keller, High-power ultrafast thin disk laser oscillators and their potential for sub-100-femtosecond pulse generation, Appl. Phys. B 97, 281 (2009).
  * Baer  _et al._ [2010] C. R. E. Baer, C. Kränkel, C. J. Saraceno, O. H. Heckl, M. Golling, R. Peters, K. Petermann, T. Südmeyer, G. Huber, and U. Keller, Femtosecond thin-disk laser with 141 w of average power, Opt. Lett. 35, 2302 (2010).
  * van de Wetering  _et al._ [2023] J. J. van de Wetering, S. M. Hooker, and R. Walczak, Stability of the modulator in a plasma-modulated plasma accelerator, Phys. Rev. E 108, 015204 (2023).
  * Fedeli  _et al._ [2022] L. Fedeli, A. Huebl, F. Boillod-Cerneux, T. Clark, K. Gott, C. Hillairet, S. Jaure, A. Leblanc, R. Lehe, A. Myers, C. Piechurski, M. Sato, N. Zaim, W. Zhang, J. Vay, and H. Vincenti, Pushing the frontier in the design of laser-based electron accelerators with groundbreaking mesh-refined particle-in-cell simulations on exascale-class supercomputers, in  _2022 SC22: International Conference for High Performance Computing, Networking, Storage and Analysis (SC) (SC)_ (IEEE Computer Society, Los Alamitos, CA, USA, 2022) pp. 25–36.
  * Esarey  _et al._ [2009] E. Esarey, C. B. Schroeder, and W. P. Leemans, Physics of laser-driven plasma-based electron accelerators, Rev. Mod. Phys. 81, 1229 (2009).
  * Kurtz and Streifer [1969] C. Kurtz and W. Streifer, Guided waves in inhomogeneous focusing media part i: Formulation, solution for quadratic inhomogeneity, IEEE Transactions on Microwave Theory and Techniques 17, 11 (1969).
  * Firth [1977] W. Firth, Propagation of laser beams through inhomogeneous media, Optics Communications 22, 226 (1977).
  * Esarey and Leemans [1999] E. Esarey and W. P. Leemans, Nonparaxial propagation of ultrashort laser pulses in plasma channels, Phys. Rev. E 59, 1082 (1999).
  * Sadler  _et al._ [2020] J. D. Sadler, C. Arran, H. Li, and K. A. Flippo, Overcoming the dephasing limit in multiple-pulse laser wakefield acceleration, Phys. Rev. Accel. Beams 23, 021303 (2020).
  * Schroeder  _et al._ [2010] C. B. Schroeder, E. Esarey, C. G. R. Geddes, C. Benedetti, and W. P. Leemans, Physics considerations for laser-plasma linear colliders, Phys. Rev. ST Accel. Beams 13, 101301 (2010).
  * Panofsky and Wenzel [2004] W. K. H. Panofsky and W. A. Wenzel, Some Considerations Concerning the Transverse Deflection of Charged Particles in Radio‐Frequency Fields, Review of Scientific Instruments 27, 967 (2004), https://pubs.aip.org/aip/rsi/article-pdf/27/11/967/8340739/967_1_online.pdf .
  * Gorbunov and Kirsanov [1987] L. M. Gorbunov and V. I. Kirsanov, Excitation of plasma waves by an electromagnetic wave packet, Sov. Phys. JETP 66, 290 (1987).
  * Andreev  _et al._ [1998] N. Andreev, E. Chizhonkov, A. Frolov, and L. Gorbunov, On laser wakefield acceleration in plasma channels, Nuclear Instruments and Methods in Physics Research Section A: Accelerators, Spectrometers, Detectors and Associated Equipment 410, 469 (1998).
  * von Boetticher  _et al._ [2023] A. von Boetticher, R. Walczak, and S. M. Hooker, Modulational instability in large-amplitude linear laser wakefields, Phys. Rev. E 107, L023201 (2023).
  * Picksley  _et al._ [2020] A. Picksley, A. Alejo, R. J. Shalloo, C. Arran, A. von Boetticher, L. Corner, J. A. Holloway, J. Jonnerby, O. Jakobsson, C. Thornton, R. Walczak, and S. M. Hooker, Meter-scale conditioned hydrodynamic optical-field-ionized plasma channels, Phys. Rev. E 102, 053201 (2020).
  * Feder  _et al._ [2020] L. Feder, B. Miao, J. E. Shrock, A. Goffin, and H. M. Milchberg, Self-waveguiding of relativistic laser pulses in neutral gas channels, Physical Review Research 2, 043173 (2020).
  * Miao  _et al._ [2020] B. Miao, L. Feder, J. E. Shrock, A. Goffin, and H. M. Milchberg, Optical Guiding in Meter-Scale Plasma Waveguides, Physical Review Letters 125, 074801 (2020).
  * Watson [1966] G. N. Watson,  _A treatise on the theory of Bessel functions_ , 2nd ed. (Cambridge University Press, 1966).
  * Gradshteyn and Ryzhik [2007] I. S. Gradshteyn and I. M. Ryzhik,  _Table of integrals, series, and products_ , seventh ed. (Elsevier/Academic Press, Amsterdam, 2007) pp. xlviii+1171, translated from the Russian, Translation edited and with a preface by Alan Jeffrey and Daniel Zwillinger, With one CD-ROM (Windows, Macintosh and UNIX).
  * Esarey  _et al._ [1993] E. Esarey, P. Sprangle, J. Krall, A. Ting, and G. Joyce, Optically guided laser wake‐field acceleration*, Physics of Fluids B: Plasma Physics 5, 2690 (1993), https://doi.org/10.1063/1.860707 .
  * Andreev  _et al._ [1994] N. E. Andreev, L. M. Gorbunov, V. I. Kirsanov, A. A. Pogosova, and R. R. Ramazashvili, The theory of laser self-resonant wake field excitation, Physica Scripta 49, 101 (1994).
  * Mora and Antonsen [1997] P. Mora and T. M. Antonsen, Jr., Kinetic modeling of intense, short laser pulses propagating in tenuous plasmas, Physics of Plasmas 4, 217 (1997), https://doi.org/10.1063/1.872134 .
  * Zhu  _et al._ [2012] W. Zhu, J. P. Palastro, and T. M. Antonsen, Studies of spectral modification and limitations of the modified paraxial equation in laser wakefield simulations, Physics of Plasmas 19, 033105 (2012), https://doi.org/10.1063/1.3691837 .
  * Brinkmeyer [1979] E. Brinkmeyer, Spot size of graded-index single-mode fibers: profile-independent representation and new determination method, Appl. Opt. 18, 932 (1979).
  * Mori [1997] W. Mori, The physics of the nonlinear optics of plasmas at relativistic intensities for short-pulse lasers, IEEE Journal of Quantum Electronics 33, 1942 (1997).
  * Streeter  _et al._ [2018] M. J. V. Streeter, S. Kneip, M. S. Bloom, R. A. Bendoyro, O. Chekhlov, A. E. Dangor, A. Döpp, C. J. Hooker, J. Holloway, J. Jiang, N. C. Lopes, H. Nakamura, P. A. Norreys, C. A. J. Palmer, P. P. Rajeev, J. Schreiber, D. R. Symes, M. Wing, S. P. D. Mangles, and Z. Najmudin, Observation of laser power amplification in a self-injecting laser wakefield accelerator, Phys. Rev. Lett. 120, 254801 (2018).
  * [40] https://www.archer2.ac.uk.