Semi-analytical description of the modulator section of the coherent electron cooling

In the coherent electron cooling, the modern hadron beam cooling technique, each hadron receives an individual kick from the electric field of the amplified electron density perturbation created in the modulator by this hadron in a co-propagating electron beam. We developed a method for computing the dynamics of these density perturbations in an infinite electron plasma with any equilibrium velocity distribution -- a possible model for the modulator. We derived analytical expressions for the dynamics of the density perturbations in the Fourier-Laplace domain for a variety of 1D, 2D, and 3D equilibrium distributions of the electron beam. To obtain the space-time dynamics, we employed the fast Fourier transform (FFT) algorithm. We also found an analytical solution in the space-time domain for the 1D Cauchy equilibrium distribution, which serves as a benchmark for our general approach based on numerical evaluation of the integral transforms and as a fast alternative to the numerical computations. We tested the method for various distributions and initial conditions.


I. INTRODUCTION
A few of years ago, a novel hadron beam cooling technique capable to deal with the accelerators operating in the range of few TeVs, the Coherent electron Cooling (CeC), was proposed [1]. Currently, a test facility is under construction at Brookhaven National Laboratory. For the present status of the developments of the machine, we refer to [2]. The CeC is the modern realization of the stochastic electron cooling, wherein the electron beam serves as a pick-up and a kicker. It records the information about the hadron beam via the electron density perturbations resulting from the shielding of the hadrons. Then, these perturbations are amplified in the free electron laser (FEL) section, and then, in the kicker, every hadron experience the electric field produced by its own amplified perturbation receiving kicks. Before the kicker, in the dispersion section, each hadron is emplaced such that these kicks accelerate or decelerate it depending on its energy deviation, thereby reducing the energy spread of the hadron beam. To analyze the performance of the CeC, all the sections of the device must be studied in detail; in particular, the shielding of a hadron in an electron beam should be computed with high precision.
In this article, we offer a theoretical description of the modulator section of the coherent electron cooling, i.e., we address the problem of dynamical shielding of a charged particle in an electron beam. The simplest interpretation of this problem is a screening of a stationary particle in an infinite plasma (by plasma we * Electronic address: andrey.elizarov@stonybrook.edu † Electronic address: vl@bnl.gov mean a collisionless single-species electron plasma), i.e., the well-known Debye screening. The next step is to consider a moving ion in an infinite plasma. This problem was studied recently and the density perturbation for the Lorentz distribution was expressed as a one-dimensional integral [3]. The most advanced approach to resolving this problem is a general method for a shielding in a finite electron beam that takes into account the focusing field and the space-charge effects; this methodology was proposed last year by the authors of the present paper [4] and is under development now. There are also simulations of this effect using the PIC ("particle-in-cell") method [5]. In this article, we present a solution of this problem via the Fourier and Laplace transforms for the 1D, 2D, and 3D infinite plasmas and various equilibrium distributions. We derived the expressions for the solution in the Fourier-Laplace space, then inverted them numerically. For the 1D plasma with the Cauchy equilibrium distribution, we found a fully analytical solution in the space-time domain, which gives an opportunity to test the semi-analytical ones involving numerical evaluations of the inverse integral transforms and to perform many-particle computations much faster than with non-exact solutions. This method can also work with empirical equilibrium distributions. As a fast and robust solution, it has its own practical value and it will also serve as a testing ground for the PIC simulations and the general solution mentioned for a realistic case of a finite beam.

II. THE VLASOV-MAXWELL SYSTEM
Generally, the shielding of a charged particle in a plasma is described by the Vlasov-Maxwell system of equations [6], i.e., the dynamics of the electron den-sity is governed by the Vlasov equation and the electromagnetic field by the Maxwell equations. We first describe the system in a co-moving frame of reference, then derive a formal solution via the integral transforms, then introduce convenient dimensionless variables, and finally write a solution for a particle moving along a straight line.

A. General formulation for an infinite plasma
We consider the Vlasov-Maxwell system for the 1D, 2D and 3D plasmas simultaneously, which means that x is a one-, two-or three-dimensional vector depending on the dimensionality of the plasma we are considering and by x we denote its absolute value, even for the 1D case; the same conventions are applied for the dimensionless vectors that we will introduce in subsection II C. For the electron phase-space density f ( x, p, t), the Hamiltonian H, and the electric potential U ( x, t), we have: and the charge density is en( x). For this system, we consider the test charge problem with an external timedependent density d( x, t). We assume that is an equilibrium electron density, and f 1 = f 1 ( x, p, t) is an unknown perturbation resulting from the interaction with the test charge. The linearized Maxwell-Vlasov system looks as follows: B. Solving the Vlasov-Maxwell system via the integral transforms The Poisson equation (5) can be solved via the Fourier transform: whereŨ ( k, t),ñ 1 ( k, t), andd( k, t) are the Fourier images of the corresponding functions. Using this solution, we transform the equation (4) to [7]: where,Ñ 1 ( k, s) and LF d( x, t) are, respectively, the Laplace-Fourier images of n 1 ( x, t) and d( x, t): We introduce one more notation: Denoting the inverse Fourier and Laplace transforms, respectively, by F −1 and L −1 , we obtain the following expression: for the details on definitions of the integral transforms and our notations, see Appendix A. Generally, the expression (11) can be complex. Looking back to our initial equations and assuming complex f 1 , we note that the equation with Imf 1 corresponds to the equation without an external charge, while the equation with Ref 1 is the one with it, consequently, Imf 1 = 0 and it is confirmed by further computations. Hence, the expression (11) is real, as it should be. Even though the Poisson equations and their Green's functions differ for the 1D, 2D, and 3D cases, their solutions in the Fourier domain (6) and the expression (11) for n 1 ( x, t) have the same form.
In proceeding further, we need to specify the dimension of the problem, the external charge density d( x, t), and the equilibrium distribution, but first we introduce dimensionless variables.

C. Introducing dimensionless variables
We define the dimensionless variables, denoting them using the sans-serif font, as follows: where are, respectively, the root-mean-square velocity, the plasma frequency, and the Debye radius. The equilibrium density is normalized via: We introduce the dimensionless equilibrium density f 0 ( v) by the relation: wherein all the dimensional constants are gathered into f d and d stands for the dimensionality of the space, and can be 1, 2, or 3. We have the following dimensionalities for other quantities: We note that f d and v rms are not the same for the different equilibrium densities and must be computed via (13) and (16); for the non-integrable densities, the values have to be chosen voluntarily; among the densities we consider, only the Cauchy one is of that type. Using the dimensionless units, we rewrite formula (11) as follows: where LF d( x, t) and LF kt (tf 0 ( v)) are the dimensionless analogs of (9) and (10), respectively, 1 is a dimensionless factor, and L −1 , F −1 are the inverse Laplace and Fourier transforms for the dimensionless variables.

D. The external point charge
We assume that the charge's trajectory is unaffected by the space charge fields and consider the charge moving along a straight line y (t) = x 0 + v 0 t, we have: this assumption is reasonable for a hadron moving in an electron beam, as the hadron's mass is much larger than the electron's. For simplicity, we assume Z = 1 and the final density for the non-unitary charge can be recovered just by multiplying it by Z. Using the dimensionless units introduced, for any number of the spatial dimensions, we have: Finally, we can write the expression for the electron density perturbation resulting from the interaction with the external charge moving along a straight line y (t) = x 0 + v 0 t, valid in 1D, 2D, and 3D spaces: In the next section, we consider this solution for some particular equilibrium densities f 0 ( v); for each case, we just need to compute LF kt (tf 0 ( v)) and f −1 d v −d rms and insert them into (21).

III. APPLICATION TO THE PARTICULAR EQUILIBRIUM DISTRIBUTIONS
Generally, the equilibrium distribution f 0 ( v) has to be a solution of the unperturbed Vlasov equation, i.e., it has to be a function of the unperturbed Hamiltonian, in our dimensionless units it is v 2 , thus we consider the following functions: They correspond to the Kapchinskij-Vladimirskij (KV), water-bag (WB), normal (or Maxwell), and Cauchy (or Lorentz) equilibrium distributions, Θ(v) stands for the Heaviside step function, d is the dimensionality of the space, and v is a one-, two-, or three-dimensional vector. However, for the case of an infinite plasma that we are considering here, any function of velocity is a solution of the unperturbed Vlasov equation; thus, all our formulae can be easily generalized for the equilibrium distributions of the form: corresponding to an anisotropic plasma, where a i , i = 1, ..., d are dimensionless constants characterizing the plasma's temperatures. The changes should be applied only to the expression for LF kt (tf 0 ( v)), i.e., k i should be substituted with k i /a i for i = 1, ..., d, and the whole expression should be divided by For the 1D Cauchy distribution, the inverse Laplace and Fourier transforms in (21) can be evaluated analytically, while for the other distributions, the numerical techniques should be applied. We describe the 1D Cauchy case in detail and just quote the results for the other distributions starting with the KV and WB, which have different expressions for LF kt (tf 0 ( v)) in spaces of different dimensionalities. We conclude this section with the Cauchy and Maxwell distributions that have the same expressions for this quantity in all cases. Although, below we present v rms computed via (13), we considered dimensionless equilibrium distributions, f 0 ( v), corresponding to with this v rms , in all cases, , where β and H c are dimensional constants that can be used for fitting the experimental distributions. Computing LF kt (tf 0 ( v)) via (10), we obtain the following expression: or, using the dimensionless variables: Then, we insert the expression for LF kt (tf 0 ( v)) into the formula (21) and obtain: For all distributions we are considering, excepting the 1D Cauchy, the inverse integral transforms in the corresponding expressions for n 1 ( x, t) have to be inverted numerically, while, for the 1D Cauchy, they can be computed analytically giving the following expression: where and E 1 (z) and Ei(z) are the exponential integral functions [8] that can be computed via the series expansions, for details, see Appendix B 1. The whole expression (29) is real even though it contains complex numbers.
In Fig. 1, we show the densities obtained via the exact formula (29) and the ones obtained by the discussed in section IV numerical inversion of the integral transforms for x 0 = 0 and v 0 = 0.2, 1.0, 10.0. We note the perfect agreement of the exact solution and the numerical one. The solution has several interesting features, i.e., starting from some time, the left tail of the density has negative values, meaning that there is an accumulation of the charge of the same sign as that of the external perturbation, its maximum is oscillating, and the shape of the peak depends on the charge's velocity, being spiky for small velocities, widening as it increases, and, for large velocities, a discontinuity of the density's shape derivative appears in the right tail. All these features are equally well captured by the numerical computations and the analytical formula. We will comment further on the parameters' values in subsection IV B.

B. The 1D KV and WB distributions
For the 1D KV distribution, we have f 0 ( v) = ρH c β Hc δ βv 2 − H c . Using the dimensionless variables, we obtain: For the 1D WB, we have:

C. The 2D KV and WB distributions
In the 2D case, we have for the KV distribution: and, for the 2D WB: (41)

D. The 3D KV and WB distributions
The expressions are slightly bulkier in the 3D case; we obtain for the 3D KV: whereÎ and k 3 is a third component of k; for the 3D WB, we obtain: the integral in (49) has to be computed numerically.

E. The general Cauchy distribution
For the Cauchy distribution, we obtained the expressions valid in 1D, 2D, and 3D cases:

F. The normal distribution
For the Maxwell distribution, we also found universal formulas valid in 1D, 2D, and 3D cases: where Erfc(z) is the complementary error function [8], for its definition and some computational details, see Appendix B 2.

IV. NUMERICAL METHODS AND RESULTS
In this section, we briefly discuss the numerical methods we employed and present our results.

A. A few remarks on integral transforms inversion
To evaluate the expression (21) for a particular distribution, we first need to compute LF kt (tf 0 ( v)) via the formulas presented in the previous section; the expressions are either elementary functions or include special functions or a one-dimensional integral, all these things can be computed straightforwardly. Next step is an evaluation of the inverse Fourier and Laplace transforms. It is well-known that the inverse Fourier transform can be approximated by the discreet Fourier transform and then computed using the FFT algorithm, for details, see Appendix A. In this algorithm, the domain of interest of the resulting function is divided into N = 2 q segments. The inverse Laplace transform can be expressed via the inverse Fourier transform: which can be evaluated in a way that we described above, σ is a real constant greater than the real parts of all singularities off . In the upcoming section, we graphically present our results using the dimensionless units. We note that the dimensionless values for the different distributions are not always comparable to each other, since the values for v rms can be different; the corresponding conversion factors should be applied.

B. Results and discussion
In this section, we discuss the results obtained numerically and shown in Fig. 2, 3, 4, and 5; the velocity of the external charge, v 0 , is measured in units of v rms corresponding to the electron's equilibrium distribution and the initial position of the charge, x 0 , is measured in units of r D . The possible space-time ranges differ for the different distributions and are limited by the required precision and the number of points N = 2 q in the FFT algorithm. For each plot, we increased q until the values stabilized; the values used are shown in the legends in each plot. The most well-behaving case corresponds to the Cauchy distribution, the KV and normal distributions require greater values of q. For the 1D Cauchy distribution, the numerical results were already shown in Fig. 1.
In Fig. 2, we show the densities computed numerically for all 1D distributions for v 0 = 1. For the 1D KV distribution, we see that beam's response is a delta function-like peak, for the WB, the density is very spiked and asymmetric. For all distributions, excepting the KV, an accumulation of the charge of the same sign as the external charge occurs. For the normal equilibrium distribution, the perturbation is skewed and spiked resembling the shape of the α-stable distribution. For the Cauchy case, the perturbation is also skewed and spiked for v 0 ≈ 10 −1 , as illustrated in Fig. 1. While all other distributions exhibit a symmetric peak around the external charge for v 0 = 0, the KV distribution has two peaks that spread out with time, as shown in Fig. 3; with the increase of the velocity, the relative sizes of the peaks change, and, for v 0 = 1, the left peak almost disappears and the right one looks almost like delta-function, as evidenced in Fig. 2.
In Fig. 4, the lines of equal densities, for a certain set of times, for all 2D distributions considered, are shown for v 0 = 4. For the 2D KV distribution, we see the spreading out delta function-like "fronts", similar to the 1D case for v 0 = 0. For the 2D WB distribution, the lines are triangular with a peak following the charge. For the 2D Cauchy distribution, outer lines are almost circular; for the normal distribution, they have a bit more complicated shape. For the smaller velocities, the profiles are less directed toward the charge.
In Fig. 5, the lines of equal densities in certain planes, i.e., in three planes, each of which is parallel to the two out of the three coordinate axes, are illustrated. The shape of the lines for every distribution has the same features as the ones for 1D and 2D cases.
We emphasize that our 1D and 2D cases are not simple reductions of the 3D case, but correspond to the 1D and 2D theories of electrodynamics with the different from the 3D case Green's functions; the expressions for LF kt (tf 0 ( v)) also have different forms for the same distribution, but for the different dimensionalities. Although, in the method we consider in the present article, the Green's functions are not used and the solutions of the Poisson equation in the Fourier space were employed, which have the same form (6) in all cases, the general method for a finite beam [4] uses Green's function explicitly and its singularities is the main difficulty there. For the 1D case, the Green's function is not singular and the general method can be implemented without accounting for the singularities and the results can be compared with the exact solution (29) for the infinite beam; after a certain change of variables, the general method is also able to deal with the infinite beam case. This will be a reliable test of the method itself; and the 2D and 3D results of the present paper will help to develop and test good ways to handle the singularities. The general method is currently under development and we will elaborate on this in our future publications. The similarity of the 1D, 2D and 3D cases is also a sign of the applicability of the exact 1D results for estimating the shielding in the real device.

C. The code
The method we discussed herein was implemented as an object-oriented program in C++. The solution is stored as a multidimensional array over some grid in space-time; for further usage, it can be evaluated for any point using interpolation. The program is easily expendable for other external charge densities and equilibrium distributions and, in particular, can deal with the empirical ones. The visualization is also very flexible: it is possible to adjust time values, the number of equal density lines, set the particular values of interest, and look at different projections and cross-sections of the 2D and 3D densities.

D. Application to the Proof-of-principle experiment
As it was mentioned in the introduction, the proof-ofprinciple (PoP) experiment is planned at Brookhaven National Laboratory and the corresponding facility is cur-rently under construction [2]. In this subsection, we describe how the results can be applied to the modulator of the real device. To recover the dimensional quantities, we need the Debye radius and the plasma frequency; for the PoP, we have r D = 4.65 · 10 −5 m and ω p = 6.436 · 10 9 s −1 . We obtain for the dimensional density perturbation: where d is a spatial dimension of the problem, for the real 3D case, d = 3. For the PoP experiment, the modulator is constructed such that the interaction time is about one half of the plasma period, it depends on the hadron's velocity, as the modulator length is constant. The velocity is measured in units of v rms , in the PoP experiment, we have v rms = 3.0 · 10 5 m s . Our computations, shown in Fig. 2, demonstrate that extending the modulator up to a few plasma oscillations can significantly increase the density perturbation, i.e., its maximum will be up to four times greater. Further increasing of the modulation time doesn't increase the perturbation, as the modulator saturates, as shown in Fig. 2. The amplification of the perturbation in the FEL section is limited by the FEL saturation. For the model-independent description of the FEL saturation and its application to the theory of CeC, see [11] and [12], respectively. These considerations pro-vide limitations on a possible amplified perturbation that we can get, which, in their turn, determine the performance of the CeC device.

V. CONCLUSION AND FUTURE PLANS
In the present article, we considered a possible way to model the modulator section of the coherent electron cooling, i.e., we developed a method for evaluating the dynamical shielding of an external charge in an infinite electron plasma; for the certain case, we found analytical solution. The software package we developed gives reliable results for a variety of equilibrium distributions and initial conditions. We plan to use it in the analysis of the next section of the CeC, the FEL section, wherein the electron density perturbation from the modulator evolves in a free electron laser [9,10]. The results obtained can also be used as a testing ground for the more general method for a finite beam and for the PIC simulations.
Recently, we proposed a full theoretical model of the CeC [10]. All sections, i.e., the modulator, the FEL amplifier, and the kicker, were described using the inverse integral transforms. The kicker can be described in a very similar way to the modulator. For the details on the FEL section, we refer to [9,10]. The methods for the inverse integral transforms inversion that we developed and tested in the present article open an opportunity to implement this model and get a reliable and fast complete numerical model of the CeC.

VI. ACKNOWLEDGMENTS
Various communications with G. Wang, S. Webb, I. Pogorelov, and A. Fedotov are gratefully acknowledged. We thank A. Woodhead for proofreading.