Numerical calculation of dynamical friction in electron cooling systems, including magnetic field perturbations and finite time effects

The orders-of-magnitude higher luminosities required by future electron–ion collider concepts require a dissipative force to counteract the numerous factors acting to gradually increase the phase space volume of relativistic ion beams. High-energy electron cooling systems could provide the necessary dissipation via dynamical friction, but will have to be designed for new parameter regimes. It is expected that magnetic field errors, finite interaction time and other effects will reduce the dynamical friction and hence increase the cooling time, so improved understanding of the underlying dynamics is important. We present a generalized form of the classical field-free friction force equation, which conveniently captures some of these effects. Previous work (Bell et al 2008 J. Comput. Phys. 227 8714) shows both numerical and conceptual subtleties associated with undersampling of strong collisions, and we present a rigorous mathematical treatment of such difficulties, based on the use of a modified Pareto distribution for the electron–ion impact parameters. We also present a very efficient numerical algorithm for calculating the dynamical friction on a single ion in the field free case. For the case of arbitrary magnetic field errors, we present numerical simulation results, showing agreement with our generalized friction force formula.


3
Although numerical simulation can provide an answer for a particular set of parameters, it is important to know how friction force generally depends on the frequency and amplitude of the error fields. Such a formula is presented in section 3. A modification of the classic equation, this formula is motivated by and agrees with the results of simulations presented here.
Our previous work [14] has shown the importance of finite time effects in simulating or otherwise calculating the dynamical friction force. There we obtained an approximate formula for the effective values of minimal and maximal impact parameters reduced by the finite time effects. Section 4 of the present paper goes further by providing an efficient algorithm for numerical calculation of the dynamical friction in a finite time interval. This new algorithm correctly samples the multidimensional space to obtain accurate results with the minimum possible computational effort. Also, the inherent asymmetry of finite time electron trajectories is included in the calculation. This algorithm also calculates the contributions to the friction force from different impact parameters, and estimates the computation error. This method, however, does not capture the effect of external electromagnetic fields.
Initial study [14] indicated that direct numerical simulations have shown a discrepancy with tested analytic formulae because of the so-called under-sampling effect. The fact that there are only a few strong collisions, which however contribute significantly to the average, makes the friction force hard to estimate accurately and shows a close analogy with the Pareto distribution. In section 5, we present a rigorous statistical treatment of electron trajectories with small impact parameters and show that for a realistic number of collisions, the change in the ion's velocity can significantly deviate from its average. The whole distribution must be considered to obtain an accurate estimate of the cooling rate. The quantiles of this distribution for different numbers of collisions is numerically calculated.
In this paper, we use abbreviations and notation common to probability theory books. We use P(E) to denote a probability of event E, which is often in the form of an inequality. The cumulative distribution function (CDF) of a scalar random variate (RV) ξ is defined as P(x) = P(ξ < x), and the corresponding probability density function (PDF) is defined as p(x) = d dx P(x). For the average and variance of ξ , we use the operator notation E(ξ ) and Var(ξ ), which in the 1D case are defined as Similarly, P(E|C), E(ξ |C) and Var(ξ |C) denote the probability, average and variance under the condition of event C. In multifold integrals, we often place the infinitesimal right after the integral sign to indicate the correspondence between the integration variable and its range. For example, we may use the notation

Assumptions of the classic formula for the friction force
We start with a sketch of the derivation of the classical formula and its assumptions. The details of this derivation can be found in [6]. This will serve as a starting point for the discussion; the results obtained in the following sections of this paper will rely on much weaker assumptions.
Following [6], let us consider an ion moving through some volume V uniformly filled with electrons. The velocities of the electrons have Gaussian distribution where σ x , σ y and σ z can be different. Throughout this paper, the velocities are given in the frame of the electron beam, in other words, in the frame in which the average electron velocity is (0, 0, 0). The interaction is non-relativistic, and the particles interact only via the Coulomb force.
We assume that the collisions involving three bodies are negligibly rare, and the stopping effect is the result of individual binary collisions, which occur independently. The friction force is defined as the average rate of momentum transfer over a time interval of duration T , where v ion (x, v) denotes the change of the ion's velocity v ion due to an individual collision with an electron which has initial velocity v and position x. Here, E (· · ·) denotes the averaging operator over x in some volume V ⊂ R 3 and all possible velocities v ∈ R 3 , and n e is the density of the electrons. Using the law of momentum conservation where v is the change in the electron's velocity during the interaction, we transform equation (2) to Expanding the averaging operator, we obtain As it will be evident from further derivation, the spatial 3D integral diverges if it is taken over all R 3 . Therefore, to proceed further, we must restrict the volume over which this integral is taken. It is convenient to choose the volume depending on v, and describe it in a cylindrical coordinate system chosen such that the z-axis is parallel to (v − v ion ).
We assume that the volume is axisymmetric and ϕ ranges in [0, 2π ]. The component of F perpendicular to (v − v ion ) is zero because of the symmetry 4 , therefore we could replace v with v || , which is the component of v parallel to (v − v ion ) and independent of ϕ. For now, we assume that the range of z is [L 1 , L 2 ] and the range of r is [0, ρ max ]. The volume defined by these limits is shown in figure 1. Equation (6) transformed into the cylindrical coordinate system becomes The expression in brackets is generally treated under the assumption that the domain and interaction time are 'large'. In mathematical terms, this means taking limits T → ∞, L 1 → −∞, L 2 → ∞ and ρ max → ∞. As we will show next, the result depends on the order in which these limits are taken. Therefore, this operation requires rigorous mathematical treatment. To our knowledge, such treatment has not yet been published. Indeed, if the first limit is taken in time, then the result is zero: In the applications relevant to the electron cooling, the interaction time is short and comparable to the plasma period. Hence, the Debye shielding does not occur, and effectively the domain is unbounded. Therefore, for this application, it is appropriate to take the limits first in L 1 → −∞, L 2 → ∞ and ρ max → ∞. This indeed leads to a finite non-zero limit, This limit, however, does not linearly depend on T . Moreover, if the limit T → ∞ is taken of expression (9), then the result is infinity according to [14]. In [6], the momentum transfer in a single electron-ion collision is estimated by This estimate is obtained under the assumption that the electron starts infinitely far from the ion, the particles interact for infinite time, and the trajectory of the electron is only slightly deflected from a straight line. Taking into account the rate of the collisions, which is proportional to n e (v − v ion ) || , the classic estimate on (9) is obtained, The frame of reference is chosen such that the ion is at rest.
Such a derivation would be justified if the duration of the collisions were very short and much smaller than the average time between the collisions. In this case, however, each collision takes an infinite amount of time, so a more careful treatment is required.
To obtain the classic result, we need to choose differently the volume V that initially contains the electrons. It must be at some distance from the ion, as shown in figure 2, and it is defined as in the following. For simplicity, let us consider the system of coordinates in which the ion is at rest, and we only consider electrons that have velocity v. Imagine that the particles evolve for some time T 1 and then they end up in some different volume V . The location of the volume V is chosen so that V and V are roughly symmetric around the ion, so the volume V is roughly at the distance |v|T 1 /2 from the ion. This fulfills the assumption of the classic approach that the collisions are symmetric. Now the limit T 1 → ∞ can be taken, and equation (10) is justified for each particle inside the volume. In the resulting expression, the limit T 2 → ∞ can be taken, which corresponds to increasing volume V . This will result in a factor n e (v − v ion ) || , and the classic estimate (11) is obtained.
The limit in T 1 and T 2 can be taken simultaneously if some dependence between T 1 and T 2 is introduced such that T 1 T 2 ; for example, T 1 = T 2 2 /T * for some constant T * . Admittedly, this is an artificial setup. However, it has been verified empirically that (11) can lead to a reasonable estimate, which has been used for decades in a wide range of physical systems. In the remainder of this section and the next section, we proceed with the calculations using this estimate. The more natural setup with a finite volume chosen around the ion, as shown in figure 1, will be considered in section 4, where we will present the results of the simulations that do not rely on the numerous assumptions of the classic approach.
The change in electron velocity in a single ion-electron interaction v || has to be integrated over all impact parameters r with the Jacobian factor r . Therefore, the expression to be integrated is r v || . Its dependence on r is shown in figure 3. Evidently, for a wide range of impact parameters r , r v || depends on r approximately as 1/r . Using these estimates, equation (7) is transformed into  To avoid the logarithmic singularities, the range of r had to be limited with ρ min and ρ max . The value of ρ min is chosen such that the scattering angle is 90 • . It can be calculated using The integral over r can be taken, and we obtain the classic formula for the friction force, The choice of ρ max is problematic. It is often set to be the Debye shielding length λ D , which is only an order-of-magnitude estimate. In fact, the Debye shielding occurs gradually, and the classic approach can only provide a crude estimate. In addition, for parameters related to high-energy electron cooling, the interaction time is so short that the Debye shielding does not occur. We will address these limitations of the classical formula in section 4, where we will present an algorithm that computes (7) exactly for finite T , ρ max , L 1 and L 2 . In the following section, we will address a different limitation and estimate the effect of electromagnetic fields on the friction force.

Friction force in the presence of an electromagnetic field
We found that the classical formula (14) can be modified to capture the effect of an error field. The friction force predicted by this formula was compared to numerical simulations (figure 4) Electron density n e 9.5 × 10 13 (19). and found to be in good agreement in a parameter regime typical of an electron cooling section (table 1). We used the VORPAL framework [18] to model electron-ion interactions using a semianalytic binary Coulomb collision (BCC) algorithm [14]. When interaction times were long enough that Debye shielding was not negligible, the electrostatic particle-in-cell (PIC) algorithm for electron-electron interactions was included via second-order operator splitting techniques.
In the laboratory frame, the error field is time-independent magnetic induction B varying only along the direction of electron motion. We chose a coordinate system such that the z-axis is aligned with this direction. The magnitude of the magnetic field imperfections is measured by M defined as To clearly see the friction force dependence on the spectrum of |B|, we include only one wavelength, Transforming this field into the frame of the beam traveling along the z-direction with a relativistic speed [19], we obtain where E 0 = γβcB 0 and γ is the relativistic Lorentz factor. Because in the beam frame the velocity of the electrons are non-relativistic, the Coulomb force dominates the Lorentz force.
Since the volume is small, the dependence of E on z is neglected. Such a field makes particles move in circles with radius ρ osc at speed v osc , Evidently, the collisions with impact parameters smaller than 2ρ osc are affected by the circular motion, which essentially increases the velocity of the transverse motion, thus decreasing the friction force. For the collision with impact parameters larger than 2ρ osc , the trajectories of the electrons are only slightly perturbed helices. A perturbative approach can be used to estimate the strength of the collision in a way similar to the classical approach. Thus the contribution of the weak collisions is not affected by the present field. This intuition is formalized by the following formula, where the large impact parameter F LIP,z and small impact parameter F SIP,z parts of the force are defined as Equation (20) is obtained from (14) by replacing ρ min with 2ρ osc . Equation (21) is obtained from (14) by replacing ρ max with 2ρ osc and using modified velocity distribution with PDF adjusted by adding v osc to the transverse drift, Note that if p(v) were in place ofp(v), the integrals would factor out and the logarithms would combine, reducing (19) to (14). The Coulomb logarithm represents the integral over impact parameter, and splitting it into two parts means that we treat the collisions with impact parameters smaller than 2ρ osc differently from those with the larger ones. The collisions with the smaller impact parameters are affected by the oscillatory motion, which effectively folds into the temperature, and the PDFp( → v ) gets adjusted accordingly. The collisions with the impact parameters larger than 2ρ osc are not significantly affected by the oscillatory motion. Therefore the first integral contains the original PDF. To test this formula, its predictions were compared with first-principle numerical simulations. The friction forces were calculated for several wavelengths and three different field strengths, and good agreement was found with parallel VORPAL simulations (see figure 4). The error bars are of the size of ±1 rms calculated among eight simulated ions. The magnitudes of F LIP and F SIP relative to each other are shown in figure 5. They are comparable, therefore neither part (20) nor part (21) is redundant.

Efficient calculation of the friction force in the field free case
We found that, in the absence of the external magnetic field, the friction force can be calculated very efficiently without making any assumption on ρ min and ρ max . In contrast with the classic formula, the electron motion is not assumed to be only slightly deflected from straight lines; instead the trajectories of the electrons are calculated exactly by solving a two-body problem. Since there are no problems with integrating over a singularity, we do not have to impose any restriction on the smallest impact parameter. Also, there are no assumptions that collisions happen for infinite time, instead the collisions are modeled in a finite time interval [0, T ].
The developed algorithm calculates the contribution to the friction force from individual electron-ion collisions, bins them by different impact parameters, and estimates the error by calculating the rms among collisions within each bin. The values of these contributions are important for cooling in accelerators as well as choosing the domain size for BCC simulations. The results of two simulations with 10 5 and 10 6 collisions are plotted in figure 6 with error bars of ±1 rms.
Because rare collisions have very significant contributions, a particle-based simulation produces noisy results, unless the number of particles is large. Even if the number of particles is large, such simulations spend a large fraction of the total time calculating small, predictable contributions from weak collisions. Our algorithm uses the importance sampling method instead, which gives tremendous advantages over the straightforward approach. For the parameters listed in table 1, the friction force can be calculated with 1% precision in less than 15 s on a desktop with a single 1.2 GHz processor, while it took 2.75 h on 16 nodes of Franklin 5 to obtain a single data point, shown in figure 4. The latter approach has its advantages however: it is suitable to model the effects of the error field, while the former has not been generalized yet to do this. The general idea of the importance sampling method is the following. For a general function H and a CDF P, one can use the Monte-Carlo method to evaluate the integral, where ξ n are identically distributed independent random variates (RV), with the CDF P(x). The Monte-Carlo method can yield a poor convergence if some very rare realizations may contribute a large value to the sum. In such a case, however, the convergence can be improved using the importance sampling method, a general technique in which the sampling happens from a different distribution that generates the important addends more often. The integral (23) can be transformed to where Q(x) = P(ξ n < x) is the CDF of some distribution that generates more realizations in the interesting region. The last integral can be approximated numerically, At the beginning of the simulation, the electron's initial velocity is a random vector sampled from a non-isotropic 3D Gaussian distribution. The electrons are placed at distance h from the plane, which is perpendicular to the electron's initial velocity v and contains the ion (see figure 7).
For this particular distribution, we found that (12) can be efficiently calculated by parameterizing uniformly distributed ). Corresponding CDFs can be transformed as The uniform distribution of φ is left unchanged. This transformation gives many more close collisions with small impact parameters and therefore greatly improves the convergence of the Monte-Carlo method.

Undersampling of close collisions
We have observed that the friction force calculated by VORPAL simulations is often less than the one predicted by the analytic formula (12). We have found that this happens because in numerical simulations the strong collisions occur very infrequently, and a given simulation may not have a single one of them, yet a sufficient number of such collisions must be sampled to accurately estimate the average friction force. We call this effect undersampling of close collisions.
For simplicity, we assume that all electrons have the same velocity v and they start |v − v ion |T /2 away from the plane containing ions and perpendicular to v (see figure 7). Therefore their trajectories are 'symmetric'. The ion velocity change in the kth collision depends only on the impact parameter ρ k , which is considered to be a random number with CDF independent of k, The change in ion velocity due to one such interaction is Similarly to section 2, the coordinate system is chosen such that v is parallel to the unit vector e z .
Here, we view the friction resulting from N collisions as a sum of random numbers N k=1 K (ρ k ) and develop an algorithm that calculates the quantiles of this sum's distribution. This allows us to estimate the number of macro-particles sufficient to obtain accurate results in a firstprinciple particle simulation. Also, it suggests that, for the number of collisions typical for the electron cooling section, the change in the ion's velocity in one pass can deviate from its average significantly.
Similarly to the Pareto distribution, here we have to deal with a sum of variates that take extremely large values (compared to their average) with extremely low probabilities. Given such a distribution, unless the number of addends is extremely large, the CDF of the sum of such RVs cannot be approximated efficiently by the Central Limit Theorem, on which we rely in Monte-Carlo-type simulations including PIC or BCC.
The idea of this approach is borrowed from the literature about the Pareto distribution [20]. We use the fact that, for some small number m, the distribution of the sum of N − m smallest addends can be approximated well by a Gaussian distribution. The CDF of the whole sum is found by integrating the CDF of this Gaussian over the space of realizations of the m largest addends. In our case, the largest addends correspond to the collisions with the smallest impact parameters.
In the following, ρ k represents the impact parameter of the kth collision, and ρ (k) is the sequence of the same numbers sorted in ascending order. We calculate the CDF, Here, ρ (1) . . . ρ (k) , . . . represent the sequence of ρ 1 . . . ρ k , . . . in ascending order. The m-fold integral is calculated using the Monte-Carlo method. The integrand is approximated using the Central Limit Theorem, where α(r ) = E ((K (ρ 1 )|ρ 1 > r )) and σ (r ) = Var ((K (ρ 1 )|ρ 1 > r )) (32) are calculated using the classic approximation (11). The results from numerical calculations are shown in figure 8. These numbers correspond to the case of perfectly cold electrons and, hence, overestimate the friction force, which is acceptable for our discussion of the effects of rare collisions with small impact parameters. For N < 2 × 10 5 , all quantiles increase as the probability of a realization in the important region of small impact parameters is increasing; the 80% interval is very large and the algorithm does not seem to converge. Only after N becomes larger than 2 × 10 5 does the distribution start to behave like a Gaussian, and the quantiles start to converge. The precise calculation of friction force is possible only for 10 8 -10 9 realizations per macroparticle. For 5 × 10 4 collisions, which is typical for a real cooling section, the 80% confidence interval is so large that the friction force cannot be calculated in a single simulation with a realistic number of electrons.

Conclusions
We have reviewed the classic formula for the dynamical friction force on an ion moving through an electron distribution, clarifying the strong underlying assumptions. As explained in this 15 paper, for the short interaction times and electron densities to be found in proposed high-energy electron cooling systems, the standard approximations are not accurate.
We have presented a new algorithm to calculate the friction force very efficiently. The algorithm exploits the fact that collisions with different impact parameters have dramatically different contributions, and many do not have to be calculated at the same level of detail. This algorithm takes into account the finite time of the interaction and does not rely on any assumptions regarding the maximal or minimal impact parameters. The efficiency of the algorithm enables rapid estimates of the error in computation of the friction force.
We have provided the first rigorous mathematical treatment of the effects of rare, strong collisions (i.e. those with small impact parameter). Contrary to the physical intuition developed from years of using the standard formula for dynamical friction, which indicates that small impact parameters can be ignored, we have shown that the friction force on an ion during a single pass through an electron cooling section can significantly deviate from its average due to rare, strong collisions, as shown in figure 8. Hence, any analytic formula that calculates the average friction force provides a limited view of the full dynamics.
This also shows a conceptual limitation of single-pass simulations of an electron cooling section. In a certain parameter regime, in order to estimate the friction force, a number of singlepass simulations have to be averaged, or some other techniques of noise reduction must be used.
We found a generalization of the classical formula that takes into account the effect of the error field. Such a formula has been verified to be in good agreement with the high-precision massive numerical simulations based on a minimal set of assumptions.