High-quality positron acceleration in beam-driven plasma accelerators

Acceleration of positron beams in plasma-based accelerators is a highly challenging task. To realize a plasma-based linear collider, acceleration of a positron bunch with high-efficiency is required, while maintaining both a low emittance and a subpercent-level energy spread. Recently, a plasma-based positron acceleration scheme was proposed in which a wake suitable for the acceleration and transport of positrons is produced in a plasma column by means of an electron drive beam [Diederichs et al., Phys. Rev. Accel. Beams 22, 081301 (2019)PRABCJ2469-988810.1103/PhysRevAccelBeams.22.081301]. In this article, we present a study of beam loading for a positron beam in this type of wake. We demonstrate via particle-in-cell simulations that acceleration of high-quality positron beams is possible, and we discuss a possible path to achieve collider-relevant parameters.


I. INTRODUCTION
Plasma-based particle accelerators potentially enable compact linear electron-positron colliders due to their large acceleration gradients [1]. In a plasma wakefield accelerator (PWFA), an ultrarelativistic, high-charge density particle beam expels all plasma electrons from its propagation axis and an ion cavity is formed [2,3]. The cavity, also referred to as bubble or blowout, features a region with a large, longitudinally accelerating gradient and a transversely linear restoring force for relativistic electrons. Whereas high-energy gain, high-efficiency [4,5], and stably beam-loaded [6] electron acceleration has been demonstrated experimentally in PWFAs, stable and quality preserving positron acceleration remains a challenge. Identifying a positron acceleration scheme that fulfills the requirements imposed by a particle collider, namely the stable and efficient acceleration of high-charge positron bunches, while maintaining both a low emittance and a low energy spread, has been an outstanding challenge, and previously proposed positron acceleration concepts were not able to meet all the necessary requirements. For instance, utilizing hollow core electron drive beams showed only a per-mille-level driver-to-witness energy conversion efficiency [7]. PWFAs driven by a positron beam have been investigated in Ref. [8]. While this scheme demonstrated high-efficiency acceleration of the positron witness beam, the nonlinear nature of the transverse focusing fields, and their variation as the drive beam evolves renders the preservation of the witness beam emittance challenging. Hollow core plasma channels have been proposed as potential plasma target candidates for positron acceleration [9,10]. However, owing to the lack of any focusing field for the beam in a hollow channel, this scheme suffers from severe beam breakup instability [9,11].
In a recent article, a novel method for positron acceleration was proposed that uses an electron beam as driver and a plasma column as the acceleration medium [12]. For a plasma column with a column radius smaller than the blowout radius, the transverse wakefields are altered, resulting in an elongation of the background plasma trajectories returning toward the axis. This creates a long, high-density electron filament, leading to the formation of a wake phase region which is suitable for acceleration and transport of positron beams. Despite the nonlinear nature of the transverse wakefields, it was shown that quasimatched propagation of positron beams was possible. Due to the non-uniformity of the accelerating field created in these structures, the energy-spread was found to be at the percent-level, which is too high for application in a plasma-based linear collider. Another study has investigated beam loading of simple Gaussian beams in these plasma structures [13]. Despite achieving higher efficiency than the one reported in Ref. [12], the emittance was not preserved.
In this article, we investigate beam loading of a positron bunch in the nonlinear wake formed in a plasma column with the goal of minimizing the energy spread of the bunch, while maintaining both a low emittance and a high charge. Beam loading has been first described for linear wakes in Ref. [14]. In the nonlinear blowout regime, beam loading of electron beams was studied in Ref. [15], where an analytical expression for the longitudinal witness beam current profile that eliminates the energy spread was obtained. Owing to the different nature of the wakefield structure, this type of analytic result is not valid in the case of the nonlinear positron accelerating fields considered in this study. Here, beam loading is studied by means of a numerical algorithm that reconstructs, slice-by-slice and self-consistently, the longitudinal current profile of an optimal witness beam which flattens the accelerating fields within the bunch. We further discuss the transport of the positron witness bunch and its optimization with the goal of minimizing the energy spread and preserving the emittance, both crucial parameters for the employment of this acceleration scheme in a future plasma-based linear collider.
Lastly, we assess a possible path to achieve colliderrelevant parameters.

II. NONLINEAR WAKEFIELDS FOR POSITRON ACCELERATION
The generation of positron beam focusing and accelerating wakes using plasma columns was first described in [12]. Using an electron drive beam and a plasma column with a radius smaller than the blowout radius leads to the formation of a wide longitudinal electron filament behind the blowout bubble. This elongated region of high electron density provides accelerating and focusing fields for positron beams. This is illustrated in Figs. 1 and 2, which show two-dimensional maps of the accelerating field E z =E 0 and focusing field ðE x − B y Þ=E 0 , respectively. The fields are normalized to the cold, nonrelativistic wave-breaking limit E 0 ¼ ω p mc=e, where c denotes the speed of light in vacuum, ω p ¼ ð4πn 0 e 2 =mÞ 1=2 the plasma frequency, n 0 the background plasma density, and e and m the electron charge and mass, respectively. In this example, we consider a plasma column with a radius k p R p ¼ 2.5 and a Gaussian electron drive beam with sizes k p σ ðdÞ and peak current I kA is the Alfvén current. The modeling was performed using the quasistatic particle-in-cell (PIC) code HiPACE [16]. To reduce the high computational cost of the modeling imposed by the required numerical resolution, the wakefields were computed using an axisymmetric cylindrical solver based on the one implemented in the quasistatic version of the code INF&RNO [17], while the particles are advanced in full 3D. Denoting by k p ¼ ω p =c the plasma wave number, the dimensions of the computational domain are 12 × 12 × 20k −3 p in the coordinates x × y × ζ, where x, and y are the transverse coordinates, and ζ ¼ z − ct is the longitudinal co-moving coordinate, with z and t being the longitudinal coordinate and the time, respectively. The resolution is 0.0056 × 0.0056 × 0.0075k −3 p . The background electron plasma was modeled with 25 constant weight particles per cell. The drive beam was sampled with 10 6 constant-weight particles.
The positron focusing and accelerating phase is located between −14 ≲ k p ζ ≲ −9. The accelerating field has its peak at k p ζ ≈ −11.5. Unlike in the blowout regime case, E z has a transverse dependence. The inset of Fig. 1 shows E z =E 0 along the transverse coordinate x at three different longitudinal locations denoted by the dashed (k p ζ ¼ −12.5), solid (k p ζ ¼ −11.5), and dotted (k p ζ ¼ −10.5) lines in Fig. 1.  In all three locations, E z ðxÞ has an on-axis maximum and decays for increasing distances from the propagation axis. Notably, the transverse gradient of the accelerating field is smaller further behind the driver. The nonuniformity of E z will lead to a ζ-dependent uncorrelated slice energy spread since particles that remain closer to the axis will experience a larger accelerating gradient compared to the ones further off axis. This effect will be investigated more thoroughly in Sec. III C.
The transverse behavior of the focusing field, ðE x − B y Þ=E 0 , is depicted in the inset of Fig. 2, where we show transverse lineouts of the focusing wakefields for the same three longitudinal locations used in Fig. 1. We see that the transverse wakefild decays almost linearly for increasing distances from the propagation axis. The field decrease is smaller further behind the driver. As shown in Ref. [12], the field becomes almost a step-function when sufficiently loaded by a positron bunch.

III. SELF-CONSISTENT BEAM LOADING TO MINIMIZE THE ENERGY-SPREAD
In many beam-driven plasma wakefield accelerator applications, both the driver and the witness beams are usually highly relativistic and evolve on a much longer time scale than the background plasma. In this case the quasistatic approximation [18], which allows treatment of the plasma and the relativistic beams in a separate manner, can be used. In the quasistatic approximation, the wakefields generated by a given beam are determined by initializing a slice of unperturbed plasma ahead of the beam and then follow its evolution as the slice is pushed through the beam from head to tail along the negative ζ direction (here ζ can be interpreted as a fast "time" that parametrizes plasmarelated quantities), while the beam is assumed to be frozen. This implies that to calculate the fields at some longitudinal position ζ, only the information upstream of this point is required.
We used this feature of the quasistatic solution to design an algorithm that recursively constructs, slice-by-slice and starting from the head, the optimal current profile of a witness bunch such that the accelerating field along the bunch is constant and equal to a set value. This leads to a reduced energy spread of the accelerated particles. The algorithm is described in detail in the Appendix. We considered a (radially symmetric) bunch initially described as where g k ðζÞ and g ⊥ ðζ; rÞ denote the longitudinal and transverse density profiles, respectively. We require that, for any ζ, where ζ head is the location of the bunch head, so that the bunch current density profile only depends on g k ðζÞ. For simplicity, we first consider bunches with transverse profiles that are radially Gaussian and longitudinally uniform, i.e., g ⊥ ðζ; rÞ ¼ exp½−r 2 =ð2σ 2 r Þ, where σ r is the (longitudinally constant) rms bunch size. At every longitudinal location ζ (bunch slice), the algorithm performs an iterative search for the optimal bunch current, determined via g k ðζÞ, that flattens the accelerating field in that particular slice. The procedure is repeated recursively for all the slices going from the head to the tail of the bunch. Note that, besides a constant accelerating field along the bunch, other field configurations yielding an energy chirp during acceleration are possible. In order for a solution to be found, the positron bunch has to be located in a phase of the wake where ∂ ζ E z < 0. To take into account the fact that, in general, E z varies in the transverse plane across the beam, the figure of merit considered by the algorithm is a transversally weighted accelerating field hE z i, defined as In case of a transversally uniform accelerating field, e.g., as in the blowout regime, the averaged accelerating field simply reduces to the on-axis accelerating field.

A. Optimization of the witness bunch position
The choice of the location of the witness bunch head, ζ head , sets the amplitude of the accelerating gradient and determines the shape of bunch current profile. In the following, we study the effect of different witness head positions for a bunch in the wake described in Sec. II. To fulfill the requirement that the bunch head has to be located in a wake phase such that ∂ ζ E z < 0, and to achieve a reasonable acceleration gradient, we chose −11.5 ≲ k p ζ head ≲ −10. Also, we consider a witness bunch an emittance such that k p ϵ x ¼ 0.05, and a bunch size k p σ r ¼ 0.0163. Numerical results for the current profiles and their corresponding loaded averaged accelerating fields, hE z i, for four values of the witness bunch head position are depicted in Fig. 3. Interestingly, placing the bunch head in a more forward position in the wake, corresponding to a lower accelerating gradient, does not necessarily increase the charge of the witness bunch. This can be seen in Table I, where we show the witness charge, Q w , as a function of the bunch head position.Values of the charge have been computed assuming a background density of n 0 ¼ 5 × 10 17 cm −3 . For this density the charge of the drive beam is Q d ¼ 1.5 nC.
The driver-to-beam efficiency, η, can be calculated from the charge of the witness beam, its energy gain rate, E þ w , the charge of the drive beam, and its energy loss rate, For the chosen density the driver energy loss rate is Values of the energy gain for the witness bunch and the efficiency as a function of the witness head position are given in Table I.
The results show that the efficiency peaks around the position −10.8 ≲ k p ζ head ≲ −10.6. As shown in Sec. II, the accelerating field is transversely flatter for more negative head positions, therefore the case k p ζ head ¼ −10.8 is preferable since the choice of this witness position will result in a smaller energy-spread, while maintaining close to maximum efficiency. We recall that for this witness position the charge of the bunch is 52 pC and the efficiency η ≈ 3%. This is less than what was achieved with a simple Gaussian density profile in [12], which featured a witness bunch charge of Q w ¼ 84 pC and an efficiency of η ≈ 4.8%.
However, energy-spread minimization was not taken into consideration in that study, which lead to an energy-spread on the few-percent-level.
Choosing bunch head positions that are closer to the driver, i.e., k p ζ head ≥ −10.2, yields complex (e.g., multipeaked) bunch current profiles. In this case, the positron beam significantly alters the background plasma electron trajectories, resulting in the formation of a second on-axis electron density peak behind the blowout region. This, in principle, allows for the loading of a second positron beam or an increase of the length of the first. This can be seen in Fig. 3. In fact, for k p ζ head ¼ −10.4 (dotted line) we see that hE z i has a local maximum behind the bunch which is higher than the value within the bunch, allowing for further beam loading. We did not investigate further such forward starting positions because we consider the resulting complex bunch structures difficult to realize experimentally.

B. Minimizing the correlated energy-spread
Using the weighted accelerating field hE z i from Eq. (2) as the figure of merit in the proposed algorithm yields a bunch current profile that eliminates the correlated energyspread only under the assumption that the bunch size does not change during acceleration. However, this assumption is generally not true. First, if the spot size is not matched to the focusing field at some position along the bunch due to, e.g., the slice-dependent nature of the transverse wakefields, it will evolve until it is matched. Second, due to the acceleration, the matched spot size adiabatically decreases with increased particle energy. Both effects must be taken into account in order to eliminate the correlated energy spread entirely. Eliminating the mismatch requires performing a slice-by-slice matching of the beam, i.e., introducing a slice-dependent bunch size, σ r ðζÞ. Note that this also leads to a ζ-dependence of g ⊥ . In our algorithm, calculation of the self-consistent, slice-dependent bunch size can be done numerically while the optimal bunch is generated. As a desirable side effect, the slice-by-slice matching also minimizes the emittance growth [19]. To take into account the change of σ r ðζÞ due to acceleration, the averaged spot size over the acceleration distance should be used in calculating hE z i in Eq. (2). We recall that for the here considered steplike wakes with a field strength of α the matching condition for a given emittance ϵ x is σ 3 r;matched ≃ 1.72ϵ 2 x =ðk p αγÞ and so the matched spot size is expected to scale with the energy as σ r;matched ∝ ffiffiffiffiffiffiffi 1=γ 3 p , where γ is the bunch relativistic factor [12].
The averaged bunch size over the acceleration distance can then be estimated as  where γ init and γ final refer to the initial and final bunch energy, respectively. Note that to calculateσ r ðζÞ, the final beam energy γ final is required. We also notice that the inclusion of the slice-by-slice matching and of the energyaveraged bunch size when computing the optimal bunch profiles do not significantly alter the current profiles and charges discussed in Sec. III A. Changes to the optimal beam loading algorithm including slice-by-slice matching and averaged spot size are described in the Appendix. The efficacy of the slice-by-slice matching and inclusion of the average spot size is demonstrated in Fig. 4, where we show the mean energy of each slice for a positron witness bunch that accelerates from 1 GeV to ≈5.5 GeV in a distance of 15 cm. The blue line refers to algorithm flattening hE z i with a longitudinal uniform σ r . The red line and the green line refer to additionally applying sliceby-slice matching and averaging of the bunch spot size over the acceleration distance, respectively. In this example, the location of the bunch head was k p ζ head ¼ −10.8 and the bunch had an initial emittance such that k p ϵ x ¼ 0.05. All the other parameters were as before (see Sec. II). In order to mitigate the computational cost, these results were obtained with a frozen field approximation (i.e., the particles of the witness bunch are pushed in a nonevolving wakefield). This approach has shown both reasonable agreement of the energy spread and the emittance of the witness bunch with full quasistatic PIC simulations. The agreement is facilitated by the slice-by-slice matching, which mitigates the witness beam evolution. We see that the bunch obtained without slice-by-slice matching and using the initial longitudinal uniform σ r (blue line in Fig. 4) shows a range of mean energy variation of ΔE ≈ 100 MeV. Using the sliceby-slice matching (red line) reduces the amplitude of the variation to ΔE ≈ 60 MeV. Finally, by using the energyaveraged bunch size σ r ðζÞ in the calculation of the optimal bunch (green line) the correlated energy spread is essentially removed (ΔE ≈ 3 MeV).

C. Minimizing the uncorrelated energy-spread
Whereas the correlated energy-spread can be completely eliminated, the uncorrelated energy-spread can only be reduced, as it arises from the transverse nonuniformity of E z . We did not identify a strategy to reduce the transverse gradient of E z by loading the wake with a positron bunch. However, we explored two possible solutions to minimize the impact of such gradient and reduce the uncorrelated energy spread. First, one can position the witness bunch in a region of the wake where E z is transversally as flat as possible, and second, one can use a transversally smaller witness beam.
As described in Sec. II, E z flattens transversally further behind the driver (i.e., for more negative ζ). Therefore, it is favorable to choose the starting position of the bunch furthest behind the driver, which still has a reasonable efficiency. According to this criterion and the results from Sec. III A, the optimal starting position is k p ζ head ¼ −10.8.
In the following, we study the dependence of the uncorrelated energy-spread on the witness bunch emittance. A matched bunch with a smaller emittance will have a smaller transverse extent and will sample a smaller domain of E z and, hence, it will acquire a smaller uncorrelated energy spread. For a flat beam, and assuming that in the vicinity of the axis the accelerating field can be modeled as E z ðxÞ ¼ E z;0 − βjxj, where β describes the transverse gradient of E z (see inset of Fig. 1), then, from geometric considerations, we expect the relative slice (uncorrelated) energy spread at saturation to scale as σ γ =γ ∼ βσ r =E z;0 , and so σ γ =γ → 0 in the limit of a small bunch. We note that it is not possible with our current numerical tools to model collider-relevant low-emittance witness beams, owing to the required high resolution and associated computational costs. To overcome this limitation, we use a reduced model to assess the scaling of the energy-spread for these conditions. Since the correlated energy-spread can be eliminated with the procedure discussed in the previous section, we consider a single slice of the beam in the reduced model. We chose the slice of the peak current of the positron bunch, which we have found to reasonably represent the total energy-spread of the bunch. Using the previous example with a starting position of k p ζ head ¼ −10.8, the peak of the current is located at k p ζ peak ¼ −11.45. We reuse the simulation, which included the slice-by-slice matching and the averaging over the acceleration distance. Assuming the same density of n 0 ¼ 5 × 10 17 cm −3 as for the efficiency consideration, the emittance of the beam is ϵ x ¼ 0.05k −1 p ¼ 0.38 μm. In the reduced model, we generate test particles, which we advance with a second-order-accurate particle pusher in the radial fields provided by the simulation. High-resolution simulations with the cylindrically symmetric PIC code INF&RNO indicate that the focusing field converges toward a step function [12]. Likewise, we model the focusing field in the reduced model with a piecewise constant function, ðE x − B y Þ=E 0 ¼ −αsignðxÞ, where α ¼ 0.6 for our example. We have found the model in reasonable agreement with HiPACE simulations in terms of energy-spread, emittance, and bunch size evolution. This is shown for the energyspread in Fig. 5. The black dashed line and the blue solid line describe the energy-spread at an emittance of ϵ x ¼ 0.38 μm obtained from the HiPACE simulation and the reduced model, respectively. The energy-spread of the peak-current slice of the beam is ≈0.65% for both the simulation (dashed line) and the reduced model (blue line). The final energy-spread and the emittance growth of the whole bunch in the PIC simulation are ≈0.7% and ≈2% (both not shown in Fig. 5), respectively. Under the assumption that a smaller emittance beam with the same charge does not significantly change the wake structure, we can decrease the beam emittance in the reduced model to previously numerically inaccessible values. The results are shown in Fig. 5. The results of the reduced model indicate that for emittances well below 0.1 μm, we can achieve energy spreads below 0.1%. The red and green line denote the energy spreads for initial emittances of ϵ x ¼ 0.19 μm and ϵ x ¼ 0.08 μm, respectively. Their corresponding final energy spreads are 0.3% and 0.1%. This indicates the path to possible collider-relevant parameters. However, this model does not capture the change of the wake structure due to a reduced witness bunch spot size. Eventually, when the on-axis density of the positron bunch exceeds the density of the background electrons, we expect a significant disruption of the positron accelerating wake structure. Additionally, a finite initial background plasma temperature can smooth the piecewise constant focusing field, possibly affecting the results presented here. These effects will be the topic of further research and require extensive development of simulation tools to enable detailed studies.

IV. CONCLUSION
High-quality positron acceleration with sub-percentlevel energy spread is possible in beam driven plasma wakefield accelerators. Utilizing an electron drive beam and a narrow plasma column allows for high-charge, and low-emittance positron beams. By shaping the longitudinal density profile of a transversally Gaussian witness beam, the energy-spread can be controlled and kept at the subpercent-level. Thereby, correlated energy spread can be completely eliminated. The uncorrelated energy spread scales with the transverse beam spot size. Our results indicate that using collider-relevant beam emittances might yield energy spreads as low as 0.1%. Further research will aim to strengthen this result. Additionally, the efficiency might be increased by proper shaping the drive beam [20,21], by optimizing the transverse plasma profile [12] or by using the here proposed technique to generate longitudinally chirped bunches. Extending these results to higher efficiencies will pave the path to a plasma-based collider. The algorithm used in this study calculates, by exploiting the quasistatic approximation, the longitudinal current profile of a witness bunch that maintains the average accelerating gradient over the full bunch length. The average accelerating gradient is set at the bunch head, hE z;head i. The bunch is constructed recursively by stacking infinitesimal longitudinal slices of charge, one after the other, starting from the head and going toward the tail of the bunch. For each slice, the calculation of the optimal current is done using an optimized bisection procedure.

ACKNOWLEDGMENTS
The steps of the algorithm are as follows. First, for any generic longitudinal slice i (i ¼ 0 represents the bunch head, slices are counted starting from the head), the FIG. 5. Relative slice energy spread vs acceleration distance. Advancing test particles in an approximated step function yields similar energy-spread in comparison with the HiPACE simulation.
The results indicate that emittances smaller than 0.1 μm induce an energy spread below 0.1%. algorithm computes the weighted accelerating field right behind the current slice assuming zero charge in the slice. We denote this quantity by hE z;i i. We then check if jhE z;i ij > jhE z;head ij. The absolute value is used so the algorithm works for both electron and positron witness bunches. If this condition is not fulfilled then no further beam loading is possible and the recursive procedure terminates (i.e., the bunch tail is reached). On then other hand, if the condition is satisfied, then beam loading is possible and the algorithm initializes the optimized bisection procedure to determine the current in the slice. We recall that the current is set via the g k function in Eq. (1). In order for the bisection procedure to converge, values of the current lower (g k;min ) and higher (g k;max ) than the optimal one need to be determined. Since we know that with no charge in the ith slice we have jhE z;i ij > jhE z;head ij, then we can set g k;min ¼ 0. Determining g k;max requires a trial and error procedure where, starting from, e.g., g k;max ¼ 1.2g k;i−1 , the value of the current in the slice is progressively increased in a geometric way (i.e., typically multiplying the current by a factor 10) until overloading of the wake is reached, i.e., until the condition jhE z;i ij < jhE z;head ij is satisfied. Note that every time the value of the current in the slice is changed, a solution of the quasistatic field equations for the slice is required in order to determine the current value of the weighted accelerating field behind the slice. Once g k;min and g k;max are known, the optimized bisection procedure begins. A new value of the current is computed according to g k ¼ w g g k;min þ ð1 − w g Þg k;max ; where w g ¼ ðjhE z;head ij−jhE z;min ijÞ=ðjhE z;max ij−jhE z;min ijÞ, and where hE z;min i and hE z;max i are the averaged field values behind the slice correspond to g k;max and g k;min , respectively. The bisection procedure terminates, and the algorithm advances to the next slice (i þ 1), if the averaged field computed with g k converges to hE z;head i within a predetermined tolerance, otherwise g k;min and g k;max are updated and a new optimized bisection is performed. We note that by using Eq. (A1) instead of the classical bisection procedure, i.e., g k ¼ 0.5ðg k;min þ g k;max Þ, the number of iterations required to reach convergence is significantly reduced.
Algorithm modifications for slice-by-slice matching and average bunch size Incorporating the slice-by-slice matching procedure into the algorithm requires the following modification. At each slice i, the matched spot size σ r;matched needs to be determined. This is done by exploiting a fixed-point method. We generate a Gaussian test particle distribution with some rms size, σ r ðiÞ. As an initial guess, the spot size from the previous slice, σ r ði − 1Þ, is used. Then, the test particles are evolved in time without acceleration in the focusing field given by ðE x − B y Þði − 1Þ by using a second order accurate particle pusher until the second order spatial moment of the distribution has saturated. The value of the moment is used to set a new value for σ r ðiÞ, and the whole process is repeated until the sequence of values of σ r ðiÞ has converged. Note that the focusing field of the slice i − 1 is used to compute σ r ðiÞ under the assumption that the longitudinal resolution is high enough that the focusing field changes only marginally between two adjacent slices. To further take into account the spot size reduction due to acceleration of the particles, the averaged matched spot size σ r;matched can be calculated via Eq. (4). Finally, σ r;matched or σ r;matched can be used to calculate the average accelerating field hE z i by Eq. (2). It should be noted that σ r;matched is only used to calculate hE z i, the bunch is still generated with a spot size of σ r;matched .