Stochastic effects in real and simulated charged particle beams

The Vlasov equation embodies the smooth field approximation of the self-consistent equation of motion for charged particle beams. This framework is fundamentally altered if we include the fluctuating forces that originate from the actual charge granularity. We thereby perform the transition from a reversible description to a statistical mechanics' description covering also the irreversible aspects of beam dynamics. Taking into account contributions from fluctuating forces is mandatory if we want to describe effects like intrabeam scattering or temperature balancing within beams. Furthermore, the appearance of ``discreteness errors'' in computer simulations of beams can be modeled as ``exact'' beam dynamics that is being modified by fluctuating ``error forces''. It will be shown that the related emittance increase depends on two distinct quantities: the magnitude of the fluctuating forces embodied in a friction coefficient $\gamma$, and the correlation time dependent average temperature anisotropy. These analytical results are verified by various computer simulations.

The Vlasov equation embodies the smooth field approximation of the self-consistent equation of motion for charged particle beams. This framework is fundamentally altered if we include the fluctuating forces that originate from the actual charge granularity. We thereby perform the transition from a reversible description to a statistical mechanics' description covering also the irreversible aspects of beam dynamics. Taking into account contributions from fluctuating forces is mandatory if we want to describe effects like intrabeam scattering or temperature balancing within beams. Furthermore, the appearance of "discreteness errors" in computer simulations of beams can be modeled as "exact" beam dynamics that is being modified by fluctuating "error forces". It will be shown that the related emittance increase depends on two distinct quantities: the magnitude of the fluctuating forces embodied in a friction coefficient γ, and the correlation time dependent average temperature anisotropy. These analytical results are verified by various computer simulations.

I. INTRODUCTION
Analytical approaches to the dynamics of charged particle beams that are based on the Liouville -or equivalently on the Vlasov -equation do not include effects due to the actual charge granularity. A variety of beam phenomena are adequately described by this continuous description. As the first example, we cite the pioneering work of I.M. Kapchinskij and V.V. Vladimirskij [1] covering the description of beam transport under space charge conditions. As a second example, we may quote the well-understood transient effects that occur if a beam is launched with a non-self-consistent phase space density profile [2][3][4]. Furthermore, the various kinds of parametric resonances and instabilities that may occur in the course of beam propagation through focusing lattices and storage rings have been successfully tackled on the basis of a perturbation analysis of the Vlasov equation [5][6][7][8][9].
Despite all these achievements, there is still an important class of beam phenomena the analysis of which leads beyond the scope of the Vlasov approach. Due to the invariance of Vlasov's equation with respect to time reversal [10], we must realize that it restricts the analysis to only reversible aspects of beam dynamics. However, a reversible, continuous description of beam dynamics no longer applies if the individual interactions of the point charges must be taken into account. Effects of elastic Coulomb scattering like the well-known phenomenon of intrabeam scattering [11] observed for intense beams that circulate in storage rings, or the process of temperature balancing within a charged particle beam -commonly referred to as beam equipartitioning -fall into this category. In order to include these irreversible effects into our analytical description of beams, the Vlasov approach must be generalized appropriately [12][13][14][15]. This will be achieved by switching from a deterministic to a statistical treatment of beam dynamics, namely by separating the actual forces that act on the beam particles into a smooth and a fluctuating component. We will review this transition in detail in Sec. II.
Owing to the fact that an analytical solution for the problem of particles interacting by Coulomb forces does not exists, computer simulations have become the tool of choice for the study of charged particle beams. In these studies, the actual beam is represented by an ensemble of simulation particles. A simulation thus means to numerically integrate the coupled set of equations of motion constituted by the particle ensemble and the beam optical lattice. Although the equations of motion of individual particles are invariant with respect to time reversal, the evolution of the particle ensemble is inevitably rendered irreversible because of the limited accuracy of numerical methods. Therefore, a simulation based on individual particles can never be a strict realization of a solution of the associated Vlasov equation.
The idea pursued in this article is to describe these "computer noise" effects analogous to random force effects emerging within the granular charge distribution of a "real" beam. We can then interpret the beam simulation results within the framework of the generalized Vlasov approach, which will be reviewed in Secs. III to VIII. This allows us to separate effects caused by the specific realization of the computer simulation from the "real beam" physics. The onset of irreversibility in a computer simulation of a charged particle ensemble will be demonstrated in Sec. IX A. We numerically transform a beam forward a specific amount of time, followed by the backward transformation to its starting point. The reversible transient effect of "initial emittance growth" -occurring for beams with non-self-consistent phase space densities, as described by the Vlasov equation -is rendered irreversible because of the accumulated action of "error" forces that inevitably accompany our numerical calculations.
In a second simulation example presented in Sec. IX B, the joint occurrence of reversible and irreversible effects within simulated charged particle beams is visualized. For a specific time span after a numerical time reversal, the beam evolution behaves reversible. After this, the irreversible "computer noise" effects prevail, indicated by a sharp change of the sign of the emittance growth rate.
In Sec. IX C, we will analyze the numerical emittance growth factors obtained for different focusing lattices, matching conditions, and number of simulation particles. It will be shown that the specific emittance growth rates emerging in these simulations can indeed be explained within the frame-work of the generalized Vlasov approach. We also investigate the scaling of these emittance growth factors with the number of particles used in the simulation in order to distinguish "computer noise" related effects from those occurring within a real beam.

II. LANGEVIN EQUATION
We start our analysis reviewing the single particle equation of motion for a set of charged particles interacting through Coulomb forces within the co-moving beam frame with m the particle mass and q its charge, F ext denoting the external force field and E sc the total electric self-field generated by all other particles. If a total of N particles of the same species is given, the N-body distribution function contains the complete information on the state of the system. Eq. (1) together with the knowledge of (2) defines a "reversible" system, hence a system that is completely determined and does not contain any sources for loss of information. Without such losses, the system can be transformed to any instant of time back and forth. In principle, all effects occurring in charged particle beams can be extracted from the time integration of Eq. (1). Nevertheless, this picture is not adequate for the description of real N-body systems if N is very large, since the condition that the initial ρ is precisely known can never be fulfilled. In addition, the detailed knowledge of ρ(t) is usually not useful. Therefore, a statistical description of the time evolution of the particle ensemble is appropriate. This description must be consistent with exact solutions of Eq. (1) for a large number of particles N.
On the single particle level, a statistical description means to replace the exact, fine-grained force contained in Eq. (1) by its smoothed coarse-grained average force. The fine-grained aspect of the particle motion is then modeled by an additional fluctuating force F L . As pointed out by Jowett [16], this concept constitutes "an attempt to describe the effects of the neglected microscopic degrees of freedom". In order not to introduce a systematic error into the statistical description of the N-particle ensemble, this force must vanish on the ensemble average In a statistical description, we must not conceive F L (v, t) to be an ordinary vector function but a quantity that has only statistically defined properties. Fluctuating forces of this nature are usually referred to as "Langevin forces" [17]. In performing the transition from an "exact" fine-grained description of the evolution of ρ(t) according to Eq. (1) to a statistical description of this evolution, not only the fluctuating Langevin force F L (v, t) but also a force that is referred before after closest encounter and to as "dynamical friction" force F fr (v, t) must be introduced. For repelling forces, the mechanism of dynamical friction is sketched in Fig. 1. We observe that the deceleration of the leftmost particle in horizontal direction before its closest encounter with the other particles is greater than its acceleration afterwards. This means that a net deceleration, hence a friction occurs. As is easily verified, the same is true for attracting forces.
Owing to the statistical description of the N-particle ensemble, the self-field appears now as a smooth function of x and t that is equivalent to an external force field. The stochastic counterpart of the deterministic single particle equation of motion (1) can now be written as containing the smooth part of the self-force E sm sc (x, t), the dynamical friction force F fr (v, t), and the fluctuating Langevin force F L (v, t). As usual, we made the reasonable assumption that the stochastic effects in our statistical description are independent of the "external" force functions F ext (x, t) and qE sm sc (x, t). This means that the Langevin force F L as well as the friction force F fr do not depend on the position x in real space.
Each particle encounters a specific realization of the Langevin force F L (v, t). These forces are defined by their statistical properties only, a direct integration of Eq. (4) is thus not possible. On the other hand, a deterministic equation of motion for the phase space probability density f (x, v, t) can be derived on the basis of Eq. (4). This topic will be the subject of the next section.
The statistical mechanics equation (4) is supposed to provide an equivalent description of the dynamics of an N particle ensemble as the "exact mechanics" equation (1). Therefore, the magnitudes of Langevin and friction forces contained in Eq. (4) are completely determined by the force fluctuations that follow from the actual charge granularity, as described by Eqs. (1) and (2). This is no longer true if Eq. (4) is used as the analytical basis to interpret results of computer simula-tions of N particle ensembles. If we model the impact of the generally limited accuracy of numerical methods by Langevin and friction force terms that act on the simulation particles in addition to "true" forces experienced by the "real" beam, the magnitude of the stochastic forces depends on the specific nature of the simulation code. In other words, the results of beam simulations must be regarded as solutions of Eq. (4) with the magnitudes of F fr and F L being determined by the particular scheme of simulation, defined by the time step width of the numerical integration, the number of simulation particles, the computer's word length, and others [18].
In recent studies [19], the close relation between results of molecular dynamics (MD) simulations and coarse-grained dissipative particle dynamics has been worked out. Similarly, an approach based on the Fokker-Planck equation in order to explain numerical emittance growth effects observed in particle-in-cell (PIC) simulations has been presented earlier [14,15]. In the following sections, we will review and further extend this analysis.

III. FOKKER-PLANCK EQUATION
We define q ≡ (x, v) as the position vector in the 6dimensional µ-phase space. If the function f (q, t) represents a normalized phase space probability density, f dq provides us with the probability of finding a particle inside a volume dq around the phase space point q at time t. In these terms, the generalization of Eq. (4) can be written aṡ with smooth functions K i (q, t) and the random variables Γ i (q, t) vanishing on the ensemble average. We now assume the random variables Γ i (q, t) to be Gaussian-distributed and their time correlation proportional to the δ-function Under these conditions, the Kramers-Moyal expansion for ∂ f (q, t)/∂t terminates after the second term [20][21][22]. The expansion with only the first and second term is usually called Fokker-Planck equation with the Fokker-Planck operator L FP given by We observe that the coefficients Q i j are determined by the amplitude of the δ-correlated noise functions Γ i according to (6), whereas the K i are defined by Eq. (5). Consequently, Eq. (7) represents the deterministic equation of motion for the probability density f (q, t). It is uniquely determined by the coupled set of Langevin equations (5) provided that (6) holds.
In terms of the special Langevin equation (4), the Fokker-Planck operator L FP reduces to with F tot,i defined as the sum of all non-Langevin forces and the diffusion coefficients D ii defined by The off-diagonal terms of the diffusion matrix D i j vanish since the Langevin forces in Eq. (4) are not correlated for different degrees of freedom. We further note that the friction forces F fr,i must always be decelerating. This means that F fr,i must change sign if v i does, hence must be an odd function of v i . With regard to Eq. (9), it follows that the diffusion coefficients of Eq. (8) must be even functions of the v i A Fokker-Planck equation that describes the evolution of the probability density f appertaining to the stochastic motion of particles in external force fields if often referred to as Kramers equation. As will be shown in the next section investigating equilibrium solutions of Eq. (7) with the Fokker-Planck operator (8), the diffusion coefficients D ii (v i , t) are uniquely determined by the friction forces F fr,i (v i , t).

IV. FOKKER-PLANCK COEFFICIENTS UNDER TIME REVERSAL
If we perform a transformation that reverses the direction of time flow, the positions x i and hence all quantities that only depend on the positions do not change sign. In contrast, the velocities v i do change sign, which means that quantities depending on the v i may change sign under time reversal. We may thus separate the components of the Fokker-Planck operator (8) with respect to their behavior under time reversal The "reversible" operator L rev is defined to consist of those components of (8) that change sign under time reversal The smooth self-field E sm sc is obtained from the real space projection of the probability density f (q, t) via Poisson's equation. The components that do not change sign constitute L ir Here we made use of (10), which states that under time reversal F fr,i changes sign, whereas D ii does not change sign. The external forces F ext,i have been assumed to be not velocity dependent.
Since ∂ f /∂t changes sign on time reversal, a Fokker-Planck equation with only L rev remains unchanged if the direction of time flow is reversed. It therefore describes the reversible transformation of the probability density function f (x, v, t). This means that earlier states are fully restored if a reversed time integration of Eq. (7) with L FP ≡ L rev is carried outjust like a movie that is reversed at some instant of time t 0 . Correspondingly, L ir describes exactly those effects that do not depend on the direction of the time flow. In other words, it describes the irreversible aspects of the particle motion. With L ir = 0, Eq. (7) is commonly referred to as Vlasov equation.

V. EQUILIBRIUM DISTRIBUTIONS IN AUTONOMOUS SYSTEMS
If the external force F ext (x) contained in Eq. (8) is not explicitly time dependent, a stationary solution L FP f st = 0 may exist. If it exists, it can always be written in the form Obviously, all irreversible currents must vanish for f = f st to be stationary. With L ir,i given by (12), this means, explicitly, Eq. (15) states that for given φ st , the diffusion function D ii (v i ) is uniquely determined by the friction force function F fr,iand vice versa. This mutual dependency of the diffusion effects -driving a system away from its steady stateand damping effects that cause the decay of these deviations makes up the physical substance of "fluctuation-dissipation theorems". In agreement with Eq. (10), we express the friction force function F fr,i , and the diffusion function D ii as odd and even power series in v i , respectively Here we assumed the coefficients a k , b k not to depend on x and the degree of freedom i -in agreement with the precondition that the stochastic effects are not influenced by the external forces.
With (16), we find that Eq. (15) can only be fulfilled if φ st is a quadratic function of the v i . Therefore, φ st may always be separated as the angle brackets denoting the respective averages over the phase space density function: a = a f dτ. The quantity v 2 i thus embodies the ensemble average of the squares of all particle velocities, also referred to as the second moment of the velocity v i of the equilibrium distribution f st . In a state of equilibrium these moments must agree for each degree of freedom, hence can be identified with the equilibrium temperature T eq according to with k denoting Boltzmann's constant.
Inserting (17) into the Fokker-Planck equation (7,8), the generalized potential ψ st (x) follows from In final form, the equilibrium probability density of the Fokker-Planck equation (7) reads We summarize that the equilibrium distribution (20) follows directly from the assumption that the stochastic component of the particle motion is caused by Gaussian-distributed Langevin forces with a time correlation function proportional to the δ-function, regardless of the dependency of the friction and the diffusion coefficients on the v i . For a given temperature T eq , the spatial probability function following from ψ st is uniquely determined by the external force F ext , and the stationary self-field E sm sc . Together with the unique velocity distribution, the entire phase space probability density function is uniquely determined, which means that no other equilibrium distribution of (7) exists -in contrast to Vlasov systems where friction as well as diffusion effects vanish. If the external force function F ext (x) does allows for an equilibrium, and if the friction is not negligible, arbitrary non-equilibrium density functions always settle down to a unique equilibrium. This is what we observe in long-term simulations of charged particle beams [4,23]. Regardless of our initial phase space filling, we always end up with a Gaussian velocity distribution if no resonance effects are involved.

VI. BOLTZMANN ENTROPY GROWTH
As a simple ansatz, we truncate the power series (16) for the friction force F fr,i (v i ) after the cubic term in v i With (15) and (17) we then immediately obtain the diffusion coefficient In order to quantify the impact of friction and diffusion on the phase space probability density function f (x, v, t) we now define the negative Boltzmann entropy S (t) as [15,24] We easily convince ourselves that S (t) remains unchanged if only the reversible part (11) of the Fokker-Planck operator L FP drives the time evolution of f If a phase space density f does not represent an equilibrium state, the irreversible part (12) of L FP must therefore account for a non-constant S (t). We approximate an arbitrary non-equilibrium phase space probability density distribution with g(x, t) the non-equilibrium real space probability density, T i (t) the non-equilibrium "temperature" appertaining to the ith degree of freedom, and v inc i the incoherent contribution to the total particle velocity v i In energy units, the temperature kT i (t) = m v 2 i inc is defined as the incoherent kinetic beam energy. If a phase space distribution f is in non-equilibrium state, the total kinetic energy of a beam particle consists of both, a coherent as well as an incoherent part. The temperature kT i (t) is thus obtained by subtracting the coherent kinetic energy from the total kinetic beam energy. kT i (t) then evaluates to with ε i (t) denoting the root mean square (RMS) emittance of the particle ensemble in the i-th degree of freedom. With the non-equilibrium density function (24), and L ir according to (12), we obtain for the time derivative of the entropy (23) [15] Inserting our approximations for the friction force (21) and the diffusion coefficient (22), this expression simplifies to We note that the cubic term of the friction force ansatz (21) does not modify the form of Eq. (27). The equilibrium temperature T eq has been defined in Eq. (18) for autonomous systems. Before applying Eq. (27) to non-autonomous systems, we must discuss how to define appropriately the "equilibrium temperature" T eq in these systems. This will be the subject of the following section.

VII. EQUILIBRIUM TEMPERATURE IN NON-AUTONOMOUS SYSTEMS
For equilibrium distributions in autonomous systems as discussed in Sec. V, the equilibrium temperature T eq contained in (27) is a constant of motion. For real focusing systems, i.e. non-autonomous systems with the external force F ext (x, t) being explicitly time dependent in the beam system, such a constant does not exist. Under these circumstances, an equilibrium temperature must be defined analogically as the mean temperature averaged over a correlation time interval δτ. For t ≥ δτ this means The length of the time interval δτ must depend on the amplitude of the Langevin forces acting within our system, i.e. on the effective friction coefficient γ. For γ → 0, hence for reversible systems, the memory on earlier states is never lost. Therefore the average temperature depends on the whole time interval elapsed since launching of the beam On the other hand, if γ is very large, the memory on earlier states is rapidly lost, which means that δτ → 0. The equilibrium temperature is then given by the instantaneous temperature T (t), namely the arithmetic mean of the temperatures In real systems -as well as in numerical simulations that take into account the space charge forces -γ is usually small but finite. As a consequence, the equilibrium temperature definition (28) must be used in our analytical description. Inserting (28) into (27), the entropy change over the time span δτ follows by integration with the dimensionless temperature anisotropy coefficient A δτ (t) defined as The total entropy change at multiples of the correlation time interval t = Mδτ is then given by the sum over all elementary intervals Eq. (31) states that the entropy growth a particle ensemble experiences is determined by both, the friction coefficient γ, and the average temperature imbalance according to Eq. (30). It represents what we expect for a Markov process: the total entropy increase is given by a sum of elementary increases that do not depend on each other. The temperature imbalance is determined by the optical properties of the focusing lattice, the matching conditions of the beam, and the finite correlation time δτ > 0 that measures the typical duration of an elementary scattering event. δτ thus defines the time span over which the instantaneous temperatures T i (t) must be averaged in order to obtain the reference temperature (28). The actual value of γ can be estimated for charged particle beams by averaging the elementary process of a binary Coulomb scattering event over the impact parameters, and subsequently over the beam's velocity distribution [25].
For a beam propagation that is being influenced by a large number of internal scattering events, irreversibility just means that the time-reversed evolution is highly improbable. This perception of irreversible behavior of a dynamical system cannot be applied directly to numerical noise effects in computer simulations. For this case, a more appropriate interpretation of the entropy has been given by Shannon [26], in the context of his founding of "information theory". Within this framework, the entropy measures the amount of missing information about the state of the system in question. The entropy growth thus quantifies the amount of information an observer looses during the system's time evolution. In case of our simulations, the necessarily limited accuracy of numerical methods accounts for the mechanism that causes a loss of information. We can therefore no longer regard γ to represent a "real" physical process if we make use of Eq. (31) in order to interpret beam simulation results. Unfortunately, it appears to be impossible to analytically estimate the "computer noise" induced γ in so far as it depends on the particular realization of the simulation. The actual value of γ must therefore be extracted from the simulation data itself by comparing the obtained emittance growth factors for different average temperature imbalances. To this end, the relation between entropy change (31) and irreversible emittance growth must be established. This will be worked out in the following section. Afterwards, the discussion to what extend Eq. (31) can also be used to interpret computer simulation results will be presented in Sec. IX.

VIII. EMITTANCE GROWTH ASSOCIATED WITH L ir
A second order moment analysis of Eq. (7) yields the following set of coupled moment equations [14] for each phase Calculating the time derivative of the RMS emittance (25), and inserting the above derivatives of the second moments, we find that three distinct sources for the RMS emittance change can be distinguished namely the external field contribution, the contribution related to the smooth space charge fields and the contribution due to the Langevin forces described by the irreversible part (12) of the Fokker-Planck operator.
Collecting together the terms containing the external force F ext , we find that its contribution to the change of the RMS emittance vanishes if it is a linear function of the spatial coordinates This applies to all our computer simulations addressed in Sec. IX. The terms containing the smooth space charge field E sm sc sum up to Writing this equation for all three spatial degrees of freedom, the electric field terms together form the physical quantity of "free field energy", i.e. the difference of the actual field energy W and the field energy W u of the equivalent uniform charge density [3,27,28] In Sec. IX, we will show by numerical simulation that the exchange of RMS emittance and "free field energy" is indeed a reversible process. The third contribution to the change of the RMS emittance emerges from the irreversible Fokker-Planck operator (12) We restrict ourselves to a friction force (21) linear in the v i . The related diffusion coefficient D ii then follows as Since γ 1 is always positive, the equation states that the RMS emittance ε i (t) increases, as long as the beam temperature T i lies below the equilibrium temperature T eq -and vice versa. This is what we expect for a temperature balancing process. We observe that the right hand sides of Eqs. (36) and (27) agree for γ = γ 1 . This means that the growth of the Boltzmann entropy (23) is related to the irreversible emittance growth according to With Eq. (31) and ε = 3 √ ε x ε y ε z , the related irreversible growth of the total RMS emittance at t = Mδτ is given by The validity of Eq. (37) will be verified by numerical simulations in Sec. IX C. Beforehand, we demonstrate the emerging of irreversibility due to "computer noise" effects that necessarily accompany numerical simulations of charged particle beams.

A. Initial Emittance Change
In this section, we present results of a simulation code that numerically integrates the reversible equations of motion (1), starting from a sharply known initial phase space filling (2).
As the first example, we simulate the transient effect of "initial" emittance change -a well understood phenomenon in the realm of charged particle beams [2][3][4]. It occurs if a beam is injected into an ion optical system in a non-equilibrium state. As a consequence, the phase space density f changes rapidly until an average equilibrium state is reached. As part of this process, the charge density profile adjusts itself to the external forces, which means that the electrostatic field energy constituted by the initial charge density is modified. This process is accompanied by a change of the RMS emittances. With the emittance definition of Eq. (25), the "exchange" of emittance and field energy is described quantitatively by Eq. (33). We notice that this equation is invariant with respect to the reversal of time -in agreement with the fact that it follows directly from the Vlasov equation ∂ f /∂t = L rev f . The conclusion that the initial charge density adjustment is indeed a reversible process is verified by the simulation results displayed in Fig. 2. Because of a peaked initial charge density defined by the initial phase space filling, field energy is released immediately after launching the beam, followed by a damped oscillation around the field energy of the equilibrium density profile. This process is accompanied by a variation of the RMS emittance according to Eq. (33). Nonetheless, if the direction of the beam transformation is reversed after 5 periods, the initial non-equilibrium state is recovered. This is no longer true if the forward transformation exceeds a certain amount of periods. Fig. 3 shows the emittance variations obtained from the similar simulation as in the previous case, but with the forward and the subsequent backward transformation now extending over 20 focusing periods. Obviously, the RMS emittance does no longer return to its initial value, but keeps on oscillating around the level associated with the self-consistent state. After having been transformed over a certain time span, the simulated beam has evolved in a way that cannot be reversed anymore. This behavior can be explained if we interpret the numerical inaccuracies that inevitably accompany all our simulations to arise from forces F fr and F L of Eq. (4) that additionally act on the simulation particles. As outlined in Sec. III, these forces induce a nonvanishing irreversible component of the Fokker-Planck operator L ir 0 in the equation of motion for the phase space density function f . Accordingly, we obtain a gradual loss of "memory" of previous states of the phase space density f (x, v, t). This loss of information during the beam transformation gradually renders all emittance growth effects irreversible, regardless of their specific nature being reversible or irreversible. In other words, after having passed 20 cells, in our particular simulation example the beam does not "remember" anymore that a reversible emittance growth has taken place right after the start of the simulation. As the consequence, the specific initial beam state that is assigned to a lower emittance cannot be recovered. Instead, the more probable self-consistent state associated with the increased emittance is kept. B. Emittance Growth due to Anisotropic Focusing Fig. 4 displays the emittance growth curve obtained for a simulation of beam transport through a fictitious lattice that focuses the beam anisotropically in all three spatial directions. This anisotropic focusing enforces a non-vanishing beam temperature anisotropy coefficient (30) throughout the lattice. As stated in the previous subsection, numerical simulations are always accompanied by "discreteness errors", which may be described by additional Langevin and friction force terms in the single particle equation of motion (4). With regard to our stochastic description of beams, this means that both quantities, the finite temperature anisotropy A δτ as well as a nonvanishing effective friction coefficient γ 1 induce a specific irreversible growth rate of the beam emittance according to Eq. (37). The amplitudes of the Langevin forces emerging in beam simulations are larger than the corresponding charge granularity forces occurring within a real beam. We thus encounter a larger friction coefficient γ 1 within our simulation procedure than we would expect from an analytical estimation of γ 1 for a real beam [25]. Consequently, the numerically obtained emittance growth rate must be larger than the intrabeam scattering related growth rate for a real beam. In order to verify that the emittance growth rate obtained in our simulation of Fig. 4 is indeed caused by the action of Langevin forces, we again numerically reverse the direction of the time integration of the single particle equations of motion. For a small number of focusing periods -in this particular case for about 6 periods -the beam evolution "behaves" reversible, visualized by the "roll back" of the emittance curve that covers exactly the forward transformation graph. Having exceeded this time span, the beam's evolution in rendered irreversible, indicated by the sharp change of sign of the slope of the emittance curve. This behavior is directly related to the presence of Langevin forces which, as noted earlier, induce a non-vanishing irreversible part of the Fokker-Planck operator (12). As shown in Sec. IV, this operator describes exactly those aspects of the evolution of a beam dynamical system that do not change if the direction of time flow is reversed. This is what we observe in Fig. 4 after the reversible phase of the back transformation. The emittance growth rate persisting during the irreversible phase of the back transformation agrees exactly with the growth rate obtained along the preceding forward transformation.

C. Scaling Law for the Friction Coefficient
We conclude this article by presenting simulation results aimed at investigating the range of validity of Eq. (37). As shown in Sec. VIII, this equation relates the logarithm of the irreversible emittance growth to the product of the effective friction coefficient γ 1 and the average temperature anisotropy the beam experiences between 0 and t. Accordingly, the simulations comprise examples of beam tracking through various focusing systems and beam matching conditions, each of them inducing a specific temperature anisotropy A δτ along the beam line. The correlation time parameter δτ contained herein has been introduced in Eq. (28) as constituent part of the definition of the "equilibrium temperature" for non-autonomous systems. For a given friction parameter γ 1 , this correlation time δτ must be adjusted appropriately. Furthermore, similar simulations are performed with only the numbers of simulation particles being varied. We hereby modify the friction coefficient γ 1 in order to verify Eq. (37) independently for different noise levels.
Unlike the beam entropy, the beam's RMS emittance can be calculated directly from positions and velocities of all simulation particles. Unfortunately, as stated by Eq. (32), different mechanisms that all modify the RMS emittance in the course of the beam propagation must be distinguished in order to isolate the noise-related emittance growth effects. Because of the linear external focusing forces applied throughout the simulations presented here, growth effects due to a non-linear focusing force may not appear. The actually obtained emittance changes must therefore be attributed to either a variation of the beam's "free field energy" as described by Eq. (33), or the action of Langevin forces according to Eq. (37).
In our simulations, equivalent ensembles of macro-particles representing equivalent beams are tracked through three distinct fictitious external focusing geometries. At a time, the external focusing forces define (i) a system that isotropically focuses the beam in all three spatial directions, with the focusing forces acting continuously along the beam line; (ii) a system that again isotropically focuses the beam in all three spatial directions, but with focusing forces now acting periodically along the beam line; and (iii) a system that focuses the beam anisotropically in the three spatial directions, with focusing forces acting also periodically along the beam line.
For the isotropic focusing systems, the simulations are launched with both isotropic as well as anisotropic mismatch conditions. The strength of a particular mismatch will be quantified by the dimensionless "mismatch factor" ∆ [29], defined as with α cs i , β cs i , and γ cs i denoting the Courant-Snyder [30] functions for the actual and the matched beam, respectively.
In Fig. 5, we plot the RMS emittance growth factors obtained from the simultaneous numerical integration of N coupled single particle equations of motion (1) with N = 5000 macro-particles representing the beam. We observe that the emittance growth rates even for strong mismatch are much smaller if the beam stays isotropic -as given for isotropic external focusing in case of isotropic mismatch. On the other hand, if the beam's phase space symmetry is rendered anisotropic because of anisotropic focusing or mismatch, the emittance growth factors come out considerably larger. This behavior can be explained considering Fig. 6. It shows the evolution of the respective temperature anisotropy coefficients A δτ (t), defined by Eq. (30). As is easily understood by their definition, the anisotropy coefficients are much smaller for isotropic beam mismatch oscillations, compared to anisotropic oscillations of comparable strength ∆. According to Eq. (37), the anisotropy coefficients are directly related to the irreversible part of the RMS emittance growth. We therefore expect the actual emittance growth rates to follow the magnitude of the temperature anisotropy if emittance changes due to variations of the "free field energy" can be neglected. This is indeed what we observe in our simulations. As stated above, the friction factor γ 1 contained in Eq. (37) provides us with a global measure for the magnitude of the Langevin forces that act within the ensemble of beam particles. We conclude that γ 1 should be similar for equivalent ensembles of macro-particles, as defined in our simulations.
We now estimate the normalized friction coefficient γ 1 τ that is consistent with the emittance growth effects experienced in the simulations for the defined variety of temperature anisotropies. According to Eq. (37), γ 1 determines the amount of irreversible emittance growth that is obtained for a given temperature anisotropy. Since the actual emittance growth factors resulting from the evolution of the simulated particle ensemble reflect a mixture of reversible as well as irreversible effects, we cannot extract γ 1 directly from the simulation data. Instead, we plot in Fig. 7 the related factorsγ, defined byγ Under the condition that emittance changes due to variations of the "free field energy" can be neglected, hence that only the Langevin forces account for emittance changes,γ is identical with the global friction factor γ 1 . Otherwise, only the time average ofγ provides us with an approximation of γ 1 . With regard to Fig. 7 observed for the cases of isotropic beam symmetry indicate that "free field energy" contributions to the emittance change according to Eq. (33) take place. Following from a linear perturbation analysis of the envelope equations [31], the phase advance σ env,H of an isotropic beam "breathing" mode can be expressed in terms of the zero current single particle phase advance σ 0 and the related phase advance σ occurring in presence of space charge forces as As displayed in Fig. 5 and Fig. 7, we observe 5 beam core oscillations in about 17 focusing periods, which corresponds to a phase advance of about σ env,H = 106 • per cell for this mode. This number is in excellent agreement with Eq. (40), predicting a value of σ env,H = 105 • for our simulation parameters. The statement that the friction coefficient γ 1 is related to the magnitude of stochastic forces is confirmed by simulations performed with different number of particles while leaving all other simulation parameters invariant. Fig. 8 shows the normalized friction coefficientsγτ obtained from tracking of 10000 simulation particles. A comparison with the corresponding Fig. 7 -displaying the results for 5000 particlesshows that theγτ values are reduced. Finally, Fig. 9 shows the normalized friction coefficientsγτ obtained for 2500 simulation particles. Under these circumstances theγτ values come out larger compared to the previously cited cases, as expected. As a rough estimate, we find that γ 1 scales with the inverse square root of the number N of simulation particles γ 1 ∝ N −1/2 .

X. CONCLUSIONS
In this paper, we have reviewed the analytical description of emittance growth effects that are caused by the actual granularity of the charge distribution of particle beams -commonly referred to as intrabeam scattering. Within the same framework, the description of noise phenomena that necessarily accompany computer simulations of beams has been outlined. A formula has been derived that relates the expected emittance growth rate to both the magnitude of the noise force, and the average temperature anisotropy the beam experiences within a correlation time span δτ. We have presented computer simulations of beam transport through various focusing lattices and matching conditions, each of them enforcing a specific beam temperature anisotropy. It has been shown that the numerically obtained emittance growth rates can indeed be explained on the basis of this formula. The magnitude of the noise force has been shown to depend significantly on the number of simulation particles. As a consequence, we may identify the related emittance growth rates as "computer noise" artifacts.
Discreteness errors inevitably emerge in computer simulations of dynamical systems. Therefore, the actual time evolution of the simulated system always encloses irreversible aspects -even if the actually coded equations of motion are strictly reversible. In that sense, the simulation results can be regarded as exact solutions of a modified dynamical system that always comprises Langevin force terms. The magnitude of these "computer noise" related forces strongly depends on the specific realization of the simulation. The simulation may thus pretend effects that either do not occur at all within the "real" system in question, or that occur at different time scales.
The crucial point for the correct interpretation of computer simulations of beam dynamical systems is to keep in mind that the appearance of noise-related emittance growth depends on both the magnitude of the noise forces as well as the time averaged temperature anisotropy. Therefore, macroscopic emittance growth effects do not appear in cases where the temperature anisotropy within the system is negligible -even if the "computer noise" related forces are large. On the other hand, even if strictly periodic solutions of non-autonomous Vlasov-Poisson systems exist, the results of computer simulations of such systems will never be strictly periodic. Only if we take into account these subtle differences, we avoid misinterpretations of our computer simulation results.