Plasma-based acceleration (PBA) uses an intense laser pulse [1] or a charged particle beam [2] to excite a plasma wake, which can be utilized to accelerate electrons and positrons with high acceleration gradients [3][4][5][6][7][8]. The acceleration gradients facilities that are built for conducting PWFA research, such as Facilities for Accelerator Science and Experimental Test (FACET) II [22], Advanced Proton Driven Plasma Wakefield Acceleration Experiment (AWAKE) [23], Future Oriented Wakefield Accelerator Research and Development at FLASH (FLASHForward) [24] and EuPRAXIA [25]. In PWFA, when the highly relativistic drive beam passes through the plasma and its self-field is intense enough to expel all the plasma electrons away from the axis, a plasma bubble filled with plasma ions can be formed and moves along with the drive beam (which is the so-called blowout regime) [16]. As a result, the trailing beam will continuously gain energy until the drive beam exhausts its energy and no longer excites the plasma bubble.
In the blowout regime, when the trailing beam is loaded into the plasma wake, the longitudinal electric field of the wake will be modified. When the trailing beam is properly loaded (optimal beam-loading [18]), the longitudinal electric field felt by the trailing beam is locally flattened so that all the contained particles can be accelerated at the same rate resulting in the smallest increase in the energy spread as required by most accelerator applications. This beam-loading effect plays an important role on the beam quality and has been actively studied [18,21,[26][27][28][29][30]. Scaling laws for beam-loading are always useful as the guidance to design the PBA stage efficiently. There were two scaling laws proposed for a laserdriven stage. The number of particles loaded into a 3D bubble wake excited by a laser driver was found to scale with the normalized volume of the bubble or the square root of the laser power [27]. A similar scaling law but with a distinct parameter space was also offered by [26]. However, these scaling laws did not give the exact coefficient and the proper place for loading the trailing beam. In [18], an analytical theory was proposed for beam-loading effect in the blowout regime to maintain the energy spread of the trailing beam. The charge, the shape and the placing of the trailing beam can be estimated for both a laser-driven stage and a beam-driven stage via this theory. However, when designing a two-bunch PWFA stage, the theory provided by [18] is still not easy to use because it lacks the parameters for the drive beam. In addition, this analytical theory was obtained based on the assumption that the maximal normalized bubble radius is much larger than 1. Due to the limitation on the beam peak currents at present PWFA facilities, most PWFA experiments are conducted at a smaller maximal bubble radius, and no analytical model exists to predict their performances. Therefore, we here take a numerical approach to provide fitting formulas for the optimal beam-loading in a data-driven way that will help the design of two-bunch PWFA experiments. The fitting formulas consider parameters for both drive beam and trailing beam. In section 2, the method to find the optimal beam-loading in a two-bunch PWFA stage is discussed. Subsequently, two fitting formulas are given in section 3. Specifically speaking, their availability for trailing beams with a longitudinal flat-top profile or a longitudinal trapezoidal profile are discussed in section 3.3. In section 4, the scaling law for charges of drive beams and trailing beams under the optimal beam-loading is derived from the fitting formulator. In section 5, the optimal plasma densities for the maximum accelerated charge and maximum acceleration efficiency under the optimal beam-loading are discussed. In the last section, we summarize the results presented in this paper.
In a two-bunch PWFA stage, when the blowout occurs, the beam energy spread is mainly affected by the longitudinal wakefield [16]. Thus, having the longitudinal wakefield within the trailing beam as flat as possible is the most effective method to preserve beam energy spread. Parameters including beam charge Q, rms beam length σ z , rms beam spot size σ r , beam separation d and plasma density n p are usually considered in a two-bunch PWFA design. For tri-Gaussian beams, the beam separation is defined as the distance between the center of the drive beam and that of the trailing beam. Electron beams with a tri-Gaussian profile have
, where ξ = ctz is the co-moving coordinate, x and y are the transverse coordinates, and the beam peak density is
where N b is the total number of electrons in the beam [31]. In this paper, we adopt normalized units. The beam density is normalized to the plasma density n p and the charge density is normalized to en p where e is the electron charge. The length is normalized to the plasma skin depth k -1 p ≡ c/ω p , where c is the speed of light and ω p = 4πe 2 n p /m e is the plasma frequency where m e is the electron mass. The electric field is normalized to m e cω p /e. By using normalized units, we can drop the dependency of plasma density to simplify the model. Actually, engineering formulas that take the plasma density into account (described in section 5) can be easily obtained from our fitting formulas in normalized units.
In the blowout regime, if the bubble radius R b is much larger than the rms beam spot size σ r , any variation within the beam spot size for the same charge per unit length Λ = n b σ 2 r will hardly change the wake [32]. In other words, the acceleration structure is determined by Λ as long as R b ≫ σ r and the beam length is fixed. Therefore, we assume the beam has a very small spot size like a delta-function, in which case the dependency of the beam spot size is neglected. The deltafunction-like beam is implemented in the simulation code QuickPIC [33,34] by modifying the subroutine to directly initialize the beam density on the axis, which indicates that the beam has a spot size equal to the transverse cell size as shown in figure 1(a). In this simulation, the simulation box has the size of 8.0 × 8.0 × 10.0 (x, y, ξ) with 512 × 512 × 512 cells. The drive beam has Λ d = 0.2 while the trailing beam has Λ t = 0.16. The length of the drive beam and that of the trailing beam are σ zd = 1.0 and σ zt = 0.25, respectively. The beam separation is d = 4.0. Figure 1(b) shows the comparison of the on-axis E z lineouts from the wake driven by one cell wide beams and beams with σ r = 0.1, and they are almost identical. Therefore, we can ignore the beam spot size and find Λ t for the optimal beam-loading with given Λ d , σ zd , σ zt and d. The goal of the optimization is to achieve the minimal energy spread for the trailing beam in the blowout regime, which requires the trailing beam feels the E z that is as flat as possible in the longitudinal direction. We use the following objective function for the optimization,
where ξ s (ξ e ) is the head (tail) location of the trailing beam, λ bt (ξ) = ´ρbt (x, y, ξ)dxdy is the normalized charge per unit length of the trailing beam, ρ bt is the normalized charge density of the trailing beam and Λ t is the peak value of λ bt (ξ). F(Λ t ) is the mean square deviation of weighted on-axis E z , where the density profile of the trailing beam is used as the weight. This is a single-objective optimization [35] process because we aim to find the minimum of F(Λ t ) while changing Λ t . By doing several tests, we find the optimization is a typical convex optimization [36], in which for any two points Λ t1 , Λ t2 in the domain of Λ t and m ∈ (0, 1) we have
). For a convex optimization, the local optimum is the global optimum, and the extreme value is the optimal solution [36]. Thus, the local optimization algorithm can be applied.
To achieve high performance, we optimize the F(Λ t ) with the Broyden-Fletcher-Goldfarb-Shanno (BFGS) algorithm [37], which has been extensively used to solve nonlinear optimization problems and has been considered to be the most effective of all quasi-Newton methods [38][39][40][41][42]. We set Λ t = Λ d as the initial solution for the optimization process. By assuming the wakefield does not evolve, the objective function can be evaluated from one-time-step QuickPIC simulation result (i.e. the static wakefield).
A typical optimization result is shown in figure 2. In this example, beam parameters are Λ d = 1.0, σ zd = 1.0, σ zt = 0.25 and the beam separation is d = 4.5. We plot the on-axis E z at different Λ t . The plasma and beam densities are just for illustration, and they do not vary. As shown in figure 2, with the optimal Λ t = 1.49 the trailing beam feels a more flat E z than that with the initial Λ t = 1.0 we used. The E z at the optimal beam-loading is a little overloaded compared with that of Λ t = 1.2, in which the ξ derivative of E z only has one zero point within the trailing beam. This is because the trailing beam has a Gaussian profile and the optimal beam-loading case will generate a smaller rms energy spread. To verify the result obtained from the BFGS algorithm, we manually do a parameter scanning for Λ t from 0.1 to 4.0 with a step size of 0.01. The Λ t for the optimal beam-loading agrees very well with the result from BFGS algorithm. The relative difference between them is about 0.02%. With the BFGS algorithm and QuickPIC simulation, the case shown above requires 16 evaluations by QuickPIC to find the optimal Λ t , and the total computing time is 7 min with 64 cores. We then perform longdistance accelerations. We find that the energy spread of the trailing beam is 1.69% at Λ t = 1.49 which is smaller than that with 2.35% at Λ t = 1.2 with the same initial energy (10 GeV) and the same energy gain (about 7.3 GeV). This comparison result agrees well with our optimization. We note that it used to be a common sense that the case of red line in figure 2 would have the smallest rms energy spread. This is not true because the E z for that case is monotonically decreasing while the black line in figure 2 is not. As a result, the case of the black line may let more beam particles have the same energy gain at different longitudinal locations, and finally have a smaller rms energy spread than the case of red line.
A Python program is developed to automatically optimize a large number of parameter sets of (Λ d , σ zd , σ zt , d) (see appendix A for details). In these sets of (Λ d , σ zd , σ zt , d), the Λ d has a range of [0.0144, 7.70] and the σ z for both beams has a range of [0.0952, 1.90]. These ranges basically cover the parameters of FACET [43], FACET II [22], FLASHForward [24] and other facilities [20] with a plasma density of 10 16 cm -3 .
where R bmax ≃ 2 √ Λ d gives a good estimate of the maximum bubble radius [31], in order to have the trailing beam be approximately located inside the first plasma bubble. Once the ranges of Λ d , σ zd , σ zt and d are determined, we evenly select values within the range of each parameter. In addition, we also need to ensure that settings for each Quick-PIC simulation are appropriate (see appendix B for details). In each optimization process, we dump the Λ d , σ zd , σ zt , d, the optimal Λ t , the maximum decelerating wakefield W dec inside the drive beam, the averaged accelerating wakefield felt by the trailing beam W acc = ´ξe ξs E z (ξ)λ bt (ξ)dξ/ ´ξe ξs λ bt (ξ)dξ and the transformer ratio R = |W acc /W dec |.
Data from the automatic optimization will have some bad parameter sets, i.e. the outliers. For example, some datasets have the trailing beam too far away from the drive beam so that it cannot be effectively accelerated even though the optimization process succeeds. Therefore, we use the boxplot method [44] and standard normal distribution method [45] to eliminate these outliers. We finally obtain 8537 sets of data for the optimal beam-loading database. The average time for each optimization is only around 7.6 min with 64 cores. The range of Λ d , σ zd , σ zt , d and Λ t is presented in table 1. Note that table 1 shows the global range for the beam separation. The actual range of the beam separation varies according to the beam parameters.
We use the data-driven method to solve the optimal beamloading problem in the blowout regime. To obtain explicit fitting formulas, we use the Python library scikit-learn [46] to carry out polynomial regression, which can be generalized into the linear regression [47].
During the process, the data are split into several random but with general equal-size folds. And we set some of them as the training dataset and the remaining as the test dataset. Then constructing polynomial features is demanded because the degree of polynomial features we choose directly affects the goodness of fit. Here, we use the coefficient of determination r 2 [46] to measure how well unseen test dataset tends to be predicted by the model. The closer r 2 is to 1, the better the goodness of fit is. To determine the best choice of degree, we use the k-fold cross-validation method to evaluate our model to avoid over-fitting [46]. It divides the training dataset into k subsets at once and then trains a model k times in total. In each model training, we use k -1 subsets to train the model and use the remaining one to validate the model and obtain the r 2 for each training. The averaged r 2 is obtained at the end of this loop for a particular degree. And the best degree should have the largest averaged r 2 with this k-fold cross-validation method. As a common choice, we choose k = 10 for our calculation. After determining the best degree, we use the whole training dataset to train a model (i.e. get the fitting formula) and use the test dataset to do the final evaluation.
3.2. The fitting formulas for Λ t and R By using the method described above, we can obtain the fitting formula for the optimal Λ t , which can be written as Λ t = f(Λ d , σ zd , σ zt , d). More specifically, training dataset and test dataset account for 75% and 25% of the database, respectively. When we use training dataset to perform 10-fold crossvalidation, we obtain the averaged r 2 ≃ 0.999 at degree of 3, which is larger than those at other degree values. Therefore, we use the whole training dataset to do the polynomial regression at degree of 3 and obtain r 2 ≃ 0.999 when evaluating the test dataset. This represents high prediction accuracy. The final result of the polynomial regression, i.e. the fitting formula for Λ t , is
Table 2. Fitting coefficients for the fitting formula of Λt. where the fitting coefficients are given in table 2. Besides Λ t , the transformer ratio R is also an important parameter we concern in a two-bunch PWFA stage. We consider that R is dependent on Λ d , Λ t , σ zd , σ zt and d. Following the same procedure, we can get the explicit expression of R = f(Λ d , Λ t , σ zd , σ zt , d). In this case, training dataset and test dataset comprise 80% and 20% of the whole database, respectively. We finally choose the degree of 2, with which we get the highest averaged r 2 ≃ 0.98 when performing 10-fold crossvalidation. In the final evaluation using the test dataset, we get r 2 ≃ 0.99, which represents high prediction accuracy. The fitting formula for R is
where the fitting coefficients are given in table 3.
Through the fitting formulas, we can obtain the optimal Λ t without running the optimization program. For example, for Λ d = 1.0, σ zd = 1.0, σ zt = 0.2 and d = 4.0, equation (2) gives the optimal Λ t = 1.652, while the optimization program gives Λ t = 1.644. The results agree well with each other. When calculating the transformer ratio R using the fitting formula, we first need to obtain the optimal Λ t through equation (2), and then substitute the optimal Λ t into equation ( 3) to obtain R. This gives R = 0.622 in this case, while the optimization program gives R = 0.622. They still agree very well with each other. In figure 3, we compare more results from the optimization program with the results given by the fitting formulas. The green solid line in figure 3(a) plots the optimal Λ t versus d with Λ d = 1.0, σ zd = 1.0 and σ zt = 0.2 by using the fitting formula equation (2). The blue cross points are the results from the optimization program and they agree very well with the fitting results. The pink and black solid lines in figure 3(a) have different Λ d but the same σ zd and σ zt , and they agree very well with the results from the optimization program. We also change σ zt and Λ d while keeping σ zd and still find good agreements between the fitting results (dashed lines) and the optimization results (dot points) as shown in figure 3(a). Furthermore, we calculate the transformer ratio R from the fitting formula equation ( 3) and show the results in figure 3(b), which has another axis of R than figure 3(a). The fitting results also agree very well with the optimization results. From the results shown in figure 3, we can also find that for given Λ d , σ zd and σ zt , the bigger the d is, the smaller the Λ t is and the higher the R is, which agrees with the understanding of beam-loading in the nonlinear plasma wake [18]. The applicable parameter range for these two fitting formulas is listed in table 1. In addition, the beam energy had better to be larger than 100 MeV when using these fitting formulas.
We also test the availability of the fitting formulas for trailing beams with a longitudinal flat-top profile or a longitudinal trapezoidal profile. We pick up three tri-Gaussian cases with the same drive beam parameters and the same σ zt = 0.190 but different d. We plot the on-axis E z of the plasma wake in figure 4(a). In these simulations, the drive beam has Λ d = 0.918 and σ zd = 0.952, and its beam center is located at ξ d = 3.33. For each simulation as shown in figure 4 4(a), all these three cases have reached the optimal beam-loading. When switching them to the longitudinal flat-top profile, we keep Λ t and the total particle number the same as those of tri-Gaussian trailing beams. Therefore, the flat-top beam length should be l zF = √ 2πσ zt . We load these flat-top beams with their heads at a distance √ 2σ zt in front of ξ 1,2,3 in order to maintain the transformer ratio (as suggested in [18]). As shown in figure 4(b), the beam-loading effect of flat-top trailing beams mimics that of tri-Gaussian trailing beams. In [18], it is shown that the trapezoidal trailing beams can perfectly flatten the E z . For trapezoidal trailing beams, we still keep the total particle number and maximal Λ t the same as those of tri-Gaussian trailing beams. The trapezoidal beam also has a sharp edge as the flat-top beam. Thus, we load trapezoidal beams at ξ1,2,3 = ξ 1,2,3 -√ 2σ zt . The slope of the trapezoidal profile a equals to E z where the beam-loading starts [18], which roughly equals to the averaged accelerating wakefield of the tri-Gaussian beam. For three trapezoidal trailing beams plotted in figure 4(c), we have a 1 = -0.539, beam length l z1 = 0.562 at ξ1 , a 2 = -0.709, l z2 = 0.609 at ξ2 and a 3 = -0.932, l z3 = 0.716 at ξ3 , where the beam length is derived from the total charge of the beam. As shown in figure 4(c), E z is almost flattened and the transformer ratio is well maintained. Therefore, through proper beam parameter transformations, the fitting formulas of tri-Gaussian beams can still give a good estimation for flat-top or trapezoidal trailing beams.
Not only can the fitting formulas be used to find particular beam parameters for the optimal beam-loading, they can also unveil many physics features under the optimal beam-loading. One of the features is the relation between the charge of the drive beam and that of the trailing beam under the optimal beam-loading. The beam charge is proportional to Λσ z . Therefore, by multiplying σ zt on both sides of equation ( 2) and rearranging the right hand side of the equation, we can find the relation between Λ t σ zt and Λ d σ zd as where
According to equation ( 4), with σ zd = 1.0, σ zt = 0.2 and d = 4.0 as an example, we can get Λ t σ zt = -0.0004781(Λ d σ zd ) 3 -0.01231(Λ d σ zd ) 2 + 0.3835(Λ d σ zd ) -0.04041, which is plotted as the blue line in figure 5(a). From the plot, we can find that Λ t σ zt almost increases linearly with Λ d σ zd . This is because the high order terms are much less than the Λ d σ zd term in this example. Therefore, equation ( 4) can be reduced to Λ t σ zt = D(Λ d σ zd ) + G. This means that once the optimal beam-loading is reached, it is always satisfied when increasing charges of both beams with the same ratio D. In figure 5(a), we plot three other lines with different d or σ zt . And they all obey the simple scaling law Λ t σ zt = D(Λ d σ zd ) + G, where D and G depend on σ zd , σ zt and d. If G is much less than DΛ d σ zd , we can further neglect G and Λ t σ zt will become proportional to Λ d σ zd . This means that we can change charges of both beams at the same rate without breaking the optimal beam-loading condition. In addition, with equation (3) we can also calculate the transformer ratio R for the lines in figure 5(a), which is shown in figure 5(b). This will bring much convenience for designing a two-bunch PWFA stage.
So far, we are using the normalized units for each parameter. This means the physics features we obtained in the last section is only available for a fixed plasma density. However, we are also interested in how the beam parameter varies as the plasma density changes under the optimal beam-loading. This can be obtained by switching the normalized units in the equation back to the original ones. We have the charge of the drive (trailing
Therefore, equations ( 2) and ( 3) can be converted to equations that have the plasma density as an additional variable (see appendix C and equation ( 6) for details). Here we will focus on how the plasma density will affect equation (4). We convert equation ( 4) into an engineering formula
where
The coefficients in equation ( 5) are listed in table 4.
The equation (5) shows the relation between the charge of the trailing beam and the plasma density under the optimal beam-loading. For example, when Q d = 1.5 nC, L d = 60 µm, L t = 12 µm and l = 300 µm, we can obtain
p [10 16 cm -3 ] -0.3230, which is plotted as the blue line in figure 6(a).
The plot shows an interesting feature that the charge of the trailing beam has a maximal value when the plasma density varies under the optimal beam-loading. For this case, the (5). trailing beam reaches its maximal charge Q tmax = 0.440 nC at an optimal plasma density n pQ = 2.946 × 10 15 cm -3 , which is marked as the blue dashed line in figure 6(a). We plot the plasma wake and the on-axis E z for the same case in figure 6(b)(1), and we can see that the trailing beam reaches the optimal beam-loading. In figure 6(a), we also plot equation (5) as three other lines with different Q d or L t . Parameters for each case are listed in table 5. They all show that there is an optimal plasma density n pQ (marked as the dashed lines) for obtaining the maximum accelerated charge Q tmax . The values of n pQ and Q tmax for each case are also listed in table 5. For each case, we plot the plasma wake and the on-axis E z at the optimal plasma density in figure 6(b).
Although the trailing beam reaches its maximum charge, the transformer ratio in each case is low (less than 1) as shown in figure 6(b). In other words, the acceleration efficiency is low for these cases. Actually, it is easy to find how the acceleration efficiency varies with regard to the plasma density. The acceleration efficiency can be calculated through
By switching the units back to the Fitting coefficients in equation (6).
original ones in equation ( 3) and substituting equation ( 5) into it, we can have an engineering formula of R that depends on
where
The coefficients are given in table 6.
Then by substituting equation ( 5) together with equation ( 6) into the equation of η, we can have
where X 1 = OH, X 2 = (OM + YH), X 3 = (OP + YM + ZH), X 4 = (OS + YP + ZM + CH), X 5 = (YS + ZP + CM + KH), X 6 = (ZS + CP + KM + WH), X 7 = (CS + KP + WM + TH), X 8 = (KS + WP + TM), X 9 = (WS + TP) and X 10 = TS.
In figure 6(c), we plot η versus n p with four sets of Q d , L d , L t and l, which are the same as those in figure 6(a). There is also an optimal plasma density (marked as the dot-dashed lines) for obtaining the maximum η under the optimal beam-loading. Note that η becomes negative at lower n p because the beam separation is so small that the trailing beam is located in the decelerating phase in the plasma wake. Table 5 also lists the optimal plasma density n pη for the maximum acceleration efficiency η max and Q t at η max . Figure 6(d) shows the plasma wake and the on-axis E z at the optimal n p for the maximum acceleration efficiency for each case in figure 6(c). We can see that trailing beams are all located at the back of the bubble, which ensures that the transformer ratio is close to or larger than 1. By comparing figures 6(a) and (c), we can see that the optimal plasma densities for maximum accelerated charge and maximum acceleration efficiency are usually different. This means that for given Q d , L d , L t and l, we have to make a compromise between having the maximum accelerated charge and having the maximum acceleration efficiency when choosing the plasma density. In order to do that, for example, we can choose the value in the middle of two optimal plasma densities. In addition, the curves shown in figure 6(a) also indicate that the optimal beam-loading condition cannot hold for fixed beam parameters at different plasma densities. Therefore, additional energy spread will be induced in the region where the plasma density varies (e.g. the plasma density ramps).
By using the BFGS optimization method and the quasi-static code QuickPIC, we obtain a large amount of optimal beamloading cases of two-bunch PWFA in a wide parameter range. Then we derive two fitting formulas from these data by using the polynomial regression with 10-fold cross-validation method. One fitting formula can find the optimal Λ t under the optimal beam-loading condition with given Λ d , σ zd , σ zt and d. The other one can find the transformer ratio with given Λ d , σ zd , Λ t , σ zt and d under the optimal beam-loading condition. We use the normalized units in these two fitting formulas that makes them not have the dependency of the plasma density. One can easily transform the fitting formula into an engineering equation that has the plasma density as a variable (shown as equation (C1) and equation ( 6)). The fitting formulas agree with the simulation results very well. It is a very efficient tool for obtaining the optimal beam-loading parameters when designing a PWFA stage using two tri-Gaussian electron beams in the blowout regime. We also test the fitting formulas with trailing beam that has a flat-top or trapezoidal longitudinal profile. The fitting formulas can still give a good estimation after the simple parameter transformation between different longitudinal profiles.
We explore new physics features of the optimal beamloading based on the fitting formulas. One feature is that once the optimal beam-loading is reached, it is always satisfied when we increase the charges of drive beam and trailing beam at the same ratio. This ratio is dependent on the length of drive and trailing beams and the beam separation. Another physics feature is that under the optimal beam-loading condition there are two optimal plasma densities for the maximum accelerated charge and the maximum acceleration efficiency with given parameters of the drive beam, the length of the trailing beam and the beam separation. These two features provide an important guidance for the two-bunch PWFA design.