Influence of proton bunch parameters on a proton-driven plasma wakefield acceleration experiment

We use particle-in-cell (PIC) simulations to study the effects of variations of the incoming 400 GeV proton bunch parameters on the amplitude and phase of the wakefields resulting from a seeded self-modulation (SSM) process. We find that these effects are largest during the growth of the SSM, i.e. over the first five to six meters of plasma with an electron density of $7 \times10^{14}$ cm$^{-3}$. However, for variations of any single parameter by $\pm$5%, effects after the SSM saturation point are small. In particular, the phase variations correspond to much less than a quarter wakefield period, making deterministic injection of electrons (or positrons) into the accelerating and focusing phase of the wakefields in principle possible. We use the wakefields from the simulations and a simple test electron model to estimate the same effects on the maximum final energies of electrons injected along the plasma, which are found to be below the initial variations of $\pm$5%. This analysis includes the dephasing of the electrons with respect to the wakefields that is expected during the growth of the SSM. Based on a PIC simulation, we also determine the injection position along the bunch and along the plasma leading to the largest energy gain. For the parameters taken here (ratio of peak beam density to plasma density $n_{b0}/n_0 \approx 0.003$), we find that the optimum position along the proton bunch is at $\xi \approx -1.5 \; \sigma_{zb}$, and that the optimal range for injection along the plasma (for a highest final energy of $\sim$1.6 GeV after 10 m) is 5-6 m.

We use particle-in-cell (PIC) simulations to study the effects of variations of the incoming 400 GeV proton bunch parameters on the amplitude and phase of the wakefields resulting from a seeded selfmodulation (SSM) process. We find that these effects are largest during the growth of the SSM, i.e. over the first five to six meters of plasma with an electron density of 7×10 14 cm −3 . However, for variations of any single parameter by ±5%, effects after the SSM saturation point are small. In particular, the phase variations correspond to much less than a quarter wakefield period, making deterministic injection of electrons (or positrons) into the accelerating and focusing phase of the wakefields in principle possible. We use the wakefields from the simulations and a simple test electron model to determine the injection position along the bunch and along the plasma leading to the largest energy gain. This analysis includes the dephasing of the electrons with respect to the wakefields that is expected during the growth of the SSM. We find that the optimum position along the proton bunch is at ξ ≈ −1.5 σ zb , and that the optimal range for injection along the plasma (for a highest final energy of ∼1.6 GeV after 10 m) is 5-6 m. The latter result is obtained from a PIC simulation that tests different injection points and is also used to validate the model mentioned above.

I. INTRODUCTION
The AWAKE experiment intends to demonstrate the concept of proton-driven plasma wakefield acceleration using 400 GeV proton bunches supplied by the Super Proton Synchrotron (SPS) at CERN to accelerate externally injected electrons [1]. The concept underlying AWAKE is one of several that have been proposed for plasma-based acceleration, which could pave the way towards higher collision energies than what conventional accelerator technology can provide. An estimate for the maximum acceleration gradient supported by plasma is given by the non-relativistic wavebreaking field where c is the speed of light, m e is the electron mass, e is the elementary charge, n 0 is the plasma electron density, and ω pe = e 2 n 0 /ε 0 m e is the electron plasma frequency and ε 0 is the vacuum permittivity. The plasma density used in AWAKE, for example, of the order of 10 14 cm −3 , yields E 0 ≈ 1 GV/m, which is approximately ten times larger than what is feasible with RF cavities at the moment [2]. For higher plasma densities (10 18 cm −3 ), however, acceleration gradients of the order of 100 GV/m could be reached. Plasma-based acceleration can be accomplished using either a laser pulse or a particle bunch as a driver. AWAKE is an instance of the latter case, which is also known as plasma wakefield acceleration [3] (PWFA). As a particle bunch propagates in plasma, the fields caused by its space charge disturb the light plasma electrons, while the more massive plasma ions can be assumed to remain immobile (at the 1/ω pe time scale). The displaced plasma electrons in the wake of the particle driver oscillate at the plasma frequency ω pe /2π, and this density oscillation is in turn associated with transverse and longitudinal fields, the wakefield. The wavelength of the resulting plasma wave (or wake) is thus related to ω pe and is called the plasma wavelength: λ pe = 2πc/ω pe . When the drive bunch is short, i.e. with a typical length L λ pe , the wake travels with the speed of the driver. A charged particle can then be accelerated by propagating with roughly the same speed as the plasma wake in a region of the wakefields that is (longitudinally) accelerating and (transversely) focusing. In the linear regime, where the beam density of the drive bunch n b is much smaller than the plasma density (n b n 0 ), the transverse and longitudinal components of the wakefield are harmonic and phase-shifted by a fourth of a period with respect to each other, as expressed by a unique relationship between both components known as the Panofsky-Wenzel theorem [4]. This means that each ideal region for acceleration, where the fields are both accelerating and focusing, is only λ pe /4 long.
In order to drive the wakefields effectively, the length of the driver should be of the order of λ pe . This is not the case in AWAKE, where the bunches delivered by the SPS are considerably longer (6-12 cm) than the plasma wavelengths in the adjustable density range (∼1-3 mm for (1-10) ×10 14 cm −3 ). This causes the long proton bunch to undergo the self-modulation instability (SMI) [5], whereby the bunch is progressively modulated into a train of shorter bunches, with lengths and separation distances of the order of λ pe , due to periodic transversely focusing and defocusing fields. This instability eventually saturates and the initial proton bunch is self-consistently transformed into a bunch train, a format that can resonantly excite the wakefield.
The onset of an instability can either be due to noise or to a seed, i.e. a signal of higher amplitude than the noise level. When the SMI starts from noise, both the phase of the wakefields along the bunch as well as their amplitude vary randomly from event to event and thus prevent reliable acceleration of injected particles. In principle, seeding the instability is a means to fix the final phase and amplitude of the wakefields once the process has saturated.
It has been shown both theoretically [6] and through numerical simulations [7] that the phase velocity of the wakefields is smaller than that of the drive bunch during the growth of the SMI. This prevents acceleration since electrons can easily find themselves in the defocusing phase of the wakefields and be lost. External injection must therefore occur near or after saturation. In addition, for this injection to succeed reliably, as is required for the application as a particle accelerator, the injected bunch must be deterministically placed in the accelerating and focusing phase of the wakefields, or within a range of λ pe /4. The wakefield phase at the point of injection along the proton bunch and along the plasma must therefore be reproducible to within a fraction of that range. This must be true even in the presence of natural fluctuations of the drive bunch and of the accelerating structure, in this case the plasma. It is therefore essential to study the effect of parameter variations on the wakefield characteristics. In this work, we assume that the plasma density and thus the frequency of the wakefields does not vary. This is an assumption that is addressed in experiments by carefully controlling the plasma parameters [8].
The effects of initial bunch parameters variations are studied through numerical particle-in-cell (PIC) simulations in two-dimensional, axisymmetric cylindrical coordinates, performed with the code OSIRIS [9,10]. The values of a set of proton bunch parameters are varied independently and the respective simulations compared to a baseline simulation with the nominal AWAKE parameters. We note here that the hose instability, which can possibly compete with the SSM [11], is not described in 2D axisymmetric geometry. We therefore assume in this work that the seed for the self-modulation process is large enough to prevent the growth of the hose instability [5,12].

II. SIMULATION AND PARAMETERS
In the simulations used for this work, a moving window approximately 30 cm long and 1.6 mm high moves at c with a proton bunch (moving at ∼ c) as the latter propagates through 10 m of plasma. The profile of the proton bunch is implemented with a sharp cut, which represents the plasma creation by the co-propagating laser pulse, i.e. the relativistic ionization front, that seeds the SSM process in the experiment. In these simulations the seeding of the self-modulation process is thus modeled by the sharp rising edge of the proton bunch. The profile is given by for ξ 0 ≤ ξ ≤ ξ s , where ξ is the beam co-moving coordinate defined as ξ = z − ct, n b0 is the peak bunch density, σ zb and σ rb are the RMS bunch length and width, respectively, ξ 0 is the position where the function crosses the ξ axis (end of the bunch at ξ 0 = − πσ rb 2 for ξ s = 0), and ξ s is the seed position along the bunch. The plasma fills the simulation window up to the ionization radius r p = 1.5 mm.
Values close to the nominal AWAKE specifications were implemented in the simulations: n 0 = 7 × 10 14 cm −3 , σ zb = 12.6 cm, and σ rb = 200 µm. However, we expect the conclusions presented here to be quite general. The peak density in Eq. 2 is calculated according to giving n b0 ≈ 1.89 × 10 12 cm −3 for the nominal proton bunch population N b = 1.5 × 10 11 . The following parameters were independently varied by ±5%: σ zb , σ rb , N b and r p . The RMS timing jitter of the proton bunch with respect to the ionizing laser pulse ∆t was also varied by ±15 ps. Note that ∆t is in practice a phase shift of the cosine in Eq. (2) with respect to the center of the profile ξ s , thus encompassing either more or less charge depending on whether the maximum of the cosine is moved to the right or left of ξ s . These parameters are representative of what can be expected in the AWAKE experiment.

III. PROPERTIES OF THE WAKEFIELDS
A reliable plasma accelerator necessarily requires both amplitude and phase stability of the wakefield in the face of natural drive bunch parameter fluctuations. Phase stability is especially critical since the accelerated electrons may otherwise slip into defocusing and decelerating regions of the wakefield and be lost before gaining a significant amount of energy.
Both the wakefield amplitude and the SSM growth rate depend on the bunch density. Wakefields driven by the self-modulated bunch can reach an amplitude of the order of E z = n b n0 E 0 . Therefore, at a given plasma density, variations of the wakefields with respect to bunch parameters are expected to follow dependencies similar to that of In order to measure the effects of the different bunch parameter variations on the wakefield amplitude, we compare the average absolute value of the oscillating field E z along the propagation distance z ( |E z | ) for Relative deviation of the average absolute value of Ez resulting from the parameter scans with respect to the baseline result (see Fig. 1 each parameter. The average is computed from E z values in the simulation window at radii smaller than a plasma skin depth, r ≤ c/ω pe (c/ω pe ≈ 201 µm for n 0 = 7 × 10 14 cm −3 ).
For example, the evolution of |E z | is shown in Fig. 1 for the baseline simulation and for variations in the bunch population N b . The average fields grow rapidly until around z = 3 m, signifying the growth of the SSM, after which the SSM process saturates and the overall amplitudes of the wakefields gradually decrease. The relative difference in |E z | with respect to the baseline simulation is shown in Fig. 2 for all the parameter variations. Before SSM saturation, i.e. where linear wakefield theory is still valid (before 4.5 m), the trend E z ∝ N b σ zb σ 2 rb seems to be respected: an increase of N b by +5% produces higher values for |E z | , for example, and the variations in σ zb and σ rb cause inversely proportional effects, with the σ rb parameter variations causing the largest effects. After z ≈ 5 m, however, the effects due to N b are the largest and the effects of both bunch dimensions become of the same order.
Since the timing jitter ∆t is associated with an increase or decrease in total charge driving the wakefield (N b ), we expect the amplitude of E z to be proportional to ∆t. This is indeed observable in Fig. 2 (green curves). The magnitude of the effect is smaller than for N b variations (blue curves), because ∆t ± 15 ps corresponds to a variation of only ±2.85% in the bunch population. With our choice of plasma radius (r p ), a ±5% variation seems to have no significant effect on the wakefield amplitude. It has been shown that a smaller plasma radius can enhance the wakefield's focusing force and hence the SMI's growth rate by hindering the plasma's shielding response to the charge in the drive bunch [13]. However, this effect only becomes prominent when r p approaches σ rb , which, despite the variations of ±5%, is not the case here.
In general, the effects of the parameter variations are maximum during the growth of the SSM, reaching a relative difference with respect to the baseline of approximately 26% at z ≈ 2.8 m for 0.95 σ rb . However, AWAKE aims to have the electrons injected only after the SSM process is completed [1], at z > 4-5 m, which is why the potential for these large variations to affect the final energy of the accelerated electrons is not critical. More relevantly, at z = 5 m the maximum observed effect is of 7.5% (for 0.95 N b ), while towards the end of the plasma (z = 10 m) the effects of variations of all parameters converge to at most ±2% [for (1 ± 0.05)N b ]. This shows that the wakefield amplitude in these simulations is weakly dependent on the initial proton bunch parameters after 7 m along the plasma (see Fig. 2).
We now turn our attention to the behavior of the wakefield phase, which, as mentioned before, is central to the AWAKE accelerator concept. Assuming that an electron is moving with a constant velocity β e = v e /c in a region of the wakefield that is accelerating and focusing (and thus λ pe /4 long in linear wakefield theory), when the phase velocity of the wakefields v φ is below or above v e , the electron will eventually slip out of this region and into an undesirable one (decelerating or defocusing). This happens at the latest when the electron and the wakefields dephase by λ pe /4 with respect to each other. The larger the difference between β e and the wakefield velocity β φ , the faster this is the case.
Numerical simulation results show that during the SSM growth the phase velocity of the wakefields varies along the plasma and along the bunch, eventually converging towards that of the driver after the SSM has saturated [6,7]. This is also shown in Fig. 3, where the evolution of the wakefield phase velocity is visualized by plotting the on-axis longitudinal field component E z in a waterfall plot along the plasma. Since the simulation window moves at c, a negative slope in this type of graph means that the phase velocity of the wakefields is subluminal, while a positive slope indicates that it is superluminal. The relativistic proton bunch (with a Lorentz factor γ b ≈ 480) moves at nearly the speed of light, so its velocity essentially corresponds to a vertical line in Fig. 3 (the slope ∆z ∆ξ ≈ −2 γ 2 for bunch particles). In order to characterize the evolution of the amplitude and phase of the wakefields, we use the longitudinal wakefields E z . We make this choice because, though the transverse wakefields drive the SSM, they must be evaluated at the proper radius (e. g. at the the bunch RMS transverse size for a Gaussian profile). Since both transverse radius and shape of the bunch change as the SSM evolves, the evaluation becomes ambiguous. In contrast, the longitudinal wakefield E z is well defined and maximum on the beam axis and is the field component that leads to particle acceleration, the point of interest of this paper. In addition, the transverse and longitudinal wakefields share a fixed phase relationship due to the Panofsky-Wenzel theorem [4], which means that the phase behavior can be measured through either component.
To illustrate this, we produce a similar plot to the one in Fig. 3, but for the product of the transverse and longitudinal force components W r and W z , which, in cylindrical coordinates, are defined as W r = q (E r − c B θ ) and W z = q E z , respectively. Here, q is the charge of the test particle, E r is the radial component of the electric field and B θ is the azimuthal component of the magnetic field. This is shown in Fig. 4, where we only consider the accelerating regions, i.e. where W z > 0, and only at r = 0.7 k −1 pe , since the transverse components of the wakefield are zero on the axis. Figures 3 and 4 show that, while the wakefields are growing (z < 5 m), they are slower than the drive beam velocity (negative slope). In the region around 1.5 σ zb or 18.9 cm behind the seed [Figs. 3(a) and 4(a)] the phase velocity of the wakefields becomes essentially equal to the driver velocity after z = 5 m (vertical slope), which at first glance would make this a suitable position for the external injection of electrons. Further behind the seed [around −2.5 σ zb , Fig. 3(b)] the phase velocity remains superluminal for z > 5 m, while earlier (for example around ξ ≈ −σ zb ) it stays subluminal (not shown). This poses an advantage experimentally in that the injection position along the bunch can be scanned so as to find the optimal phase for maximum electron acceleration including possible dephasing of the electrons with the wakefields as they are accelerated.
In order to study quantitatively the effects of the parameter variations on the phase of the wakefields, we fit the function f (ξ) = A sin [k pe (ξ − ξ s ) + φ] (expected for linear wakefields) to 2.5-λ pe -long segments (starting at ξ s ) of the stacked line-outs discussed in Fig. 3, where A and φ are the fitting parameters. The value of φ is always relative to the seed position ξ s .
The result of this analysis is shown in Fig. 5 for three different positions along the bunch. Note that the red and black curves correspond to the cases in Fig. 3. Figure 5 again shows that injecting at ξ ≈ −σ zb corresponds to injecting in wakefields that are slower than the drive bunch. In the case of a finite plasma length, the injection position along the bunch must be chosen so that the wakefield phase velocity is large enough that the accelerated particles do not dephase by more than a quarter plasma wavelength. This figure indicates that, for a long plasma, injection closer to 1.5 σ zb rather than 1.0 σ zb behind the seed may be beneficial, since a slower wakefield phase velocity leads to early dephasing.
The position ξ ≈ −1.5 σ zb was chosen for the comparison of the effects from the parameter scans. This com- parison is shown in Fig. 6, where the largest effects on the wakefield phase are again observed before the saturation of the SSM, at z = 2-3 m. Here, the largest difference is of roughly 2π/20 for 0.95 σ rb at z ≈ 2.5 m. After this point it becomes difficult to extract clear-cut numbers for the effects due to noisy data from the simulations, though both this noise and possible effects from each parameter variation are limited to within 0.4 rad (approximately λ pe /16). Moreover, the phase stops changing after z ≈ 7 m in all cases, which is also the point after which the wakefield amplitude becomes essentially independent of the proton bunch parameter variations (see Fig. 2).
This suggests that, at this plasma density, electrons injected at z ≈ 7 m or further remain in phase with the wakefields for a long distance and can therefore be accelerated to high energies in wakefields with a constant phase. A subluminal wakefield velocity at the injection point (for example at z = 5-7 m) may be convenient because the electrons are injected at low energy and must first gain some energy until their speed approaches c, at which point the phase should then stay constant.

IV. BEHAVIOR OF ACCELERATED ELECTRONS
AWAKE aims to demonstrate the acceleration of an electron bunch, and therefore it is important to study the effects of initial parameter fluctuations on the properties of these electrons and not only on the wakefields, as was done so far. The characteristics of the accelerated electron bunch are the most important experimental output, and they are non-trivially dependent on several factors besides the wakefield itself, such as the electrons' initial velocity or the injection point along the plasma. Consequently, the wakefield variations reported above are not sufficient to infer possible effects on the accelerated bunches.
A simple diagnostic was devised to determine the en- ergy gain acquired by an electron as a function of its injection point along the plasma (z inj ) and its initial position along the bunch (ξ 0 ). This algorithm is in practice a one-dimensional (1D) particle pusher: for each possible z inj along the plasma, a test particle is placed at ξ 0 along the on-axis wakefield (i.e. the data presented in Fig. 3) and propagated forwards in the wakefields. All test electrons have an initial energy corresponding to γ 0 = 39.1, or approximately 20 MeV (the maximum range of the electron injector commissioned for AWAKE [1]).
The spatial resolution of these results is limited to the resolution of the simulation box in the ξ direction (which in this case means that at most 38 evenly-spaced test electrons can be tracked for every λ pe /2), while the temporal resolution is limited to the number of simulation file dumps (in this case 300, giving a maximum resolution for z inj of 3.55 cm). In this diagnostic, the electrons are assumed to remain on the axis at all times and no transverse forces are considered. Tracking particles in axisymmetric two-dimensional space (including transverse fields) would in effect entail full-fledged PIC simulations. The longitudinal field (E z ) is maximum on the axis and decays radially. An electron performing any transverse motion about the axis is subject to weaker longitudinal forces than if it is propagating exclusively along the axis (the most effective trajectory in terms of energy gain). This approach thus provides a best case scenario for the energy gained by accelerated electrons. It nonetheless includes their dephasing with respect to the wakefields, while the simplicity of the approach means that results can be obtained quickly for many different cases, i.e. for different injection points and for different simulations, such as the parameter scans performed in this work.
The result of this diagnostic is shown in Fig. 7(a) for the baseline simulation as a scatter plot of electrons that reach the end of the plasma, with their energy (colorcoded) as a function of their injection position (ξ 0 ,z inj ).
The rest of the test electrons lose enough energy at some point along z so as to slip out of the 30-centimeter-long simulation window, and hence not reach the end of the plasma. The shadowed areas should be ignored for now.
The general features of the accelerating field [see Fig. 3(b)] are visible in the point density of Fig. 7(a). Regions with few test electrons correspond to decelerating regions. Electrons in these regions lose enough energy so as to exit the simulation box. In regions where the field is accelerating (E z < 0, for example −19.00 < ξ 0 [cm] < −18.95), all the test electrons reach the end of the plasma and the final energy decreases as electrons are injected at later z positions, as expected. Their energy also decreases because the amplitude of the wakefield decreases after z ≈ 5 m (see Fig. 3). Figure 7(a) also shows that some electrons injected in the decelerating phase of the wakefields survive energy loss and can reach large energies (scattered red dots).
In order to validate the diagnostic described above, we performed a full simulation with the baseline parameters, in which test electrons were injected at 41 equally-spaced injection points between 3.5 and 7.6 m. The electrons used in the simulation also have zero emittance and are initially uniformly distributed in space (both longitudinally and transversely). The electron data was processed so as to obtain the same type of graph as Fig. 7(a). This data is shown in Fig. 7(b) for electrons injected close to the axis (r 0 < 0.5 k −1 pe ) that reached the end of the plasma.
We would expect to observe the influence of the transverse wakefields in the final energy distribution of Fig. 7(b), which is indeed the case. The regions of electron loss in Fig. 7(b) (due to transverse forces) are much clearer than those on Fig. 7(a) (due to longitudinal dephasing). The periodic regions with the most electrons in both plots [i.e. accelerating phases in (a) and focusing phases in (b)] also appear to be shifted by around λ pe /4 with respect to each other [note the shape of the scatter plot in (b) superimposed on (a)], as would be expected from the Panofsky-Wenzel theorem [4]. Despite the lower peak values in the simulation results (which is also expected, since the 1D diagnostic represents a bestcase scenario), the overall distribution matches that of Fig. 7(a).
A further comparison of the 1D pusher with direct simulation results can be seen in Fig. 8, which shows the average along each row of both graphs in Fig. 7 as well as each row's maximum value plotted against the injection point z inj . The peak energies in the simulation (dashed blue curve) are generally lower than the 1D results (solid blue curve), and their trends do not agree before z inj = 5.5 m. Nonetheless, the trend in both curves is very similar for z inj > 5 m. The average energies, in turn, show very good agreement (red curves).
We can therefore conclude that the 1D diagnostic is an appropriate tool for a comparative analysis of the effects of the parameter variations on the final energies of electrons that are initially close to the axis.  We now wish to apply the diagnostic to the AWAKE bunch parameter scans and evaluate possible effects on the final energy of injected electrons and on the optimal injection position. In order to do this, we compared the maximum final energy attained by test electrons departing from the same wakefield period (we chose the range −18.990 ≤ ξ 0 [cm] ≤ −18.956, which is approximately λ pe /4-long and is where γ f is maximal) as well as the respective coordinates (ξ 0,max , z inj,max ) for that energy for each parameter variation. Figure 9 shows scatter plots for each one of these variables (γ f,max , ξ 0,max and z inj,max ) for σ zb , σ rb and N b , the parameters whose variation caused the largest effects.
In Fig. 9(a) we find a trend of the form γ f,max ∝ N b σ zb σ 2 rb , which is consistent with the behavior observed above for the average wakefield amplitude |E z | and with the fact that the energy gain by trailing particles is directly linked to the amplitude of the axial field component E z . The resulting maximum final energies vary at most be-tween roughly −3% and +5% according to the 1D particle pusher.
For the same parameter variations, the initial positions along the bunch that led to the maximum final energies in Fig. 9(a) did not deviate from the baseline by more than approximately λ pe /10 overall [see Fig. 9(b), noting that λ pe ≈ 0.13 cm]. We would expect any effects on ξ 0,max to be connected to the phase of the wakefields, since a wakefield shift would also translate into a displacement of the optimal position for the most energy gain at each wakefield period. In the linear and strongly-coupled regime, i.e. before saturation and for k pe |ξ| (after substituting Eq. 3) roughly between z = 3.5-5 m. This is therefore consistent with the trend displayed in Fig. 9(b) for ξ 0,max , which concerns injections around 4.3 m. If N b is increased, for example, the wakefields shift backwards (φ increases), and the optimal position ξ 0,max must shift further backwards as well, i.e. ξ 0,max decreases.
As opposed to γ f,max and ξ 0,max , the injection points for the maximum energy [shown in Fig. 9(c)] determined by the 1D diagnostic do not evince clear correlations with the proton bunch parameters. This may be due to the low resolution of the diagnostic when probing different injection points along z. Though it could be argued that the different final energies in Fig. 9(a) should be coupled with respectively proportional propagation distances, and therefore z inj , one should note that these parameter variations also influence the wakefield amplitude, and the relationship between γ f and z inj is therefore not straightforward. Nevertheless, if we use the baseline values in Figs. 9(a) and (b) to compute an average wakefield amplitude during acceleration and assume that this value does not change with the parameter variations, then the range in γ f,max would correspond to a range of around 41.5 cm in the injection point. Figure 9(c) in turn suggests that the parameter variations caused a maximum variance of about 37 cm for z inj,max .
Using the 1D particle pusher as a computationally cheap diagnostic to estimate the effects of the proton bunch fluctuations on a witness electron beam, we thus found that these fluctuations do not hinder deterministic injection along the wakefields and that the effects on the experiment's most direct output, the final energy of trailing electrons, are of the same order as the imposed fluctuations (∼ ±5%). Furthermore, the results obtained with this tool are consistent with the previous analysis of the wakefield properties (amplitude and phase).
Some further discussion of the simulation results used for the validation of the 1D diagnostic is warranted. The peak energy curve obtained from the simulation in Fig. 8 (dashed blue line) suggests that the optimum injection point lies between 5-6 m. Although this graph only represents electrons initially close to the axis (r 0 < 0.5 k −1 pe ), the optimal injection range is confirmed when the final energies of all electrons at all possible radii (up to the plasma boundary r p ) are considered, as shown in Figs. 10(a) and (b). Each data point in Fig. 10(a) consists of an average of all the simulation particles that began at a given ξ 0 over the entire plasma radius, while Fig. 10(b) shows the peak energy out of all electrons with any r 0 for each ξ 0 . Both scatter plots display the highest Lorentz factors for z inj = 5-6 m. Figure 10(b) furthermore indicates that some electrons reach high energies when injected before the saturation of the SSM (z inj < 5 m), which is also suggested by the 1D diagnostic results [note red points for z inj = 0-5 m in Fig. 7(a)]. When we decompose the data in Fig. 10(b) into electrons originating above and below a radius of 1.5 k −1 pe (approximately 0.3 mm), we find that, for injections before 5 m, the electrons far from the axis attain the highest energies [ Fig. 10(d)], while for z inj = 5-6 m it is the electrons close to the axis that gain the most energy [ Fig. 10(c)].
This difference is observable for injections that take place before the saturation of the SSM and could thus be explained by its development. In fact, the PIC simulations display different behaviors for the wakefield phase velocity along the plasma radius. This is shown in Fig. 11 for the transverse wakefield component E r − c B θ (which is responsible for the radial force W r and hence focusing and defocusing). For z < 5 m, for example, the phase velocity closer to the axis [ Fig. 11(a), at r = 1 k −1 pe ] behaves as expected during the growth of the SSM and as previously discussed for Fig. 3, while at a larger radius [ Fig. 11(b), at r = 3 k −1 pe ] the phase is approximately stable between 4 and 5.5 m. This would explain why the electrons starting before z = 5 m at smaller radii would tend to be lost (due to the rapidly changing phase and their subsequent slippage into defocusing half-periods), while the electrons further away from the axis would find a stable wakefield phase and thus gain energy over a larger distance.

V. SUMMARY
Using PIC simulations, we performed parameter scans of several initial proton bunch parameters typical of the AWAKE experiment. We varied the bunch parameters σ zb , σ rb and N b , the plasma radius r p , and the seed point timing ∆t, and studied their effect on the properties of the wakefield amplitude and phase, and on the maximum energy gain and injection position as determined by test electrons. The latter analysis was based on a simple tool to estimate the maximum energy gain by trailing electrons for different injection points along the plasma while considering their longitudinal dephasing with respect to the wakefields. The effects observed with the simplified tool were validated by a PIC simulation. We found that, as long as the injection occurs after saturation of the SSM, it is possible to place elec-trons deterministically in an accelerating region, despite bunch parameter variations.
Furthermore, different aspects towards the optimization of external injection were considered. We determined that the wakefield phase velocity is essentially the speed of light around ξ ≈ −1.5 σ zb , which makes the dephasing length between electrons and wakefields as long as possible (without a plasma density gradient). The PIC simulation designed to test different injection points and validate the simplified diagnostic mentioned above indicated that the optimum injection point, i.e. leading to the highest final energies, lies within 5-6 m into the plasma for n 0 = 7 × 10 14 cm −3 . For an injection in this range, electrons close to the axis can reach energies of the order of 1.6 GeV. Comparable final energies are also attained when injection takes place before saturation of the SSM, but by electrons far from the axis instead. In the future, we will seek further optimization towards a higher accelerated beam quality, for example by including the witness beam emittance and beam loading effects in PIC simulations of the entire injection and acceleration process. It is clear that in practice all parameters vary for each event. However, as variations may have counteracting effects, we assume that the conclusions reached through single parameter variation studies are still representative of more realistic experimental situations.