{'fraction_non_alphanumeric': 0.03934050999606595, 'fraction_numerical': 0.017309824398269016, 'mean_word_length': 3.938537619215825, 'pattern_counts': {'":': 0, '<': 3, '<?xml version=': 0, '>': 0, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 11, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
The Optimal Beam-loading in Two-bunch Nonlinear Plasma Wakefield Accelerators

Due to the highly nonlinear nature of the beam-loading, it is at present not possible to analytically determine the beam parameters needed in a two-bunch plasma wakefield accelerator for maintaining a low energy spread. Therefore in this paper, by using the Broyden-Fletcher-Goldfarb-Shanno algorithm for the parameter scanning with the code QuickPIC and the polynomial regression together with k-fold cross-validation method, we obtain two fitting formulas for calculating the parameters of tri-Gaussian electron beams when minimizing the energy spread based on the beam-loading effect in a nonlinear plasma wakefield accelerator. One formula allows the optimization of the normalized charge per unit length of a trailing beam to achieve the minimal energy spread, i.e. the optimal beam-loading. The other one directly gives the transformer ratio when the trailing beam achieves the optimal beam-loading. A simple scaling law for charges of drive beams and trailing beams is obtained from the fitting formula, which indicates that the optimal beam-loading is always achieved for a given charge ratio of the two beams when the length and separation of two beams and the plasma density are fixed. The formulas can also help obtain the optimal plasma densities for the maximum accelerated charge and the maximum acceleration efficiency under the optimal beam-loading respectively. These two fitting formulas will significantly enhance the efficiency for designing and optimizing a two-bunch plasma wakefield acceleration stage.


I. INTRODUCTION
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 inside the plasma wake can easily exceed 10 GeV/m [5,[8][9][10][11], which are orders of magnitude higher than that of conventional accelerators. Such high acceleration gradients can significantly reduce the size and the cost of accelerators. This makes the PBA a promising candidate for the future linear colliders or light sources. Recently, the particle-beam-driven plasma wakefield acceleration (PWFA) has attracted a lot of attention due to tremendous theoretical and experimental progress [5][6][7][12][13][14][15][16][17][18][19][20][21]. There have been increasing numbers of 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 laser-driven 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 Ref. [26]. However, these scaling laws did not give the exact coefficient and the proper place for loading the trailing beam. In Ref. [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 Ref. [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 Sec. II, the method to find the optimal beam-loading in a two-bunch PWFA stage is discussed. Subsequently, two fitting formulas are given in Sec. III.
Specifically speaking, their availability for trailing beams with a longitudinal flat-top profile or a longitudinal trapezoidal profile are discussed in Sec. III C. In Sec. IV, the scaling law for charges of drive beams and trailing beams under the optimal beam-loading is derived from the fitting formulator. In Sec. V, 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.

A. Optimization of Beam Parameters
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 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 Sec. V) 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 delta-function-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 Fig. 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. Fig. 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) . 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 onetime-step QuickPIC simulation result (i.e. the static wakefield).
A typical optimization result is shown in Fig. 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 Fig. 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 minutes with 64 cores. We then perform long-distance 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 Fig. 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 Fig. 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.

B. Large-range Parameter Scanning for Optimal Beam-loading
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 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 QuickPIC 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 minutes with 64 cores.
The range of Λ d , σ zd , σ zt , d and Λ t is presented in Table I. Note that Table I 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 beam-loading 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.
B. 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 cross-validation, we obtain the averaged r 2 0.999 at degree of 3, which is larger than that 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 where the fitting coefficients are given in Table II.
Besides Λ t , the transformer ratio R is also an important parameter we concern in a twobunch 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 cross-validation. 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 III.   Fig. 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 Fig. 3 (a). Furthermore, we calculate the transformer ratio R from the fitting formula Eq. (3) and show the results in Fig. 3 (b), which has another axis of R than Fig. 3 (a). The fitting results also agree very well with the optimization results.
From the results shown in Fig. 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 I. In addition, the beam energy had better to be larger than 100 MeV when using these fitting formulas.

C. Flat-top and Trapezoidal Trailing Beams
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 Fig. 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 Fig. 4 (a), a tri-Gaussian trailing beam is loaded at three distinct   Fig. 4 (a), all these three cases have reached the optimal beamloading. 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 Ref. [18]). As shown in Fig. 4 (b), the beam-loading effect of flat-top trailing beams mimics that of tri-Gaussian trailing beams. In Ref. [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 Fig. 4 (c

MAL BEAM-LOADING
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 beamloading. 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 Eq. (2) and rearranging the right hand side of the equation, we can find the relation between Λ t σ zt and Λ d σ zd as which is plotted as the blue line in Fig. 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, Eq. 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 Fig. 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 Eq. (3) we can also calculate the transformer ratio R for the lines in Fig. 5 (a), which is shown in Fig. 5 (b). This will bring much convenience for designing a two-bunch PWFA stage. where The coefficients in Eq. (5) are listed in Table IV.    Fig. 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 trailing beam reaches its maximal charge Q tmax = 0.440nC at an optimal plasma density n pQ = 2.946 × 10 15 cm −3 , which is marked as the blue dashed line in Fig. 6 (a). We plot the plasma wake and the on-axis E z for the same case in Fig. 6 (b)(1), and we can see that the trailing beam reaches the optimal beam-loading. In Fig. 6 (a), we also plot Eq. (5) as three other lines with different Q d or L t . Parameters for each case are listed in Table V. 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 V. For each case, we plot the plasma wake and the on-axis E z at the optimal plasma density in Fig. 6 (b).
Although the trailing beam reaches its maximum charge, the transformer ratio in each case is low (less than 1) as shown in Fig. 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 where O = m 20 H 2 /L 2 The coefficients are given in Table VI.
In Fig. 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 Fig. 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 V also lists the optimal plasma density n pη for the maximum acceleration efficiency η max and Q t at η max . Fig. 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 Fig. 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 Fig. 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 Fig. 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).

VI. CONCLUSION
By using the BFGS optimization method and the quasi-static code QuickPIC, we obtain a large amount of optimal beam-loading 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 Eq. (C1) and Eq. (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 beam-loading 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 round (3.5σ zd , 2), σ zd ≥ 0.3 and 4σ zt < 3.5σ zd + d, ceil (4σ zt − d), 4σ zt ≥ 1 + d or 4σ zt ≥ 3.5σ zd + d, 1, σ zd < 0.3 and 4σ zt < 1 + d, where function round (x, m) is to round x to a specific precision m in decimal digits and function ceil (x) is to return the ceiling of x as an integral. Subsequently, the center of the trailing beam C t is set to Following this, the length of box in the longitudinal direction box z can be set to box z = max {4, round (max (C t + 5σ zt , C d + 3.5σ zd ), 2)}, and the cell number in the longitudinal direction is where I = s 25 L 3 d + s 26 L 2 d L t + s 27 L 2 d l + s 28 L d L 2 t + s 29 L d L t l + s 30 L d l 2 + s 31 L 3 t + s 32 L 2 t l + s 33 L t l 2 + s 34 l 3 , J = s 9 L 2 d +s 10 L d L t +s 11 L d l +s 12 L 2 t +s 13 L t l +h 14 l 2 +s 19 Λ d L 2 d +s 20 Λ d L d L t +s 21 Λ d L d l + s 22 Λ d L 2 t + s 23 Λ d L t l + s 24 Λ d l 2 , U = s 2 L d + s 3 L t + s 4 l + s 6 Λ d L d + s 7 Λ d L t + s 8 Λ d l + s 16 Λ 2 d L d + s 17 Λ 2 d L t + s 18 Λ 2 d l, V = s 1 Λ d + s 5 Λ 2 d + s 15 Λ 3 d + s 0 and L d , L t and l are in the unit of µm. The coefficients are given in Table VII.