A stochastic-hydrodynamic model of halo formation in charged particle beams

The formation of the beam halo in charged particle accelerators is studied in the framework of a stochastic-hydrodynamic model for the collective motion of the particle beam. In such a stochastic-hydrodynamic theory the density and the phase of the charged beam obey a set of coupled nonlinear hydrodynamic equations with explicit time-reversal invariance. This leads to a linearized theory that describes the collective dynamics of the beam in terms of a classical Schr\"odinger equation. Taking into account space-charge effects, we derive a set of coupled nonlinear hydrodynamic equations. These equations define a collective dynamics of self-interacting systems much in the same spirit as in the Gross-Pitaevskii and Landau-Ginzburg theories of the collective dynamics for interacting quantum many-body systems. Self-consistent solutions of the dynamical equations lead to quasi-stationary beam configurations with enhanced transverse dispersion and transverse emittance growth. In the limit of a frozen space-charge core it is then possible to determine and study the properties of stationary, stable core-plus-halo beam distributions. In this scheme the possible reproduction of the halo after its elimination is a consequence of the stationarity of the transverse distribution which plays the role of an attractor for every other distribution.


I. INTRODUCTION
High intensity beams of charged particles, in particular in linacs, have been proposed in recent years for a wide variety of accelerator-related applications: drivers for sources of neutron spallation; production of tritium; transmutation of radioactive wastes to species with shorter lifetimes; heavy ion drivers for fusion-based production of thermonuclear energy; production of radioactive isotopes for medical use, and so on. In all these cases it is very important to keep at a low level the beam losses to the wall of the beam pipe, since even small fractional losses in a high-current machine can cause exceedingly high levels of radioactivation. It is now widely believed that one of the relevant mechanisms for these losses is the formation of a low intensity halo relatively far from the core of the beam. These halos have been either directly observed [1] or inferred from experiments [2], and have also been predicted from extensive numerical simu- * cufaro@ba.infn.it † demartino@sa.infn.it ‡ desiena@sa.infn.it § illuminati@sa.infn.it lations [3,4,5,6,7,8,9,10,11]. Often in these studies the particle-core model has been used to understand the halo formation and its extent. That notwithstanding, it is widely believed that for the next generation of high intensity machines it is necessary to obtain a more quantitative understanding not only of the physics of the halo, but also of the beam transverse distribution in general [12,13,14]. In fact "because there is not a consensus about its definition, halo remains an imprecise term" [15] so that several proposals have been put forward for its description. Charged particle beams are usually described in terms of classical deterministic dynamical systems. The standard model is that of a collisionless plasma where the corresponding dynamics is embodied in a suitable phase space (see for example [16]). Generally speaking, the beams are described in a comoving reference system so that we can confine ourselves to a non relativistic setting. In this framework, the formation of the halo has been studied mainly by means of the particle-core model initially proposed by Gluckstern, and the simulations show that the dynamical instabilities due to a parametric resonance can allow ions to escape from the core [3,4,5,6,7,8,9,10,11].
In the present paper we propose and develop a different approach. We introduce a model for the formation of the halo based on the idea that the trajectories of the charged particles are sample paths of a stochastic process, rather than the usual deterministic (differentiable) trajectories. A random collisional plasma described in terms of configurational stochastic processes seems a realistic model for the dynamics taking place in particle accelerators. Moreover, it can explain the rare escape of particles from a quasi-stable beam core, since it takes into account, statistically, the complicated inter-particle interactions that cannot be described in detail. Of course, the idea of a stochastic approach is hardly new [16,17,18], but there are several different ways to implement it. The system we want to describe is the particle beam endowed with some measure of time-reversal invariance. For such a physical situation one is faced with a many-body dynamical system which is at the same time stochastic (due to the random collisions and the fast, irregular configurational dynamics) and conservative (time-reversal invariant). One then needs a generalization of classical mechanics for conservative deterministic systems to a classical mechanics of conservative stochastic systems, i.e. a stochastic mechanics.
Despite the rather widespread misconception that stochastic processes in dynamics are always associated to the description of dissipative and irreversible behaviors, it is in fact possible for specific stochastic dynamical systems to exhibit time-reversal invariance. In general, these systems are such that a deterministic dynamics in terms of known external potentials is ascribed to a system whose kinematics, due to some intrinsic irregularities, is instead stochastic [19]. This is the case we conside for the collisional beam.
In previous papers [20,21] we analysed some basic properties of stochastic mechanics (SM) and we also introduced some basic criteria of mechanical stability in order to provide phenomenological support to the scheme of SM for classical stochastic dynamical systems. In particular, in Ref. [20] we applied SM to the description of the dynamics of charged particle beams, and we showed how the criteria of mechanical stability allow to connect the (transverse) emittance to the characteristic microscopic scales of the system and to the total number of particles in a given beam. Moreover, this method allows to implement techniques of active control for the dynamics of the beam. In Ref. [20] these techniques have been proposed to improve the beam focusing and to independently change the frequency of the betatron oscillations.
In the present paper, the SM of particle beams is used to investigate the nature, the size and the dynamical characteristics of beam halos. In particular, we determine the effects of the space-charge interaction on the formation of the halo by calculating the growth of the beam emittance, and by estimating the probability to find particles far away from the beam core. We then consider the reverse problem i.e., by imposing reasonable forms of the stationary halo distributions we derive the effective spatial potentials that are associated to such distributions.
The paper is organized as follows: in Sections II and III we summarize the results of previous papers on the model of SM for particle beams with special emphasis on the techniques of active control needed to engineer a desired beam dynamics. We then apply SM to study the formation of the beam halo. In Section IV we show how to treat with SM the effects of space charge due to the electromagnetic interactions among the particles. We introduce and solve self-consistently a set of coupled, nonlinear hydrodynamic equations for the beam density and phase, the space-charge potentials, and the external electromagnetic field. Solution of the coupled equations leads to the first main result of our paper, i.e. that the final transverse particle distribution is greatly broadened due to the presence of the space-charge potential. This broadening gives rise to a penetration of the distribution well into spatial regions very far from the core of the beam, thus realizing a type of continuous halo, i.e. without intermediate voids of particles (nodes of the transverse beam distribution). This effect is reflected in a corresponding growth of the transverse emittance that we estimate from the solution of the self-consistent equations. In Section V we study the reverse problem: given different stationary transverse distributions of the halo around the core of the beam, such as a continuous broadened distribution or a distribution with a gap (void) between the halo and the core, we determine the dynamics that leads to such configurations. This is important since it enables us to deduce the form of the effective potentials which control the stationarity of the halo distribution. Finally, in Section VI we draw our conclusions and discuss future directions of research.

II. STOCHASTIC BEAM DYNAMICS
The usual way to introduce a stochastic dynamical system is to modify the dynamics in phase space by adding a Wiener noise only in the equation for the momentum, in order to preserve the usual relations between position and velocity: where m is the mass of the system, Q(t) is the position variable, P(t) is the momentum variable, F is the external force, and B(t) is a fluctuating force modelled by a Wiener process; in general, this procedure yields a differentiable, but non Markovian process Q(t) for the position (here β denotes the autocorrelation of the fluctuating force). The standard example of this approach is that of a Brownian motion in a force field described by the Ornstein-Uhlenbeck system of stochastic differential equations (SDE) [22]: where K is the spatial gradient of the external potential, and γ is the damping constant. This system can be reduced to yield a non differentiable, Markovian process Q(t) acted on by a Wiener noise, with diffusion coefficient D, to the equation for the position. In this case we are forced to drop the equation for the momentum since Q(t) is no longer differentiable, and we can reduce the stochastic system to the single SDE where v (+) (r, t) is the so called forward velocity, and dW(t) ≡ W(t+dt)−W(t) is the increment of a standard Wiener noise W(t); as it is well known, this increment is a Gaussian process. The standard example of this reduction is the Einstein-Smoluchowski approximation of the Ornstein-Uhlenbeck process in the overdamped case. It is however possible to introduce Eq. (3) in a completely different context, as the defining equation for a theory of stochastic dynamical systems in configuration space, rather than in phase space. This amounts to consider a system whose kinematics, rather than its dynamics, is taken ab initio to be random. By doing this, no external sources of dissipation and irreversibility are introduced on the forces acting on the system. Instead, a source of randomness is assumed that perturbs only the configurations of the system. For such systems with random kinematics, the dynamics can then be introduced either by generalizing the Newton equation [22,23], or by introducing a stochastic variational principle [22,24]. The crucial aspect is that such a dynamics turns out to be both stochastic and conservative. In this scheme, since Q(t) is not differentiable, there is no velocity as a standard derivative. Instead it is possible to define in a suitable way [22] an average velocity v (+) in the forward time direction, and an average velocity in the backward time direction v (−) : these functions of r and t are different, and they both coincide with the usual velocity only if the process is differentiable, i.e. if the kinematics is deterministic. The relations between v (+) and v (−) will be introduced in the subsequent discussion. It is important to remark that v (+) (r, t) is no more an a priori given field: it plays now the role of a dynamical variable.
The stochastic dynamical scheme sketched above is known as Stochastic Mechanics (SM), and most of its applications have concerned the problem of developing a classical stochastic model for the simulation of Quantum Mechanics. Nonetheless, it is a very general model which can be applied to very different stochastic dynamical systems endowed with time-reversal invariance [25]. We will show below how one can derive from the stochastic variational principle two coupled, nonlinear hydrodynamic equations for the density and the phase of a dynamical system in configuration space. These two real, coupled equations can be recast into a single complex equation whose form is completely equivalent to that of a Schrödinger differential equation. In this sense, some authors denote classical dynamical systems described by SM as quantum-like systems, in analogy with other re-cent studies on the collective dynamics of charged particl beams [26,27].
In general, SM can be used to describe every conservative, stochastic dynamical system satisfying fairly general conditions: for instance, it is known that for any given conservative diffusion there is a correspondence between the associated diffusion process and a solution of the Schrödinger equation with Hamiltonians associated to suitable vector and scalar potentials [28]. Under some regularity conditions this correspondence is one-to-one. The usual Schrödinger equation, and hence quantum mechanics, are recovered when the diffusion coefficient D is independent from the values taken by the conservative diffusion process and coincides with /2m, where denotes the Planck constant. Here we are not interested in the stochastic modelling of quantum mechanics, but rather in the classical stochastic dynamics of particle beams. In this instance the unit of action appearing in the effective Schrödinger equation is linked to the beam emittance as will be clarified below.
We will now describe the stochastic process performed by a representative particle of the beam that oscillates, in a comoving reference frame, around the closed ideal orbit in a particle accelerator. We consider the threedimensional (3d) diffusion process Q(t), taking the values r, which describes the motion of the representative particle and whose probability density coincides with the particle density of the beam. The evolution of this process is ruled by the Itô SDE (3) where the diffusion coefficient D is supposed to be constant. The quantity α = 2mD, which has the dimensions of an action, will be connected later to the characteristic transverse emittance of the beam. Eq. (3) defines the random kinematics performed by the collective degree of freedom, and replaces the simple deterministic kinematics for the differentiable trajectory q(t).
We are in a situation in which we have both a random, diffusive kinematics and a time reversal invariance. Therefore, to counteract the dissipation, one must impose a conservative dynamics on the stochastic kinematics, at variance with the purely dissipative Fokker-Planck or Langevin dynamics. A conservative dynamics imposed on a random kinematics can be introduced by a suitable stochastic generalization of the least action principle of classical mechanics [24]. This is achieved by replacing the classical deterministic kinematics (4) with the random diffusive kinematics (3) as the independent configurational variables of the Lagrangian density entering the action functional. The equations of motion thus obtained by minimizing the action functional take the form of two coupled hydrodynamic equations that describe the dynamical evolution of the beam profile and of the associated velocity field. In the following we briefly sketch the derivation of the hydrodynamic equations, referring for details to references [22,24].
Given the SDE (3), we consider the probability density function (pdf) ρ(r, t) associated to the diffusion Q(t) so that, besides the forward velocity v (+) (r, t), we can also define the backward velocity It is also useful to introduce the current and the osmotic velocity fields, defined as: The velocities in Eq. (6) have a transparent physical meaning: the current velocity v represents the global velocity of the density profile, being the stochastic generalization of the velocity field of a classical perfect fluid.
On the other hand the osmotic velocity u is clearly of intrinsic stochastic nature, for it is a measure of the non differentiability of the stochastic trajectories, and it is related to the spatial variations of the density. In the limit of a deterministic process, i.e. of a diffusion that tends to a deterministic, differentiable trajectory q(t), the current velocity v(r, t) tends to the classical velocity field v(q(t), t), and the osmotic velocity u tends to zero (the trajectory becomes differentiable, so that the forward and backward velocities, i.e. the left and right derivatives coincide). In order to establish the stochastic generalization of the least action principle one introduces a mean classical action averaged over the probability density function of the diffusion process, in strict analogy to the classical deterministic action. In fact, the main difficulty in the stochastic case is due to the non differentiable character of the sample paths of a diffusion process which does not allow to define the time derivative of the process Q(t). Hence the definition of a Lagrangian density and of an action functional is possible only in an average sense through a suitable limit on expectations. The stochastic action is then defined as [22,24] where E( · ) = ( · )ρ(r, t) dr denotes the expectation of a function of the diffusion process Q(t), V is an external potential, and ∆Q(t) = Q(t + ∆t) − Q(t) is the finite increment of the process. It can be shown that the mean action (7) associated to the diffusive kinematics (3) can be recast in the following particularly appealing Eulerian hydrodynamic form [22]: where v and u are defined in Eq. (6). This Eulerian form of the action immediately shows that when u = 0 (and hence the trajectories are differentiable) it coincides with the action for the conservative dynamics of a classical Liouville fluid. The stochastic variational principle now follows by imposing the stationarity of the stochastic action (δA = 0) under smooth and independent variations δρ of the density, and δv of the current velocity, with vanishing boundary conditions at the initial and final times. As a first consequence we get that the current velocity has the following gradient form: which can be also taken as the definition of the phase S. The non linearly coupled Lagrange equations of motion for the density ρ and for the current velocity v, of the form (9), are a continuity equation typically associated to every diffusion process and an evolution equation for the conservative dynamics The above equations give a complete characterization of time-reversal invariant diffusion processes. The last equation is formally analogous to the Hamilton-Jacobi-Madelung (HJM) equation originally introduced in the hydrodynamic formulation of quantum mechanics by Madelung [29]. However, the physical origins of the two equations are profoundly different, as in the quantum mechanical case the diffusion coefficient is related to the fundamental Planck constant by the relation D = /2m. Notice that when D = 0 (namely when the kinematics is not diffusive) the equation (11) coincides with the usual Hamilton-Jacobi equation for a Liouville fluid. Since (9) holds, the two equations (10) and (11) can be recast in the following form which now constitutes a coupled, non linear system of partial differential equations for the couple (ρ, S) which completely determines the state of our beam. On the other hand, because of (9), this state is equivalently given by the couple (ρ, v). Eqns. (12) and (13) describe the collective behavior of the beam at each instant of time through the evolution of both its particle density and its velocity field. It is important to notice that, introducing the representation [29] ψ(r, t) = ρ(r, t) e iS(r,t)/α , with the two coupled real equations (12) and (13) are equivalent to a single complex Schrödinger equation for the function ψ, with the Planck action constant replaced by the unit of action α: In this formulation the phenomenological "wave function" ψ carries the information on the dynamics of both the beam density and the beam velocity field, since the velocity field is determined via Eq. (9) by the phase function S(r, t). This shows, as previously claimed, that our procedure, starting from a different point of view, leads to a description formally analogous to that of the socalled quantum-like approaches to beam dynamics [26]. We remark again that obviously Eq. (16) has not the same meaning as in quantum mechanics. The parameter α defined in Eq. (15) has the dimensions of an action, but in general α = . In fact, α is not a universal constant and its value cannot coincide with , as its value depends on the physical system under study. However, α plays in some sense a role similar to that of in quantum mechanics since, as we will see in the following, it identifies a lower bound for the beam emittance in phase space, and gives a measure of the position-momentum uncertainty arising from the stochastic dynamics of the beam. Thus the Schrödinger equation (16) for the stochastic mechanics of particle beams presents some features reminiscent of quantum mechanics, but at the same time is a deeply different theory with a different physical meaning. In particular, while in quantum mechanics a system of N particles is described by a wave function in a 3Ndimensional configuration space, in the scheme of SM the normalized probability density distribution |ψ(r)| 2 is a function of only the three space coordinates in physical space. If there are N particles in the beam, then the function N |ψ(r)| 2 is simply the particle density in physical space.

III. CONTROLLED BEAM STATES
In the previous section we have introduced two coupled equations that describe the dynamical behavior of the beam: the first is the Itô equation (3), or equivalently the continuity equation (either (10) or (12)); the second is the Hamilton-Jacobi-Madelung (HJM) equation (either (11) or (13)). Here we briefly summarize a general procedure, already exploited in Refs. [20,30,31], to engineer a controlled dynamics of stochastic systems. We will then apply this method to the control of the beam halo in Section V.
To this end it is useful to show, by simple substitution from (6), that Eq. (10) is equivalent to the Fokker-Planck (FP) equation formally associated to the Itô equation (3). The HJM equation (11) can also be cast in a form based on v (+) rather than on v, namely: where we have defined the function f ≡ ρ/N , and N is a constant such that f is dimensionless. On the other hand, we know from Eqns. (5) and (6) that also the forward velocity v (+) is irrotational: where the scalar functions W and S are connected by the relation and θ is an arbitrary time-dependent function. Through Eq. (19), Eqns. (17) and (18) can be cast in a form analogous to the system of Eqns. (12) and (13). In this case, the couple of functions which determine the state of the beam is (ρ, W ) (or equivalently (ρ, v (+) )). It is worth noticing that the possibility of a time reversal invariance [23] is assured by the fact that the forward velocity v (+) (r, t) is not an a priori given field, as is the case for dissipative diffusion processes of the Langevin type. Rather, it is dynamically determined at any instant of time, for any assigned initial condition, through the HJM evolution equation (11).
Let us suppose now that the beam density ρ(r, t) is given at some time t. We may think for instance of an engineered evolution from some initial density toward a final, required state with suitable properties. In particular, we could imagine to have a halo-ridden particle distribution in the beam that should be steered toward a final halo-free distribution. We know that ρ must be a solution of the FP equation (17) for some forward velocity field v (+) (r, t), which we consider here as not given a priori. Since also Eq. (19) must be satisfied, this means that -for a given ρ -we should find an irrotational v (+) in such a way that the FP equation (17) is satisfied. In other words, we are required to solve the partial differential equation for the irrotational vector field v (+) , or equivalently, by taking Eq. (19) into account, the second order partial differential equation for the scalar field W : This is not, in general, an easy task, but we will give here the solutions for a few simple cases that will be useful in the following.
In the one-dimensional (1d) case (ρ(x, t)), the FP equation defined for a < x < b (here we could possibly have a = −∞, and b = +∞) is: where J = D∂ x ρ − v (+) ρ is the probability current. It is easy to see now that the solution is where c(t) is an arbitrary function. Moreover, since the conservation of probability imposes, in particular, that J(a) = 0, it is easy to see that we must choose c(t) = 0 so that finally we have Of course, in this case we also have This 1d solution will be useful later for a simplified study of beam distributions when we take x as one of the transverse coordinates of the beam. Another case of practical interest is that of a 3d system with cylindrical symmetry around the z-axis (which will be supposed to be the beam longitudinal axis). If we denote with (r, ϕ, z) the cylindrical coordinates, where r = x 2 + y 2 , we will suppose that ρ(r, t) depends only on r and t, and that v (+) = v (+) (r, t)r is radially directed, with modulus depending only on r and t. In this case, it is straightforward to see that the equation for v (+) is In this case Finally, in the 3d stationary case the beam density ρ is independent of t. This greatly simplifies the treatment, as it is easy to check by simple substitution, and we simply find with This result corresponds to the well known fact that for stationary forward velocities the FP equation (17) always admits stationary solutions of the form Remark however that also for particular choices of a non stationary ρ(r, t) it is sometimes possible to find a stationary velocity field v (+) (r) such that the FP equation (17) is satisfied (think for instance to the Ornstein-Uhlenbeck non stationary solutions): this fact will be of use in the following Sections. Now that ρ and W (or equivalently v (+) ) are given and satisfy Eq. (17), we should also remember that they will qualify as a good description of our system only if they also are a solution of the dynamical problem, namely if they satisfy the HJM equation (18). Since W , v (+) and ρ are now fixed, this last equation must be considered as a constraint relation defining an external controlling potential V which, after straightforward calculations and in terms of the dimensionless distribution f , turns out to be of the general form: This procedure can in principle be applied also to more complicated instances, for example to engineer a beam dynamics that keeps the beam coherent even in the presence of aberrations, or to engineer some desired evolution. However the solutions will not always be available in closed form, and will in general require some approximate treatment. Here we present the simple analytic solutions for the three particular cases discussed above.
In the 1d case f (x, t), with v (+) given by Eq. (25), the controlling potential is: In the 3d, cylindrically symmetric case f (r, t), with v (+) given by (28), we have instead: Finally, in the 3d stationary case ρ(r) with v (+) given by Eq. (30), the potential is still time-dependent because of the presence of the arbitrary function θ(t), and it reads: In order to make it stationary it will be enough to require thatθ(t) be a constant E, so that In this stationary case the constant E fixes the zero of the potential energy, and the phenomenological wave function (14) takes the form typical of stationary states.

IV. SELF-CONSISTENT EQUATIONS
One of the possible mechanisms for the formation of the halo in particle beams is that due to the unavoidable presence of space charge effects. In this Section we will investigate this possibility in the framework of our hydrodynamic-stochastic model of beam dynamics. To this end, we take into account the space charge effects by coupling the hydrodynamic equations of stochastic mechanics with the Maxwell equations which describe the mutual electromagnetic interactions between the particles of the beam. We thus obtain a selfconsistent, stochastic magnetohydrodynamic system of coupled nonlinear differential equations that can be numerically solved to show the effect of the space charge.
In the following, the reference physical system will be an ensemble of N identical copies of a single charged particle embedded in a particle beam and subject to both an external and a space-charge potential. In a reference frame comoving with the beam, our system is then described by the Schrödinger equation (43), where α is the unit of action (emittance) and H the Hamiltonian operator which will be explicitly determined in the following. Since in general ψ is not normalized, we introduce the following notation for its constant norm so that, if N is the number of particles with individual charge q 0 , the space charge density of the beam will be On the other hand its electrical current density will be which vanishes when the wave function is stationary, namely when ψ(r, t) = u(r) e −iEt/α .
For a charged particle in the beam the electromagnetic field is the superposition of the space-charge potential (A sc , Φ sc ) due to the presence of ρ sc and j sc , plus the external potential (A ext , Φ ext ), and hence the Schrödinger equation takes the form Eq. (43) has then to be coupled with the Maxwell equations for the vector and scalar potentials (A sc , Φ sc ).
Since we are in a reference frame comoving with the beam, we can always assume that the wave function is of the stationary form (42). In this case we have j sc = 0, so that A sc = 0. Finally, by taking also A ext = 0, our system is reduced to two coupled, nonlinear equations for the pair (u, Φ sc ). For details see Appendix A. If the beam with space-charge interactions stays cylindrically symmetric the function u will depend only on the cylindrical radius r and our system becomes (see equations (A5) and (A6) in Appendix A): where with L the length of the beam, and Eq. (44) is the stationary Schrödinger equation in the external scalar potential, and Eq. (45) is the Poisson equation for the space-charge scalar potential. We proceed to solve numerically this nonlinear system and compare its solution with that for a purely external potential V ext which is a cylindrically symmetric, harmonic potential with a proper frequency ω in absence of space charge (see Appendix B): V ext = (mω 2 r 2 )/2, with r = x 2 + y 2 . Let us first introduce the dimensionless quantities where σ 2 = α/2mω is the variance of the ground state of the cylindrical harmonic oscillator without space charge.
Eqns. (44) and (45) can then be recast in the dimensionless form where The effect of the space charge will be finally accounted for by comparing the normalized solution w(s) with the unperturbed ground state of the cylindrical harmonic oscillator (see Appendix B): We have solved numerically the system (49), (50) by tentatively fixing one of the two free parameters b and β, and then searching by an iterative trial and error method a value of the other such that the solution shows, in a given interval of values of s, the correct infinitesimal asymptotic behavior for large values of s. We have then normalized the solutions w(s) by calculating numerically the value of B. It is clear from the definition of the dimensionless parameters B, b and β (51) that the value of β is a sort of reduced energy eigenvalue of the system, while the product which is by definition a non negative number, will play the role of the interaction strength, since it depends on the space-charge density along the linear extension of the beam. Reverting to dimensional quantities, since α 2 /2mσ 2 has the dimensions of an energy, the two relevant parameters are which are respectively the energy of the individual particle embedded in the beam, and the strength of the spacecharge interaction. In the following we will limit ourselves to discuss solutions of Eqns. (49), and (50) without nodes (a sort of ground state for the system). It is possible to see that no solution without nodes can be found for values of the space-charge strength γ beyond about 22.5 and that for values of γ ranging from 0 to 22.5, the energy β decreases monotonically from 2.0000 to −0.0894. If the unit of action α is fixed at a given value, it is apparent that the value of γ is directly proportional to the charge per unit length N q 0 /L of the beam: a small value of γ means a rarefied beam; a large value of γ indicates that the beam is intense. A halo is supposed to be present in intense beams, while in rarefied beams the behavior of every single particle tends to be affected only by the external harmonic potential. It is useful to provide a numerical estimate (in M KS units) of the relevant parameters. We are assuming a beam made of protons, so that m and q 0 are the proton mass and charge, while ǫ 0 is the vacuum permittivity. The parameter σ is determined by the strength of the external potential, and it is a measure of the transverse size of the beam when no space charge is taken into account, so that the ground state has the unperturbed form (52). From empirical data, a reasonable estimate yields σ ≈ 10 −3 m. On the other hand the value of N/L clearly depends on the particular beam we are considering: usually a value of 10 11 particles per meter is considered realistic. As for the parameter α (15) we have already suggested in Section I that its value can be connected to the beam emittance. On the ground of accepted experimental values of the beam emittance (usually measured in units of length) we can assume α/mc of the order of 10 −7 m. In fact, taking m, q 0 , σ, and N/L fixed at the above-mentioned values, and assuming a space-charge strength γ moderately large, i.e. 2.0 ≤ γ ≤ 22.5, we can determine α from (54) and we get 7.2 × 10 −7 ≤ α/mc ≤ 1.3 × 10 −7 , a number which is in good agreement with the experimentally measured values of the beam emittance. From now on we will fix α at its approximate central value: Since m and q 0 have the proton values, and σ ≈ 10 −3 m, then the assigned values of the dimensionless parameters γ and β will fix respectively, the values of N/L and E by  Eqns. (54). This implies that by changing the beam intensity (namely N/L and γ) one correspondingly changes the energy (namely E and β) of the individual particle embedded in the beam. In Fig. 1 we show the behavior of β as a function of γ. The quantities β and γ are dimensionless: the true physical quantities (energies) are obtained by multiplicating them by the unit of energy The overall effect of the space charge in this model is a conspicuous spreading of the transverse distribution of the particles in the beam with respect to the unperturbed ground state distribution (52). When γ = 0 the potential due to space charge vanishes, and the solution exactly coincides with the ground state of the cylindrical harmonic oscillator with variance σ 2 . When γ > 0 the transverse distribution begins to spread, as we show in    is the probability of finding a particle at a distance greater than cσ for systems with a given strength γ of the space charge coupling. For instance, considering the two different situations γ = 0 (no space charge) and γ = 22.5 (strong space charge) we have: We see that the probability of finding particles at a distance larger than 10σ from the core of the beam is enhanced by space charge by many orders of magnitude. This means, for example, that if in the beam there are 10 11 particles per meter, while practically no one is found beyond 10σ in absence of space charge, for very strong space-charge intensity we can find up to 10 5 particles per meter at that distance from the core.
The above analysis shows that in the hydrodynamicstochastic theory of charged beams, the space-charge potential induces a strong broadening of the unperturbed transverse density distribution of the beam, thus yielding a small, but finite probability of having particles at a distance well away from the core of the beam.
The hydrodynamic-stochastic model allows also an estimate on the growth of the emittance due to the presence of the space charge. The emittance can be calculated by exploiting a structure of uncertainty products that is inherent to the SM. In particular, the transverse emittance can be calculated as ∆x · ∆p x , where we have: More explicitly: if w(s) is a solution of the coupled dynamical equations normalized in the sense that the corresponding 3d probability density is where H(s) is the Heaviside function and This probability density function is normalized since It is easy to show that the first two moments are: so that the root mean square deviation for the position is: As for the momentum p, we recall that in SM one cannot introduce a probability distribution directly in phase space. The momentum can nevertheless be recovered from the velocity field which is well defined in SM. Since we are considering a stationary state, only the osmotic part of the velocity field is non zero. Then the momentum field reads As a consequence, the x-component of the momentum is and the first and second moments read Hence the root mean square deviation for the momentum is: We can then define the emittance E as the positionmomentum uncertainty product in phase space in the following way: (70) The emittance (70) can be thus estimated from knowledge of the numerical solutions w(s) and their first derivatives w ′ (s). It can be seen that the emittance is exactly α for γ = 0 (namely in absence of space charge), and grows with γ to a value ≈ 1.2 × α for γ ≈ 22.5. This result is consistent with the expected growth of emittance produced by space-charge effects. The dimensionless ratio ε = E/α = ∆x · ∆p x /α is plotted as a function of γ in Fig. 8. This figure and the above discussion provide evidence that in the model of SM of charged beams the constant α plays the role of a lower bound for the phase-space emittance.

V. STATIONARY HALO DISTRIBUTIONS
The coupled nonlinear equations introduced in the previous Section allow to introduce and explain a possible mechanism of halo formation due to space-charge effects. This mechanism, when considering solutions without nodes, amounts to the broadening of the transverse density distribution of the beam away from the beam core. This is a type of "continuous" halo. If we want to investigate the formation and the properties of other possible forms of halos -for example ring-like halos separated from the core by an almost void space region -we should try to analyze the solutions with nodes of our radial self-consistent equations. Being this a time consuming task, we resort for the time being to a different, simplified approach to study the structure and properties of beams of the ringed type.
Moreover, the scheme introduced in the previous Section, although in principle very general, has the drawback that it does not allow for analytical expressions of velocities and potentials, and therefore cannot suggest the engineering of controlling external potentials to remove the halo. To this end, in this Section we will develop a different approach, still based on SM, by considering some preassigned density distributions for a beam with halo, and then looking for the kind of dynamics that produces them.
Firs of all we will determine the analytical forms w(s) of the dimensionless radial distribution that can best approximate the numerical, self-consistent solutions w(s) found in the previous Section. We will do that by choosing a family of trial functions w(s) and by subsequently minimizing the mean square error (m.s.e.) with respect to a given solution w(s). Since the radial component of the cylindrically symmetric (l = 0) solutions of the harmonic oscillator (B1) contain only even powers of the radial coordinate (see Appendix B), we take as normalized trial functions with σ, a and b as free parameters, and c = ∞ 0 e −s 2 /σ 2 (1 + as 2 + bs 4 ) 2 s ds = σ 2 2 + aσ 4 + (a 2 + 2b)σ 6 + 6abσ 8 + 12b 2 σ 10 (72) as normalization constant. We then take one numerical solution w(s) of (49) and (50) for a given value of γ and consider the mean square error and numerically minimize its value by varying the parameters σ, a and b. As an example, for γ ≈ 19.3, we get The corresponding density distribution (solid line) is compared in Fig. 9 with the density distribution obtained by numerically solving the system (49) and (50) (dashed line). From this explicit form of the approximating w(s) we can now also calculate the expression of the stationary potential which controls the state. From Eq. (37) in the case of a cylindrically symmetric approximating state with we have the controlling potential Introducing the usual dimensionless quantities we may define the dimensionless controlling potential  Fig. 9 because of the large scale. We show a zoom-up of it in Fig. 11, where the node is now clearly visible. This is a feature that we will find as well in other different approaches to the description of the beam distribution, and that could be tentatively connected with the presence of a halo spatially separated from the core of the beam. In our quest for an analytic description of the halo distribution we now depart from the scheme of the selfconsistent nonlinear equations (49) and (50). Instead, while still assuming the space charge as the dominant factor, we introduce a simplified analytical model for the halo production. As a first simplification, valid for small values of the space-charge strength γ, we assume a frozen core for the beam. A frozen core is such that it produces a space-charge effect, induces a spreading of the density distribution, and is left unchanged by this spreading. In this case the potential produced by the charge distribution is a priori given by the distribution of the frozen core and it simply enters linearly in the phenomenological Schrödinger equation (16). As a consequence, the latter is now decoupled from the Maxwell equations for the electromagnetic field. The frozen charge distribution will be assumed cylindrically symmetric and transver-sally Gaussian: where N is the number of particles with elementary charge q 0 , L is the longitudinal extension of the beam, and H(·) denotes the Heaviside function. The cylindrical symmetry allows to determine via Gauss theorem the potential energy V sc (r) of a test particle of charge q 0 embedded in this charge distribution. It is a simple exercise to show that where C ≈ 0.577 is the Euler constant, and is the exponential-integral function. If we suppose that the test particle is acted upon by this space-charge potential, and by the external, cylindrical harmonic oscillator potential (B1) -which, without space charge, would keep the particle in the Gaussian state with variance σ 2 -the total potential to be inserted in the Schrödinger equation (16) reads In terms of the usual dimensionless quantities we can also write A plot of the dimensionless potential v(s) for γ = 3 is shown in Fig. 12.
For stationary solutions of the form Eq. (16) becomes and due to the cylindrical symmetry of V it reduces to Considering the dimensionless quantities we finally obtain the reduced equation Despite its simple graphical behavior displayed in Fig. 12 the potential v(s) is complicated enough to bar the hope of exactly solving Eq. (89) in a simple way. Reasonable approximate solutions can be obtained by means of the variational method of Rayleigh and Ritz. We thus look for the minimum of the average energy whose reduced and dimensionless form is The minimization can be carried out for different values of the coupling parameter γ, but the frozen core model is realistic only for small values of γ: a strong coupling with a frozen core would have the paradoxical effect of totally expelling the distribution function of the test particle from the center of the beam, while keeping the Gaussian core always concentrated in the center. To elaborate an example we have chosen γ = 3 and a test function of the form The values of the parameters that minimize the energy functional turn out to be The corresponding optimal radial density distribution s w 2 (s) is shown in Fig. 13 where it is also compared with the Gaussian solution (dashed line) corresponding to the case of a purely harmonic potential (γ = 0). Notice that this density distribution has a form similar to that of the self-consistent solutions elaborated in Section IV for values of γ of about 10 (see Fig. 2). It has a node which is not visible in Fig. 13. We zoom up its plot in Fig. 14. This behavior suggests that the halo distribution could be described, in first approximation, as a Gaussian core distribution plus a small ring of particles surrounding it and constituting the halo. Thus, in the following we will simply assume a specific form of a beam with a halo without deriving it as an effect of space charge interactions.
Since it is not clearly established that a halo can be due only to space-charge effects, starting with a realistic ring distribution and trying to understand the dynamics that can produce it could be very useful in this respect. The techniques introduced in Section III will be instrumental to this end. For a three-dimensional, cylindrically sym-metric beam we introduce the normalized radial density distribution which is composed of a Gaussian core with variance σ 2 (the simple harmonic oscillator ground state), plus a ringlike distribution whose size is fixed through the two parameters p > 0 and q ≥ 0. The parameter 0 ≤ A ≤ 1 is the relative weight of the two parts. This density distribution has the required form for suitable values of the parameters, but, at variance with the two previous examples of approximate distributions (see Figs. 11 and 14), it has no nodes. This is convenient for two main reasons: first because it is a rather general requirement for a ground state to have no nodes (this fact is a rigorous theorem for one-dimensional systems). Moreover, it has been shown in Refs. [30,32,33] that stationary distributions without nodes are also attractors for every other possible (non extremal) initial distribution: a property that will be useful in a future discussion of the possible relaxation of the system toward a stable beam halo. For this cylindrically symmetric case the expressions (30) and (37) for the velocity field and the potential become Moving to dimensionless quantities we have A three-dimensional plot of the dimensionless density distribution is given in Fig. 15 where we have chosen p = 1, q = 24 and A ≈ 0.49. In this example the importance of the halo ring has been exaggerated with respect to  the real case in order to make the effect clearly visible. The explicit analytic expression of b(s) and v(s) is rather lengthy and not particularly illuminating: we simply show the behavior of these functions, for the same values of the parameter for the previous example, respectively in Fig. 16 and Fig. 17, where they are also compared with the corresponding dimensionless quantities (−s and s 2 ) for the ground state of the harmonic oscillator. The function v(s) is the potential that in the stochastic model keeps the beam in the stationary state with a halo ring. In order to illustrate the main effects involved, let us give here the simpler formulae relative to the 1d case. From now on we will consider only 1d processes denoting by x one of the transverse space coordinates. We assume that the longitudinal and the transverse beam dynamics can be deemed independent, with the further simplification of considering decoupled evolutions along the transverse directions x and y. Under these conditions the density distribution with a halo ring now reads so that the corresponding velocities and potential are In terms of dimensionless quantities we then have for the distribution while the forward velocity and the potential are In Figrs. 18, 19, and 20 we respectively show the density distribution, the velocity field, and the potential for p = 1, q = 24 and A = 0.85. As in the three-dimensional case, the importance of the halo ring has been exaggerated with respect to the real case in order to make the effect clearly visible.

VI. CONCLUSIONS
In this paper we have presented a dynamical, stochastic approach to the description of the beam transverse distribution in the particle accelerators. In the first part we have described the collective beam dynamics in terms of time-reversal invariant diffusion processes (Nelson processes) which are obtained by a stochastic extension of the least action principle of classical mechanics. The diffusion coefficient of the process is identified with a lower bound for the emittance of the beam as discussed in the Section IV. The collective dynamics of beams is then described by two nonlinearly coupled hydrodynamic equations. This set of equations shows some formal analogies with the effective Ginzburg-Landau and Gross-Pitaevskii descriptions of the collective dynamics of self-interacting quantum many body systems in terms of classical, nonlinear Schrödinger equation. We have shown that, in the framework of SM, it is possible to have transverse distribution which show a broadening and an emittance growth, typical of the halo formation. The interest of this approach lies in the fact that its dynamical equations allow us, at least in principle, to determine the possible control potentials that are responsible for the stationarity of the halo. A few introductory examples of these potentials and distributions have been discussed in the section V. The controlling potentials can be engineered by suitable tuning of the external electromagnetic fields. In a previous paper [20] we have considered evolutions that drive the beam from a less collimated to a better collimated state, and we have furthermore shown that this goal can also be achieved without increasing the frequency of the betatron oscillations which can in fact be independently controlled during the evolution. In forthcoming papers we plan to elaborate on the extension of these techniques to the problem of engineering suitable time-dependent potentials for the control and the elimination of the beam halo. In particular, it will be shown that the transition functions of the Nelson processes can be exploited to control these evolutions. We also plan to extend the analysis of the present paper to different, non Gaussian beam transverse distributions which could better describe the existence of halo losses in terms of large ratios of the maximum displacement to the RMS size of the beam.

APPENDIX A: COUPLED FIELD EQUATIONS
Starting from the Maxwell equations for the space-charge field felt by our particle in the beam ∇ · E sc (r, t) = ρ sc (r, t) ǫ 0 ∇ × E sc (r, t) = − ∂ t B sc (r, t) ∇ · B sc (r, t) = 0 ∇ × B sc (r, t) = µ 0 j sc (r, t) + ∂ t E sc (r, t) c 2 and introducing the corresponding electro-magnetic potentials (A sc , Φ sc ) from B sc (r, t) = ∇ × A sc (r, t) E sc (r, t) = −∂ t A sc (r, t) − ∇Φ sc (r, t) with the gauge condition ∇ · A sc (r, t) + 1 c 2 ∂ t Φ sc (r, t) = 0 , through the usual procedure we get the wave equations It is also well known (see for example [16], chapter XV) that the Schrödinger equation for spinless, charged particles in an electro-magnetic field (A, Φ) is For a charged particle in the beam the electro-magnetic field is the superposition of the space-charge potential (A sc , Φ sc ) plus an external, control potential (A ext , Φ ext ), and hence we obtain the Schrödinger equation (43) that we rewrite here for convenience: It is apparent now that (A1), (A2), (A3) and (A4) constitute a self-consistent system of coupled, non linear differential equations for the fields ψ, A sc and Φ sc . Since we are in a reference frame comoving with the beam, in the following we will suppose to deal only with stationary wave functions of the form (42). In this case, as already remarked, we get j sc = 0, so that we can take the trivial solution A sc = 0 of (A2). The gauge condition (A1) then implies that ∂ t Φ sc = 0 and the wave equation (A3) reduces itself to the Poisson equation for the electrostatic potential. Finally, by supposing also A ext = 0, our system is reduced to only two coupled, non linear equations for the couple (u, Φ sc ), namely where we have defined |u(r, t)| 2 d 3 r .
Then if we pass to the potential energies by putting V ext (r, t) = q 0 Φ ext (r, t) , V sc (r, t) = q 0 Φ sc (r, t) our equations are reduced to Remark that here the wave functions ψ and u, while certainly normalizable, are not supposed in general to be already normalized. However the equations (A5) and (A6), albeit non linear, are apparently invariant for the multiplication of u by an arbitrary constant. The equations (A5) and (A6) will be used in the Section IV.

APPENDIX B: CYLINDRICAL HARMONIC OSCILLATOR
In the model discussed in this paper we suppose that the external potential V e is a cylindrically symmetric, harmonic potential with a proper frequency ω, so that the ground state without space-charge interaction will have a variance .
This means that, in the usual cylindrical coordinates {r, ϕ, z} with r = x 2 + y 2 , we have V e (r) = m 2 ω 2 r 2 = α 2 8mσ 4 r 2 (B1) and the corresponding Schrödinger equation with periodic conditions at z = ±L/2 (for a bunch of length L) has the eigenvalues E nk = (n + 1) αω + k 2 2πσ L 2 αω and the normalized eigenfunctions ψ nkl (r, ϕ, z, t) = R nl (r) Φ l (ϕ) Z k (z) e −iE nk t/α with R nl (r) = We suppose that, by neglecting the space-charge interaction and in a comoving frame of reference, the system will be correctly described by the ground state ψ 000 (r) = e −r 2 /4σ 2 σ √ 2πL (B2) associated to the eigenvalue E 00 = αω.