Simulation and Experimental Study of Proton Bunch Self-Modulation in Plasma with Linear Density Gradients

We present numerical simulations and experimental results of the self-modulation of a long proton bunch in a plasma with linear density gradients along the beam path. Simulation results agree with the experimental results reported in arXiv:2007.14894v2: with negative gradients, the charge of the modulated bunch is lower than with positive gradients. In addition, the bunch modulation frequency varies with gradient. Simulation results show that dephasing of the wakefields with respect to the relativistic protons along the plasma is the main cause for the loss of charge. The study of the modulation frequency reveals details about the evolution of the self-modulation process along the plasma. In particular for negative gradients, the modulation frequency across time-resolved images of the bunch indicates the position along the plasma where protons leave the wakefields. Simulations and experimental results are in excellent agreement.


I. INTRODUCTION
A plasma wakefield accelerator (PWFA) [2] uses a relativistic particle bunch as driver. The energy gain by a witness bunch (electrons, positrons, muons, etc.) is smaller than, or equal to the energy lost by the drive bunch. Two options exist to produce very high-energy (TeV) witness bunches carrying kJ of energy, i.e., with ≈ 10 10 particles: staging of multiple plasma sections with wakefields driven by ≈ 100 J energy level bunches [3] or using a single plasma with wakefields driven by a multi-kJ energy bunch, such as a proton bunch [4]. High-energy proton bunches available today have a root mean square (rms) length σ z of 6 to 12 cm (corresponding to an rms duration σ t of 200 to 400 ps), e.g., CERN SPS and LHC, and Brookhaven RHIC.
In a plasma with density n e0 , a bunch with a Gaussian longitudinal profile can effectively drive wakefields when the plasma wavelength λ pe is on the order of its rms length σ z as λ pe = π ne0re 1/2 ≈ σ z , where r e is the classical electron radius. In that case, the amplitude of the driven accelerating wakefields is on the order of the wave-breaking field [5] E W B (σ z ) = 2π mec 2 e 1 σz , where m e is the rest electron mass, e its charge, and c is the speed of light in vacuum. This scaling yields < 55 MV/m fields for these long proton bunches. To drive wakefields with an amplitude in the order of GV/m, a long proton bunch must undergo the process of self-modulation (SM) [4].
Self-modulation occurs when a long, relativistic, charged particle bunch propagates in a dense plasma, i.e., when σ z λ pe . To avoid filamentation [6] the bunch must have an rms width σ r < λ pe . The SM transforms the long bunch into a train of microbunches shorter than, and with a periodicity of approximately λ pe . The bunch train then resonantly drives wakefields with amplitudes that can reach a significant fraction of E W B (λ pe ). The plasma density can thus be adjusted so the amplitude of the wakefields can exceed the GV/m level in a plasma with n e0 > 10 14 cm −3 , even with a cm-long bunch.
During the SM growth of a long beam, the phase velocity of the wakefields is slower than the velocity of the drive bunch v b and evolves along the bunch ξ = z − ct and plasma z = ct [1,7,8]. The difference between the phase velocity and the bunch velocity is given by [7] |∆v| ≈ 1 2 ξ z 1/3 n b m e n e0 m p γ p where n b is the bunch density, m p is the proton mass, and γ p is the gamma factor of the protons. This slow velocity causes a phase shift between the microbunches and the wakefields during the SM growth and possibly after that. The effect of ∆v is cumulative along the bunch and thus more important at the back of the bunch than in its front. This means that particles of the drive bunch do not necessarily remain in the focusing or defocusing phase of the wakefields all along their propagation in the plasma, which affects the process of the microbunch train formation. Similarly, they do not necessarily remain in the accelerating or decelerating phase of the wakefields, which also affects the amplitude of the wakefields.
The use of a positive plasma density gradient was proposed to produce an effective positive additional phase velocity to compensate at least for some of the deleterious effects of phase velocity variations during SM growth [7,8]. The effect of positive gradients is to shorten λ pe along the bunch path, causing any given constant phase point to move toward the seeding position (Eq. 1), thereby counteracting the backwards phase shift. The effect of negative gradients is to lengthen λ pe , enhancing the phase shift. This is most important for the acceleration of electrons injected early along the plasma. In a constant density plasma, i.e., no gradient, these electrons could quickly dephase in the slow wakefields and gain little energy or even be defocused and lost. Plasma density gradients are also useful for maintaining particles of the drive bunch in phase with the wakefields, leading to trains with more microbunches [1], possibly leading to larger accelerating fields [9].
We showed the effects of plasma density gradients g on the SM of a 400 GeV proton bunch in an experimental paper [1]. These results clearly show that positive gradients along the plasma lead to a higher number of microbunches in the train with more charge per microbunch. Measurements also show that the bunch modulation frequency f mod , measured after 10 m of plasma, changes as a function of gradient value.
While the experimental results show many of the features expected in the presence of density gradients, limitations with access to the SM process (observation of the bunch only after propagation in plasma and vacuum), as well as the lack of plasma density perturbation diagnostics, restrict the ability to measure many characteristics of the SM process along the plasma and as a function of g values. However, in numerical simulations one has access to the beam, plasma, and wakefields characteristics all along the plasma.
Here we show simulation results obtained with the particle-in-cell code OSIRIS 4.4.4 [10] using parameters similar to those of the experiment. We show that there is a remarkable agreement between simulation and experimental results. We look into the details of the evolution of the proton bunch and wakefields along the plasma to explain the results obtained after the plasma. We relate the amount of charge in the core of the microbunch train to the dephasing of the wakefields with respect to the proton bunch, as well as to the amplitude of the wakefields. We show that the off-axis, time-resolved proton distribution carries information about f mod earlier along the plasma and that f mod for the microbunch train core evolves very differently along z for positive and negative gradients. We also show that, among the values we use, there is indeed a positive gradient value that leads to a larger amplitude of the transverse wakefields when com-pared to the constant density case. We show however that a linear density gradient does not eliminate the issue of decreasing amplitude of the wakefields along the plasma past the SM saturation point.

II. GENERAL SETUP
The SPS at CERN provides the proton bunch for the experiment. For the experimental results presented here, it has a population of (2.98 ± 0.16) × 10 11 protons, each with an energy of 400 GeV. The rms duration of the bunch is σ t = 230 ps (or equivalently a rms length σ z = 6.9 cm). As we use the same data as was used in [9], we also consider that the bunch is focused to an rms transverse size of σ r0 ≈ 200 µm at the entrance of the 10 m-long plasma and that it has a normalized emittance The SM process is seeded, in the experiment by a relativistic ionization front (RIF) [11][12][13], and in simulations by a step in the density profile of the proton bunch at the same position along the bunch as that of the RIF. Here, the two seeding methods are considered to be equivalent.

III. SIMULATION PARAMETERS
We perform simulations in 2D axisymmetric geometry, using the parameters in Table I chosen after suitable convergence tests. We also compared results with a wider simulation window (0.32 cm) for selected g values, and showed that for this study the halo effects [9] are negligible and we can use a narrow window to save computation time. The use of 2D simulations excludes 3D effects, which could have appeared in the experiment, such as the hosing instability [14]. However, 2D simulations reproduce the SM and many of the effects observed experimentally, while no severe hosing was observed. The simulation window is moving at c together with the proton bunch. The simulation output data is saved every ≈ 10 cm (i.e., every 36928 time steps). We propagate simulation particles in vacuum to the locations of the screens in the experiment, 2 m and 3.5 m downstream from the plasma end.
In the experiment, bunch parameters are measured as a function of time at a fixed position (that of the second screen). In the simulations, all quantities are measured at a fixed time, over the spatial extent of the proton bunch.
However, the evolution of the proton distribution over its duration or length is sufficiently small that the results do not have significant differences.

IV. EXPERIMENTAL PARAMETERS AND DIAGNOSTICS
The plasma is created by a 120 fs-long laser pulse, with a radius of ≈ 1 mm and an energy of ≈ 110 mJ, ionizing a Rb vapor and co-propagating within the proton bunch, 128 ps ahead of its density peak. The laser pulse ionizes the outermost electron of each Rb atom. It creates a RIF much shorter than the period of the wakefields (> 8 ps). The interaction of the bunch and the plasma starts suddenly at the RIF, which provides seed wakefields for the SM to grow from [12,13]. When seeded, the phase of the bunch modulation with respect to the RIF is reproducible, from event to event, within a small fraction of a modulation period [13].
The plasma density is obtained by measuring the Rb vapor density with white-light interferometry with an uncertainty of 0.5 % [15]. The Rb and plasma densities agree with each other [12]. We therefore quote the plasma density hereafter, even though the Rb density is measured.
The density gradient is created by controlling the temperature of the reservoirs that evaporate Rb into the source at the plasma entrance and exit [16,17]. The experimental density gradient is calculated by dividing the difference in densities at the entrance and exit of the plasma by its length. In all cases, we keep the density at the plasma entrance, n e0 = 1.81 × 10 14 cm −3 , constant and vary that at the exit.
After the plasma, the proton bunch propagates in vacuum towards screens where it is characterized. We image onto the entrance slit of a streak camera the backward, incoherent, optical transition radiation (OTR) produced when protons enter the metallic screen located 3.5 m downstream from the plasma exit [12,18]. The camera gives a time-resolved image of the transverse density distribution of the proton bunch in a ≈ 74 µm-wide slice of the bunch about its axis. We obtain experimental images from multiple short-time-scale (209 ps) images stitched together in time [19]. The fact that these stitched images form long trains of microbunches (see Fig. 1 (a,c,e)) as opposed to a blur, confirms that the SM process is in the seeded [13] and not in the instability regime [4].

V. DENSITY PROFILES OF THE MICROBUNCH TRAIN
We display in Fig. 1 time-resolved images of the charge density distribution from the experiment [18] and from simulations for selected g values, together with their onaxis time profiles. These images are part of a series that includes all measured g values. The experimental images for g = −2 and +2 %/m are displayed in [1] and together with the images displayed here show the same features. We emphasize here the similarities and differences between experimental and simulation images.
We generate the experimental profiles by summing counts of the images up to a transverse extent of |x| = σ r,screen = 0.536 mm from the axis (red lines on Fig. 1), where σ r,screen is the rms width of the incoming (no plasma) bunch at the screen position. We obtain this width from a Gaussian fit to the transverse profile of the part of the bunch ahead of the ionization front (t < 0 ps on Fig. 1 (a,c,e)). In simulations, we sum the charge of macroparticles in bins with size equal to that of a pixel from the experimental images, and divide by the volume in the ring corresponding to each bin (2D simulations). We calculate on-axis time profiles by summing the density profiles transversely up to σ r,screen , as for the experimental images.
It is clear that the resolution of the simulation images is better than that of the experimental ones. These are subject to the space (≈ 180 µm) and time resolution (≈ 3 ps) of the diagnostic. We use g = 0 %/m as reference to compare to the other two cases. Figure 1 shows that images (c) from the experiment and (d) from simulations with g = 0 %/m look quite similar. In both cases, the density of the microbunch train decreases from the t = 0 ps seeding position to t ≈ 200 ps. After that point, there is a series of microbunches with a small and relatively constant density when compared to that of the ones at the front. We also note that as a result of focusing, the first two or three microbunches have a charge density larger than that of the incoming bunch (t < 0 ps on experimental images, not shown on simulation images).
Images with negative gradient (g = −1 %/m, Fig. 1  (a,b)) also show similar characteristics between simulations and experiment. In this case, the microbunch train is shorter than with g = 0 %/m. Only a few microbunches remain, between t = 0 ps and t ≈ 100 ps. The defocused charge can be seen away from the axis, outside the red lines at t > 100 ps.
In the positive gradient case (g = +1 %/m, Fig. 1 (e,f)), the charge density of the microbunches remains approximately constant up to t ≈ σ t = 230 ps. In that region there are 27 microbunches at this plasma density, both in simulations and in the experiment. After that point, there is a sudden decrease of the charge density per microbunch on the experimental image and profile ( Fig. 1  (e)). In the simulation result ( Fig. 1 (f)), the charge density starts decreasing after t ≈ 300 ps and reaches a value close to zero at t ≈ 500 ps.
This quick analysis of Fig. 1 shows that there is a good general agreement between experimental and simulation results. It also shows that experimental results with possible 3D effects ( Fig. 1 (e)) differ somewhat from simulation results obtained in 2D ( Fig. 1 (f)), especially at later times along the bunch. Since the bunch train is generally longer with g > 0, the measurement is more sensitive to possible misalignment between the bunch and the streak camera slit. With the longer bunch train with more charge, the non-axisymmetric hose instability [14] is more likely to develop and steer the bunch centroid away from, and sideways along, the slit of the camera as shown by the slight wiggles of the microbunches around the axis (between t ≈ 100 ps and ≈ 250 ps) in Fig. 1 (e). These effects may explain the longer train observed in simulations with g = +1 %/m. Figure 1 shows that the effect of the plasma density gradient on the bunch charge density distribution can be threefold: change the length of the bunch train, i.e., the number of microbunches in the train; change the charge density and charge of individual microbunches; increase the likelihood of developing observable non-axisymmetric effects.

VI. CHARGE OF THE MODULATED PROTON BUNCH
In the experiment, we determined the charge fraction in the core of the microbunch train (charge within one rms width normalized to the charge in the same volume of an unmodulated proton bunch) from time-integrated images of the transverse bunch distribution for each g value. We obtain these images 2 m downstream from the plasma exit. The results were published in [1] and are reproduced in Fig. 2 (blue squares). The figure shows that there is a clear difference in the amount of charge measured with g > 0 and g ≤ 0. As mentioned in [1], charge variations can be attributed to changes in phase velocity and amplitude of the wakefields along the plasma due to the density gradient.
With simulation data we obtain the charge of the microbunch train that is within one rms width of a bunch propagating in vacuum, after 2 m of propagation from the plasma end. We then normalize it to the charge of an unmodulated bunch in the same volume. We finally add the fraction of charge ahead of the step in the den- sity profile (≈ 29 %), which is not present in simulations, but is in the experiment. Figure 2 shows that there are essentially two charge fraction values, one for g > 0 and one for g ≤ 0. The values are similar between simulation and experiment for g > 0, while the simulation values are ≈ 0.15 lower than the experimental ones for g ≤ 0.
To better understand this feature, we now investigate with numerical simulations the evolution of the bunch charge fraction along the plasma, displayed in Fig. 3. We calculate it as explained above, as a function of propagation distance in the plasma (z < 10 m) and until the screen (z ≥ 10 m) for all g values. We calculate the bunch transverse size at each z-position: . The charge fraction values at z = 12 m in Fig. 3 correspond to the ones on Fig. 2.
The figure shows that, in all cases, the charge fraction increases over the first 3 m due to the overall focusing effect of the transverse wakefields and of the plasma adiabatic response. After that there is significant loss of charge due to the SM process and microbunch train formation. Figure 3 shows that past z ≈ 5 m, the charge fraction continues to decrease for all g values due to the continuous evolution of the wakefields and proton distribution. Figure 4 (similar to Fig. 9 in [9]) shows that the mean defocusing wakefields [20] reach their maximum amplitude, i.e., saturation, between 3 and 5.5 m, sooner for g < 0 (red lines) and later for g > 0 (blue lines). For g < 0, the charge fraction decreases more than for g > 0, and it does so faster immediately after saturation, be- tween z ≈ 3 m and z ≈ 7 m, and slower afterwards. The charge fraction value for g = 0 %/m (black line) is between the ones for g = 0 %/m and ends with a value similar to those of g < 0 at the measurement location. We note here that Fig. 4 shows that g = +0.5 %/m yields the largest mean transverse wakefield amplitudes and it is the only one with a significantly larger value than g = 0 %/m. We investigate in simulations the evolution of the phase and amplitude of the wakefields to provide a more detailed explanation for the evolution of the charge fraction seen on Fig. 3. Since the transverse position of the maximum value of transverse wakefields changes along the plasma, we use the longitudinal wakefields to analyze the phase of the wakefields. To get a handle on the phase, we follow the position x 0 of the zero-crossing of the onaxis longitudinal wakefields and its consequences on the local bunch charge along the plasma, as opposed to the evolution at fixed ξ 0 positions. The rationale is that since wakefields and charge distributions shift in phase along the plasma, we follow the wakefields rather than single protons. For that purpose we choose three representative initial positions along the bunch, near the step in the density profile in the front of the bunch (ξ 0 ≈ 1 cm), after the peak of the bunch density and of the amplitude of the wakefields (ξ 0 ≈ 7 cm), and to the rear of the bunch (ξ 0 ≈ 14 cm).
To calculate the charge fraction evolution and mean amplitude of the defocusing wakefields at the different ξ 0 positions, we consider that, with a linear plasma density gradient, the plasma wavelength changes along the plasma as λ pe (z) = λ pe0 1 + g 100 z −1/2 , with λ pe0 = λ pe (z = 0 m). Figure 5 (a-c) shows that the position of x 0 shifts backwards for g = 0 %/m and shifts more the further along the bunch. They also show that it does not shift according to a local adjustment of the plasma wavelength as λ pe ∝ n e (z) −1/2 , because it is a combination of the velocity difference described by Eq. 1, the further evolution of the microbunch train after saturation, and the effect of the density gradient. The magnitude of the shift in x 0 is also different between g > 0 and g < 0, being larger in the latter.
Early along the bunch (ξ 0 = 1 cm, ≈ 4λ pe0 ), Fig. 5  (c) shows that the dephasing of the wakefields is small (< 2(λ pe0 /2), wakefields change from focusing to defocusing or vice versa over λ pe /2) and their amplitude is small (< 10 MV/m, Fig. 5 (f)) in all cases. Thus, most of the charge loss (Fig. 5 (i)) occurs after z ≈ 4 m, after the initial focusing phase, and is dominated by the continuous evolution of the wakefields past their saturation point. This early along the bunch, all quantities are similar among the various g values because the cumulative effect of the change in period of the wakefields caused by the change in plasma density is small after only ≈ 4λ pe0 ). The charge fraction remaining at the end of the plasma (≈ 70 %) is in general larger than further along the bunch (Fig. 5 (h,g)), because of the interplay between the low transverse wakefields protons are subject to, which allows them to be recaptured by the focusing wakefields after having been weakly radially pushed out by the defocusing ones.
The situation is quite different near the middle of the bunch (ξ 0 = 7 cm, ≈ 28λ pe0 ). There, the gradients cause a split (Fig. 5 (b,h)) between g > 0 for which dephasing remains small ( λ pe0 /2) and g ≤ 0 for which dephasing exceeds −λ pe0 /2 (Fig. 5 (b)) and reaches −8λ pe0 /2 (for g = −2 %/m). The maximum amplitude of the mean defocusing wakefields exceeds 50 MV/m in all cases and varies by a factor of almost two with g (Fig. 5  (e)). However, with g = +0.5, +1, and +1.5 %/m, despite the larger amplitude of the wakefields than earlier along the bunch (Fig. 5 (f)), the small dephasing leads to less charge loss. With these g values one observes an effect close to what is desired from the SM process that may be obtainable with a plasma density step [21]: small or no charge loss after saturation, signifying the driving of wakefields with constant amplitude over a long distance. That amplitude depends on the charge in each microbunch, but also on the position of the microbunches within the accelerating and decelerating wakefields. Thus, even with a similar amount of charge, amplitudes might differ. With g > 0, the charge fraction remains larger than 50 % at z = 10 m, with the most charge for g = +1 %/m. Cases with g ≤ 0 have very low charge fraction values at z = 10 m.
At the back of the bunch (ξ 0 = 14 cm, ≈ 56λ pe0 ), Fig. 5 (a) shows that the dephasing is similar to that observed in the middle of the bunch (Fig. 5 (b)), although cases with g < 0 have a larger dephasing. With g > 0, x 0 takes positive values, i.e., the wakefields become superluminal. The amplitudes of the wakefields (Fig. 5 (d)) are in general smaller than those in ξ 0 = 7 cm. With g = +0.5 %/m ( Fig. 5 (g)), the charge fraction remains high at z = 10 m and the phase (with respect to the bunch) after saturation remains relatively constant when compared to other g values. This is the result of the focusing of the charge in the region r > σ r to the axis, in addition to considering off-axis charge trapped between focusing and defocusing wakefields, but still within r < σ r , which leads to charge fraction values larger than 1. We note here that the charge at this ξ 0 = 14 cm position along the bunch (235 pC initially in the ξ 0 ± λ pe0 /2 interval, > 600 pC at the other two locations) weakly contributes to the total charge variations.
These results are consistent with what is observed on the bunch images of Fig. 1. Simulation results show that the charge preservation or loss depends on a combination of the effects of the phase and amplitude of the wakefields at the location of the protons along the bunch and all along the plasma. The case with g = +0.5 %/m shows quite a constant phase velocity along the bunch, and in general the largest amplitude of the wakefields and a high charge fraction remaining.

VII. PROTON BUNCH MODULATION FREQUENCY
In the experiment, we demonstrated that, after propagating in 10 m of plasma, the bunch modulation frequency is equal to the plasma frequency f pe as f mod = f pe = c 2 ne0re π 1/2 [12]. We used different values for the plasma density, covering a range of one order of magnitude. For each value, the plasma density was kept constant (no gradient). When using a plasma with a linear density gradient (g = 0 %/m), f mod lies between the plasma frequencies at the plasma entrance and exit [1]. Two different measurement methods were used and they showed good agreement: the CTR diagnostic [22] and applying a discrete Fourier transform (DFT) frequency analysis to time profiles of time-resolved images of the bunch such as those of Fig. 1.
Carefully observing bunch images on Fig. 1 (a-f), one notices that the charge distribution of the microbunches curves away from the beam axis, showing a general Cshape for g = −1 %/m ( Fig. 1 (a,b)). The distribution appears quite straight for g = 0 %/m ( Fig. 1 (c,d)). These curvatures are displayed also for g = −2 and +2 %/m in the inset of Fig. 7. Assuming a modulation that starts at the seeding position, such features could be an indication of f mod variations across the charge distribution at that location. Knowing that with density gradients the plasma frequency varies along the plasma, it may be expected that protons reaching a larger radial position at the screen have left the wakefields earlier along the plasma.
We performed a DFT analysis on time profiles from simulation and experimental images after propagating 10 m in plasma and 3.5 m in vacuum. In both cases, we used a time window starting 12 ps behind the seeding position to exclude from the analysis the first microbunch, which is usually longer than the following ones [23]. The window extends to 467 ps (≈ 2σ t ) behind the seeding position. The width of the DFT frequency bin for this time window is 2.2 GHz. By zero-padding the profiles, the width is reduced to 0.3 GHz, which is similar to that obtained from the accuracy of the plasma density measurement [15]. We take the frequency of the highest peak in the DFT power spectrum as f mod . Figure 6 shows f mod as measured in simulations (red symbols) and in the experiment (blue symbols), in both cases using a narrow transverse extent of |x| ≤ 0.36 mm (filled symbols), that includes only the microbunch train, and a wide one of |x| = [0.36, 1.8] mm (empty symbols) for the time profiles.
Simulation and experimental results are in very good agreement with each other. The value of f mod is proportional to f pe (z = 10 m) over the −0.5 ≤ g ≤ +0.5 %/m range in simulations and experiments. For g > +0.5 %/m, f mod saturates at ≈ 3 % above f pe0 . For g < −0.5 %/m, f mod ≈ f pe0 in the wide window case, whereas it follows f pe (z = 10 m) for the narrow window

case.
To further explore the behavior of f mod with respect to the distance from the axis, we perform a DFT analysis using multiple time profiles starting at the beam axis and with radial extents equivalent to the transverse resolution of the optical system in the experiment: 180 µm. For the figure, we mirror about the beam axis the frequency values from 2D simulations. Figure 7 (a), for the gradient value that shows the largest difference in Fig. 6 (g = −2 %/m), shows that f mod is equal to f pe (z = 10 m) near the axis (|x| ≤ 0.36 mm) and transitions to f pe0 for |x| ≥ 0.96 mm. On the contrary, Fig. 7 (b) shows that with g = 0 %/m the frequency remains close to f pe0 at all radii, in agreement with the single frequency expected. Figure 7 (c) for g = +2 %/m, shows a slight increase in frequency with radius, unlike Fig. 7 (a). There is very good agreement between the experimental and simulation results. The result for g = −2 %/m is consistent with the above hypothesis that protons reaching larger radii at the screen left the wakefields earlier along the plasma, carrying information about the plasma frequency at that location. In the other two cases (b,c), the variations are too small to draw strong conclusions from these results. The variation of f mod with radius also explains the curvature of the distributions visible on Fig. 1 images and insets of Fig. 7.
To confirm the origin of the proton distributions carrying certain frequencies, we back-track simulation particles along the plasma. We identify the macroparticles in each radial slice at the screen position (z = 13.5 m) and x (mm) calculate their mean radial position at each location along the plasma. Figure 8 shows that indeed, for g = −2 %/m, particles reaching the larger radii at the screen left the wakefields early along the plasma (z ≈ 3 m) and carry a f mod equal to f pe0 . The figure also shows in essence two populations: one that leaves the wakefields early and one that remains within the wakefields over most of the plasma length, both of them carrying the local plasma frequency as f mod , as shown on Fig 7 (a). The proton distributions, especially in the core, continue evolving even after the plasma, as they move towards or away from the axis due to their transverse momentum. We explore the difference in the frequency response among the various g values by following the evolution of f mod near the bunch axis (|x| < 180 µm) along the plasma and up to the screen position.
In Fig. 9, the modulation frequencies start around 124 GHz, higher than f pe (n e0 ) = 120.8 GHz, because the presence of the proton bunch density introduces an additional restoring force on plasma electrons, which increases their initial oscillation frequency: f pe (n e0 + n b ) = 123.3 GHz, and therefore also the initial f mod . For all g values, there is an initial decrease in f mod over the first meters of plasma caused by phase slippage during SM growth. Afterwards, as with the charge distribution and dephasing, there is again a different behavior between positive and negative gradient cases. The f mod for g > 0 is essentially constant along z, with relative changes between the minimum and maximum value of ≈ 1.5 %. The bunch train is relatively long and thus cannot adjust well to the local plasma frequency, and its effectiveness at driving wakefields decreases for large g values (see Fig. 4), despite the large charge in the train and per microbunch (see Fig. 2). It is interesting to note that the f mod = f pe for g = 0 %/m reported in [12] is reached only near the end of the plasma, when the incoming bunch is fully modulated. The g = 0 %/m case has again an evolution somewhat between that of the g > 0 and g < 0 cases. The dominant f mod for g < 0 that is first slowly decreasing along the plasma, exhibits sudden jumps to basically follow f pe (z) after some point between z = 5.5 m and z = 7 m.  Figure 9 shows limited information about the DFT power spectrum because it only displays the frequency with the highest amplitude. However, Fig. 10 shows that for (a) g = −2 and (b) g = +2 %/m the DFT power spectrum does not have a single peak along z, and amplitude variations in the wakefields and bunch density modulation, as well as structures in the transverse charge distributions caused by successive defocusing and focusing of protons create multiple spectral features. For g = −2 %/m ( Fig. 10 (a)), a component with a lower and decreasing frequency starts developing at z ≈ 5 m. Its amplitude dominates the one at the higher frequency at z = 5.5 m, which causes the sudden jumps in frequency seen in Fig. 9. Then, the peak frequency follows f pe (z). Around this point the bunch train becomes significantly shorter (see loss of charge on Fig. 3) and the driving of wakefields becomes dominated by microbunches early along the train. Because the relative dephasing between stretching wakefields and these early microbunches is smaller than with later ones that are totally defocused, their periodicity can better adjust to the local frequency at the expense of becoming shorter. Figure 11 (a) shows this effect and indicates that focusing wakefields can also bring charge back towards the axis, to regions that were previously depleted of charge, e.g., ξ ≈ 0.60, 0.85, 1.10 cm, with more charge later along the plasma than earlier. Figure 10 (b) for g = +2 %/m shows that in this case a second frequency appears at z ≈ 6 m with value f pe (z = 6 m), which remains until the end of the plasma, but with a lower amplitude than that of the initial frequency value. In that case, after the modulation has saturated and the charge between microbunches has been completely depleted, the defocusing wakefields, which have a different frequency than the train, expel charge from the microbunches and modify the envelope of the train (z = 10 m). Figure 11 (b) shows the resulting beating pattern. The experimental results also show two frequencies in the CTR spectrogram of some g = +2 %/m events [1]. Streak camera images do not have sufficient signal to noise ratio to display the actual complexity of the frequency spectrum of numerical simulations. These figures show the intricacy of the frequency spectrum of the modulated bunch, resulting from the complex evolution of the bunch train and wakefields.
The frequency analysis of the radial distribution after the plasma shows that detailed information can in general be retrieved from the experimental profiles and confirmed by simulation results, and even explained. This is possible because there is good agreement between experimental and simulation results, even though the SM process is an intricate one. The images of Fig. 1 show the distribution of particles that drove wakefields towards the end of the plasma (with the added transverse evolution in vacuum between the plasma exit and the screen). They therefore also have some of the characteristics of the wakefields themselves, as seen in simulations. In that sense, they carry information similar to that obtained from various plasma density perturbation diagnostics such as interferometry, shadowgraphy [24], and Fourier domain holography [25]. This correspondence could be established with these diagnostics, as is the plan for later experiments.

VIII. CONCLUSION
The detailed simulation and experimental results presented here regarding the charge in the microbunch train, its distribution and modulation frequency, as measured after the plasma as a function of linear plasma density gradients, are in excellent agreement with each other. They complement and explain our previously published results [1] and show that many details of the selfmodulation process are observed in the experiment. Density gradients change the phase velocity of the wakefields that we observe as a change in modulation frequency of the proton distributions measured along the plasma in simulations, and after the plasma in simulations and experiments. The effect of the density gradients is also to change the charge in the microbunch train and in the microbunches, as well as to change the amplitude of the wakefields. Time-resolved images of the proton bunch charge distribution contain information about the evolution of the self-modulation process along the plasma. For example, the modulation frequency of the defocused proton density distribution indicates the position along the plasma where protons left the wakefields. We will record similar data when optimizing a plasma density step [21], an extreme case of density gradient, for maximum plasma density modulation or maximum energy gain by externally injected electrons. This data will bring additional information about the effect of the density step on the proton bunch modulation.