Chaos and the continuum limit in nonneutral plasmas and charged particle beams

This paper examines discreteness effects in nearly collisionless N-body systems of charged particles interacting via an unscreened r^-2 force, allowing for bulk potentials admitting both regular and chaotic orbits. Both for ensembles and individual orbits, as N increases there is a smooth convergence towards a continuum limit. Discreteness effects are well modeled by Gaussian white noise with relaxation time t_R = const * (N/log L)t_D, with L the Coulomb logarithm and t_D the dynamical time scale. Discreteness effects accelerate emittance growth for initially localised clumps. However, even allowing for discreteness effects one can distinguish between orbits which, in the continuum limit, feel a regular potential, so that emittance grows as a power law in time, and chaotic orbits, where emittance grows exponentially. For sufficiently large N, one can distinguish two different `kinds' of chaos. Short range microchaos, associated with close encounters between charges, is a generic feature, yielding large positive Lyapunov exponents X_N which do not decrease with increasing N even if the bulk potential is integrable. Alternatively, there is the possibility of larger scale macrochaos, characterised by smaller Lyapunov exponents X_S, which is present only if the bulk potential is chaotic. Conventional computations of Lyapunov exponents probe X_N, leading to the oxymoronic conclusion that N-body orbits which look nearly regular and have sharply peaked Fourier spectra are `very chaotic.' However, the `range' of the microchaos, set by the typical interparticle spacing, decreases as N increases, so that, for large N, this microchaos, albeit very strong, is largely irrelevant macroscopically. A more careful numerical analysis allows one to estimate both X_N and X_S.


I. INTRODUCTION AND MOTIVATION
To what extent can a 'nearly collisionless' N -body system such as a nonneutral plasma or a charged particle beam be modeled by a smooth density distribution and a smooth bulk potential, either statistically or at the level of individual orbits? In particular, is there a well defined N → ∞ continuum limit? The idealisation of a smooth potential is extremely convenient, both conceptually and computationally. However, as noted, e.g., by beam dynamicists [1], it is not completely clear when -or ifsuch an idealisation is justified.
Even assuming that there is a well-defined continuum limit, how large must N be before the continuum approximation is justified? And to what extent can residual discreteness effects be modeled in the context of a Fokker-Planck description? The conventional Fokker-Planck description [2] [3] was formulated to extract statistical properties of orbit ensembles and distribution functions over long time scales, assuming implicitly that the bulk po-tential is regular. To what extent, then, can Langevin realisations [4] of a Fokker-Planck equation yield reliable information about individual orbits over comparatively short time scales, particularly if the orbits are chaotic?
This is an issue of practical importance for systems like high intensity charged-particle beams. In older, low intensity beams, the contribution of the space charge to the total potential is comparatively minor compared to the confining magnetic field, but in high intensity beams the space charge can become extremely important. To what extent, then, is one justified in idealising the space charge by a smooth density distribution that generates a smooth macroscopic potential? Is it, e.g., really legitimate to ignore discreteness effects entirely in a bunch comprised of ∼ 10 9 −10 10 protons with a transverse emittance of a few microns?
The situation is especially suspect in that experience with 'nearly collisionless' self-gravitating systems [5] indicates that the continuum limit must be subtle. One believes that, for N → ∞, orbits in an N -body system will in fact converge towards characteristics in the corresponding smooth potential. However, N -body orbits evolved in a density distribution corresponding to an integrable potential are typically 'very chaotic' with large positive Lyapunov exponents even though the integrable characteristics have vanishing Lyapunov exponent! For flows in smooth potentials, it is straightforward to distinguish macroscopically between phase mixing associated with regular versus chaotic orbits [6]. In particular, chaotic flows tend to mix exponentially fast, whereas regular flows only mix as a more modest power law in time. Will this distinction persist in N -body systems? One might, e.g., worry that if N -body orbits corresponding to integrable characteristics have large positive Lyapunov exponents, 'regular' flows comprised of such orbits would exhibit behaviour resembling chaotic phase mixing in smooth potentials! Understanding the bulk properties of phase mixing in N -body systems is important, e.g., in light of recent grid code simulations [7] indicating that chaotic phase mixing may be responsible for 'anomalous relaxation' observed in charged-particle beams, including the University of Maryland 'five beamlet' experiment [8].
A complete resolution of these issues will require an analysis of 'honest' (i.e., direct summation [9]) numerical integrations for systems comprised of very large particle number N , which is not yet practical. However, considerable insight may be derived by considering the simpler case of orbits and orbit ensembles evolved in 'frozen-N ' systems, i.e., time-independent N -body systems generated by randomly sampling a specified smooth density distribution. In particular, by comparing orbits in such frozen-N systems with characteristics with the same initial condition evolved in the corresponding smooth potential, one can quantify the extent to which, as a function of N , N -body orbits and smooth potential characteristics actually coincide. Such is the objective of this paper.
Section II focuses on the behaviour of orbits and orbit ensembles in frozen-N realisations of two simple potentials, one integrable and the other almost completely chaotic. The principal conclusions here are (i) that there is a well-defined macroscopic convergence towards the continuum limit, both for individual orbits and for orbit ensembles, and (ii) that discreteness effects can be extremely well-mimicked by Gaussian white noise in the context of a Fokker-Planck or Langevin description.
Section III considers the more realistic 'thermal equilibrium model' [10], well known from beam dynamics. For appropriate choices of parameter values, this model admits [11] large measures of both regular and chaotic orbits, so that one encounters a new feature [12], namely transitions between regular and chaotic behaviour triggered by discreteness effects. As for the simpler models considered in Section II, one observes clear distinctions between regular and chaotic phase mixing, although discreteness effects, again well modeled by a Fokker-Planck description, can be surprisingly important. Even when N is large, individual orbits can exhibit frequent changes in behaviour corresponding macroscopically to transitions between regularity and chaos; and the scaling implict in a Fokker-Planck description suggests that, for chaotic orbits, even a total particle number as large as N ∼ 10 9 may not be large enough to justify a continuum limit.
Section IV focuses on Lyapunov exponents and the meaning of chaos in N -body systems. The principal conclusion here is that two distinct 'types' of chaos can be present in the N -body problem, characterised by two different sets of Lyapunov exponents associated with physics on different scales. Close encounters between particles trigger microchaos, a generic feature of the N -body problem, which leads to large positive Lyapunov exponents χ N . If, however, the bulk smooth potential is chaotic, one also encounters macrochaos, which is again characterised by positive, albeit typically much smaller, Lyapunov exponents χ S . N -body realisations of integrable systems remain chaotic, even for large N , in the sense that χ N does not decrease towards zero for N → ∞. Despite this, however, microchaos becomes progressively less important macroscopically in that the range of this chaos, i.e., the scale on which the microchaos-driven exponential divergence of nearby orbits terminates, decreases with increasing N .
Section V summarises the principal conclusions and speculates on potential implications. The computations described in this Section were performed for two models which, albeit not representative of 'real' equilibrium systems, are of significant pedagogical value in illustrating the nature of the continuum limit. In particular, since one model is integrable and the other almost completely chaotic, it is simple to identify separately how discreteness effects impact regular versus chaotic orbits, an issue that becomes more difficult in 'realistic' systems which admit a complex coexistence of both regular and chaotic orbits. Model 1. A constant density triaxial ellipsoid, for which with For m ≤ 1, this corresponds to a space-charge potential with Φ 0 a constant and frequencies ω a , ω b , and ω c related to the axis values a, b, and c in terms of incomplete elliptic integrals. The system was assumed confined by an external potential Φ ext = −2Φ sc . It follows that, in the continuum limit, each charge evolves in an integrable time-independent potential of the form (modulo a constant) Attention focused primarily on the parameter values a = 1.95, b = 1.50, and c = 1.05, which, assuming units for which Q ≡ 1, implies [13] that ω a ≈ 0.4663, ω b ≈ 0.5508, and ω c ≈ 0.6753. It follows that the orbital time scale Because this potential is integrable, one knows that, in the continuum limit, only regular phase mixing is possible. However, the situation is even more exceptional: because of the harmonic character of the potential, i.e., the fact that the force is linear, every charge will orbit with the same frequencies so that, in the absence of discreteness effects, there can be no systematic phase mixing and no emittance growth in an initially localised clump. All emittance growth associated with this potential must be attributed to discreteness effects. Model 2. Perturbing Model 1 by introducing a spherically symmetric, attractive spike of charge near the origin, which yields a modified potential with ℓ = 10 −3 . Attention here focused on a central charge q = 10 −1.5 Q ≈ 0.03162, which leads to a potential for which, for orbits restricted energetically to m ≤ 1, the phase space is almost completely chaotic [14]. (The bulk properties of the potential are insensitive to the precise value of ℓ for ℓ < 10 −2 or so; but most of the chaos is lost for much larger values of ℓ.) Frozen-N charge density distributions of the form were generated by randomly sampling the smooth density distributions ρ. These correspond to N -body potentials which incorporate a tiny softening parameter e [15]. Unless otherwise stated, all Figures in this paper were generated from integrations with e = 10 −5 . Orbits were integrated in frozen-N realisations with N ≤ 10 6 using a variable time step scheme that conserved the energy of each charge to at least one part in 10 5 . Estimates of the largest (finite time) Lyapunov exponent [16] were obtained in the usual way by simultaneously tracking the evolution of a small initial perturbation, periodically renormalised at fixed time intervals ∆t.
The efficacy of phase mixing was tested by generating localised clumps of 1600 initial conditions sampling a phase space region of size |∆r| ∼ |∆v| ∼ 10 −3 and evolving these into the future in both the smooth potential and the corresponding N -body potential (7). The resulting orbital data were analysed to compute emittances as well as the total The degree to which individual N -body orbits did, or did not, 'look highly irregular' was quantified by computing the orbital complexity [17] of their Fourier spectra. This entailed determining for each orbit the quantities n x , n y , and n z , defined as the minimum number of frequencies required to capture a fixed fraction k of the power in each direction, and then assigning a total complexity n(k) = n x + n y + n z .
In order to obtain a reasonably sharp Fourier spectrum, orbital data were typically recorded at intervals δt = 0.05, a time corresponding to less than 1% of a typical orbital period.

B. Regular and chaotic phase mixing
In the continuum limit, initially localised clumps characterised by the integrable potential (4) exhibit no systematic tendency to disperse. Because each orbit executes harmonic motions with the same three frequencies, the charges remain close together, returning to their original x 0 , y 0 , and z 0 , after periods τ x , τ y , an τ z . Discreteness effects break this exact periodicity and trigger a systematic spread. This is illustrated in Figure 1, which exhibits the x and y coordinates of the same 1600 orbit ensemble with E = 1.0 at five different times, allowing for frozen-N backgrounds with N = 10 3.5 and N = 10 5 . For N = 10 3.5 , this dispersal is comparatively rapid, the charges having spread to sample the entire allowable configuration space within t = 128, a time corresponding to only 10 orbital periods or so. For the larger system with N = 10 5 , the dispersal is considerably slower, requiring a time t ∼ 512, roughly four times larger, to achieve a comparable spread.
The situation is very different for the potential (5), for which, even in the continuum limit, the particle phase space is almost completely chaotic. In this case, one observes exponentially fast chaotic phase mixing in the smooth potential, and allowing for discreteness effects only accelerates the process. This is evident from Figure  2, which exhibits the x and y coordinates for a 1600 orbit clump, again with E = 1.0, evolved in frozen-N realisations with N = 10 3.5 , N = 10 4.5 , and N = 10 5.5 , as well as (in the right hand column) the smooth potential. Even for the smooth potential, a time t ∼ 128 is sufficient for particles to sample most of the energetically accessible phase space.
The visual impression that the chaotic clump disperses far more rapidly can be quantified by computing the emittance ǫ as a function of time. The left hand column of Figure 3 exhibits ǫ(t) for the same ensemble of initial conditions used to generate Figure 1, now allowing for frozen-N backgrounds with N = 10 3.0 , 10 3.5 , 10 4.0 , 10 4.5 and and 10 5.0 . For the smallest value of N it is not completely clear whether ǫ grows exponentially or as a power law in time. However, for N ≥ 10 3.5 , the growth is distinctly subexponential. Indeed, the data for N ≥ 10 3.5 are well fit by an emittance growth law Extrapolating to the limit N → ∞ yields the expected result that there can be no systematic emittance growth. The left hand column of Figure 4 exhibits analogous results for the initial conditions used to generate Figure 2, now plotted on a logarithmic scale, allowing for N = 10 2.5 , 10 3.5 , 10 4.5 , 10 5.5 , and, in the bottom panel, the smooth potential. It is evident that, for the smooth potential, ln ǫ exhibits a roughly linear growth during the interval (say) 10 < t < 100, corresponding to an exponential growth in emittance. This is hardly surprising. The fact that individual orbits in the clump are chaotic implies that they should diverge exponentially so that, at least for small ǫ(0), one would expect ǫ to grow exponentially at a rate comparable to the value of a typical (finite time) Lyapunov exponent χ S for the smooth potential. For this ensemble, the mean exponent for the interval 0 < t < 256 assumed the value χ S = 0.056, which corresponds to the slope of the dashed line in panel (i).
Allowing for discreteness effects clearly accelerates the rate of chaotic phase mixing. For the two smaller values of N , 10 2.5 and 10 3.5 , the growth again appears exponential, albeit at a larger rate; but for the systems with N = 10 4.5 and 10 5.5 the evolution is clearly more complex. Indeed, a careful examination of the data for these two cases suggests strongly that the evolution can be decomposed into two largely distinct exponential phases, the former characterised by a growth rate ≫ χ S and the latter by a slower rate ∼ χ S . This is consistent with the analysis to be presented later in Section IV, which indicates that two sorts of chaos can act in N -body systems, large scale macrochaos characterised by a Lyapunov exponent χ S and shorter range microchaos characterised by a Lyapunov exponent χ N ≫ χ S . (For very small N , the microchaos also acts on macroscopic scales, thus overwhelming any observational effects associated with the much weaker macrochaos: hence the (near-)absence of the second exponential phase in panels (a) and (c)!) The data summarised in Figure 4 are consistent with a second exponential phase with the form of which will be motivated in Section IId. The results derived here for Model 2 which, in the continuum limit, corresponds to a chaotic potential, are likely generic for bulk density distributions corresponding to a chaotic potential. However, the results for Model 1 are special in that there is no systematic emittance growth in the continuum limit. If, as one would expect in a 'real' system, the bulk potential exhibits at least some anharmonicities, regular phase mixing will trigger linear emittance growth even in the continuum limit. In this case, allowing for discreteness effects will again accelerate the growth of emittance, but the exact form of this enhanced growth can be more complex, even though it will again be subexponential. Section III will exhibit additional examples of how discreteness can accelerate emittance growth for both regular and chaotic clumps.

C. Individual Orbits and the Continuum Limit
The preceding indicates that, as N increases, orbit ensembles evolved in frozen-N backgrounds more closely resemble orbit ensembles with the same set of initial conditions evolved in the corresponding smooth potential. This does not, however, necessarily imply that individual orbits also converge towards characteristics in the smooth potential. To what extent, then, is it true that, as N increases, individual trajectories come to more closely resemble smooth potential characteristics?
The most obvious -and compelling -check is visual: do frozen-N orbits 'look like' smooth potential characteristics when N becomes sufficiently large? As a simple, and extreme, example, consider a constant density spherical system without a central spike where, in the continuum limit, Φ reg reduces to an isotropic harmonic oscillator potential; and select an initial condition which, in the smooth potential, corresponds to a circular orbit. Results derived from integrations of such an initial condition in different frozen-N systems are exhibited in Figure 5, which shows representative frozen-N orbits generated for particle number varying between N = 10 2.5 and N = 10 5.5 , along with the smooth potential orbit. For the four smallest values of N , there is no obvious hint that the orbit 'should' be circular or that there 'should' be a net sense of circulation, although there is a crude visual sense that, as N increases, the orbit becomes 'less tangled.' However, for N = 10 4.5 one starts to discern a clear sense of net circulation, for N = 10 5.0 the orbit has clearly become centrophobic (thus suggesting that angular momentum is at least approximately conserved), and for N = 10 5.5 the orbit arguably resembles a 'distorted' or precessing circular orbit.
As is illustrated in Figure 6, the visual impression that the orbit is becoming more nearly circular can be corroborated by constructing the Fourier spectra of the orbital data. For the three smallest values of N , the quantity |x(ω)| is obviously broad band, although there is a peak at or near the circular frequency associated with the smooth potential orbit. For N = 10 4.0 and 10 4.5 the peak becomes appreciably sharper, and for N = 10 5.0 and 10 5.5 one sees only slight irregularities in the spectra which translate into the visual appearance of precession.
The conclusion is obvious: as N increases, frozen-N orbits come to more closely approximate the smooth orbit, both visually and in terms of their power spectrum. Analogous results obtain for more generic initial conditions evolved in this and other integrable potentials.
Convergence of orbits in terms of their Fourier spectra is important in justifying straightforward applications of nonlinear dynamics to many-body systems interacting via long range forces. Many physical phenomena in many-body systems, including accelerator modes [18], modulational diffusion [19], and resonant relaxation [20], are attributed to resonant couplings between, e.g., the natural frequencies of individual regular orbits and the frequency or frequencies associated with some perturbation. However, such applications can only be justified if the 'real' N -body orbits have frequencies that adequately approximate the frequencies associated with characteristics in the smooth potential.
The degree of irregularity exhibited by individual orbits can be quantified by computing Fourier complexity, as defined in Eq. (10). The results of one such investigation are summarised by the curves with diamonds in Figure 7, which exhibit n(0.95), the mean number of frequencies required to capture 95% of the total power, for collections of 100 initial conditions evolved in frozen-N backgrounds with different N . (The triangles will be discussed in the following subsection.) In each case, the initial conditions were integrated for a time T = 128, with orbital data recorded at fixed intervals δt = 0.05. The data were then Fourier analysed using an FFT solver to translate a set of 2j points into a set of j Fourier amplitudes. The upper panel was generated for orbits in the integrable Model 1, the lower panel for the strongly chaotic Model 2. In each panel, the solid curve represents the mean complexity computed in the unperturbed smooth potential. Error bars were computed by dividing the 100 initial conditions in half and analysing each half separately.
It is evident that, in both cases, n(0.95) is a decreas-ing function of N which converges towards the continuum value for N → ∞. For smaller values of N , the regular and chaotic ensembles have comparable complexities, although the regular ensemble is slightly less (∼ 20%) complex. However, for larger values of N there are clear distinctions between the regular and chaotic ensembles, the chaotic ensembles for N ≥ 10 5 being nearly twice as complex as the regular ensembles. Indeed, for other choice of regular models the complexity can be even lower: The fact that the smooth potential complexity in panel (a) is as large as it is reflects the fact that the regular orbits used to generate this Figure were relatively complex 'box' orbits, with the topology of Lissajous figures, which, even in the continuum limit, require two or three frequencies in each direction to capture as much as 95% of the power. If instead n(0.95) is computed for an ensemble of initial conditions corresponding to circular orbits, the complexity converges towards a continuum limit with n(0.95) = 3.

D. Modeling N -body orbits and flows by Gaussian white noise
Conventional wisdom holds that discreteness effects can be idealised as friction and Gaussian (nearly) white noise in the context of a Fokker-Planck description [2]. Taken literally, this suggests that individual N -body orbits can be well-mimicked by Langevin simulations. However, it is not completely clear to what extent this is really true. The original derivation of the Fokker-Planck equation (and most if not all of its tests) restricted attention to the statistical properties of orbit ensembles or distribution functions over comparatively long time scales, assuming implicitly that the bulk potential in which the particle evolves is nonchaotic. Can the friction/noise paradigm describe correctly short time behaviour and/or the behaviour of individual orbits, especially in a chaotic potential?
Granted the validity of a Fokker-Planck description, a simple rule connects N to the strength of the friction and noise. Assuming that the noise is characterised by a temperature per unit mass Θ comparable to the magnitude of the particle energy, the coefficient of dynamical friction η defines an energy relaxation time t R = η −1 . However, an evaluation of the Fokker-Planck coefficients in a binary encounter approximation leads to the prediction that t R ∝ (N/ ln Λ)t D , with ln Λ the Coulomb logarithm. Given the assumption of a nonneutral plasma, the treatment of Λ must necessarily be somewhat heuristic [21]. However, there is a general agreement that Λ should scale as some power of N , so that t R , and hence the diffusion constant D, should satisfy The obvious question, then, is whether frozen-N simula-tions with specified N can be well-mimicked by Langevin simulations with η ∝ (ln N/N ). Two practical issues arise in testing this prediction. The first is that, because of the limited range of N that can be explored, it is impractical to test the subdominant ln N dependence: for N < 10 3 or so, the very notion of a bulk potential fails; for N > 10 6 computations become prohibitively expensive. One must instead restrict attention to testing the simpler scaling relation η ∝ N −1 , i.e., for some constant p.
The second point is more serious. The usual Langevin equation reads [4] where ηdr a /dt represents a dynamical friction. F a represents Gaussian white noise, which is characterised completely by its first two moments: and with D ≡ 2ηΘ the diffusion constant entering into a Fokker-Planck description. By choosing Θ to equal the initial energy one can ensure that the average energy of the orbits remains unchanged. Such an equation is clearly unsatisfactory here. Energy is conserved absolutely for frozen-N orbits, so that one must also impose energy conservation on any scheme which aims to mimic its effects. (For very small η, energy remains almost conserved for very long times. However, comparatively small N should correspond to relatively large η, which implies large changes in energy and, as such, significant changes in the phase space regions accessible to the noisy orbit.) For this reason, the noisy integrations described here were performed using a modified energy-conserving noise [22].
This entailed (1) eliminating the dynamical friction altogether, (2) again imparting random kicks as in Eq. (17), but (3) renormalising the modified velocity at each time step by an overall factor, i.e., v(t + δt) → αv(t + δt), with α so chosen that E(t + δt) = E(t). Modulo this complication, the noise was integrated using a standard algorithm [23] based on a fourth order Runge-Kutta integration scheme with fixed time step δt. The integrations were performed for δt = 2 × 10 −4 , it having been confirmed that the statistical effects of the noise were insensitive to the precise value of δt for δt < 10 −3 .
At the level of orbit ensembles, as probed by the emittance and other bulk moments, the results of frozen-N simulations are in fact extremely well-mimicked by Langevin simulations, at least for comparatively large N . The degree to which this is true can be inferred by contrasting the right and left hand columns of Figures 3 and 4. As discussed already, the left hand columns of Figures  3 and 4 exhibit, respectively, time-dependent emittances for Models 1 and 2, allowing for several different values of N . The right hand panels exhibit data generated from Langevin integrations of the same initial conditions, allowing for amplitudes η satisfying Eq. (15) with p = 0.5, so that, e.g., N = 10 5.5 corresponds to η = 10 −5.0 . For the smallest values of N (and hence the largest values of η) -corresponding to panels (a) and (b) in Figure 3 and panels (a) -(d) in Figure 4, the agreement is not all that good. However, for larger particle number -N ≥ 10 3.5 for the regular system and N ≥ 10 4.5 for the chaotic system -, the agreement is obviously quite good.
The bottom right hand panel in Fig. 4 was generated for orbits evolved with a considerably smaller value of η, namely η = 10 −7.5 , this corresponding to the largest noise amplitude that does not alter appreciably the emittance growth observed in the smooth potential. To the extent that the scaling of Eq. (15) is in fact correct for p ≈ 0.5, the fact that larger values of η have an appreciable effect on emittance growth implies that, even over an interval as short as t = 128, corresponding to ∼ 10 orbital time scales t D , one requires N > 10 7 to justify a continuum limit! Even though the collisional relaxation time scale t R ∝ (N/ ln N )t D ≫ t D , discreteness effects can be important in a system with N ∼ 10 6.5 on a time scale as short as ∼ 10t D .
As is evident from Figure 7, this agreement also extends to the level of individual orbits. As described already, the diamond curves in this Figure were derived from frozen-N integrations. The other curves, constituted of triangles, were generated from exactly the same initial conditions, now integrated, however, in the smooth potential while being subjected to energyconserving Gaussian noise with Θ = E and variable η. To the extent that conventional Fokker-Planck theory is correct, one would anticipate a correspondence between N and η of the form given by Eq. (15). The noisy points in Figure 7 were in fact identified with the frozen-N points assuming the validity of the scaling (15) with p = 0.6. The obvious fact, then, is that, given this identification, the curves n(N ) and n(η) rather nearly coincide. Even at the level of the complexity of individual orbits, frozen-N orbits can be well-mimicked by noisy orbits with ln η + ln N = const.
Granted that discreteness effects can be mimicked by Gaussian white noise, the scaling relations (11) and (13) are easily understood. At least for a harmonic potential, it is straightforward to derive analytic solutions to the Langevin equation (16) for moments like x 2 or xv x [3]. Alternatively, it is easily seen that the Fokker-Planck equation associated with Eq. (16) implies that the clump emittance satisfies Assuming, however, that the initial emittance is extremely small, at early times one can approximate that x 2 v 2 ≈ xv 2 ; and, to the extent that the growth time is long compared with the characteristic crossing time, one can average over oscillations to set x 2 = E/ω 2 , with E the initial energy. It then follows that, for early times, Combining Eq. (19) and the analogous formulae for ǫ y and ǫ z with the scaling relation η ∝ 1/N leads immediately to Eq. (11). The same diffusive t 1/2 behaviour also arises for δr rms and δv rms . A somewhat more heuristic argument can account for the scaling (13) associated with a chaotic potential. If the initial emittance ǫ(0) = 0, it is clear that, in the absence of discreteness effects, ǫ(t) would continue to vanish identically: two smooth integrations of the same initial condition will not yield divergent orbits, even if the orbits are chaotic. However, discreteness effects act to 'kick' two nearly coincident orbits apart, at which point they will tend to diverge at a rate set by the Lyapunov exponent χ S associated with the bulk potential. Assuming, however, that the 'kicks' are random, their effects will scale as η 1/2 rather than η; but combining this with Eq. (15) implies that [24] δr rms ∝ δv rms ∝ N −1/2 exp(χ S t).

III. THE THERMAL EQUILIBRIUM MODEL
A. Defining the model Consider now a more realistic example, the thermal equilibrium model [10], which, in the continuum limit, admits large measures of both regular and chaotic orbits [11]. This model allows for a collection of N identical charged particles, interacting electrostatically, that is constrained by linear restoring forces to manifest triaxial symmetry, the focusing forces in different orthogonal directions being characterised in general by unequal frequencies. Individual particles thus have energy where x and v denote particle position and velocity, m and q denote the mass and charge, ω = (ω x , ω y , ω z ) represents the three frequencies associated with the focusing force, and φ(x) is the collective space-charge potential.
The additional assumption is that the particles can be characterised by a one-particle distribution function of the Maxwell-Boltzmann form, f ∝ exp(−E/k B T ), with k B T the temperature. This implies a bulk number density satisfying where φ(x) is defined implicitly as a function of n via the relations (in mks units) Following, e.g., [11], the problem can be cast into dimensionless form by expressing length and frequency in units of the Debye length and plasma frequency, i.e., and by introducing a dimensionless potential With appropriate rescaling, the net result is a density distribution of the form where Here Ω 2 = (ω y /ω p ) 2 and R 2 = (x/a) 2 + y 2 + (z/c) 2 , in terms of scale lengths a and c satisfying a = ω y /ω x and c = ω y /ω z . The minimum permissible focusing strength corresponds to Ω = Ω u = 1/ (1/a 2 ) + 1 + (1/c 2 ).
The experiments described here were performed assuming parameter values a 2 = 0.5, c 2 = 1.5, and Ω = 1.0001/ √ 3, for which a typical orbital time scale t D ∼ 20. These parameters represent a beam that is moderately, but not strongly, dependent on space charge. Consider, for example, a proton bunch with 1 µm root mean squared normalised emittance spanning 3 cm full 'radius'. If the bunch is described by the thermal equilibrium model, the Debye length is ∼ 2 mm and the bunch population is the ∼ 3 × 10 9 protons, this corresponding to a bunch charge ∼ 0.5 nC.
In general it does not appear possible to solve eqs. (26) and (27) analytically. This makes both the generation of N -body realisations of the density and the computation of orbits in the smooth potential much more difficult. However, these difficulties can be, and were, resolved using numerical techniques described in [11]. In principle, the accelerations for the N-body thermal equilibrium model should be given bȳ (29) with N the number of frozen particles and N m satisfying In practice, however, one cannot perform this integral, even numerically, since Φ was only evaluated on a finite grid. For this reason, the integral was first solved with limits coinciding with the grid boundaries, but then renormalised by a small constant factor so that plots of the potential and density in the smooth and frozen-N configurations overlapped perfectly.

B. Regular and chaotic phase mixing
The top four left hand panels of Fig. (8) exhibit emittance growth for an initially localised orbit ensemble evolved in frozen-N realisations of the thermal equilibrium model, selected with energy sufficiently small that the constant energy hypersurface in the smooth potential is completely regular. Since the size of the accessible phase space is roughly ten times larger than was the case for the oscillator models, the initial conditions sampled a region ten times larger, |∆r| ∼ |∆v| ∼ 10 −2 . It is evident that, as for the integrable oscillator model considered in Section II, the emittance growth is power law rather than exponential; and, at least for N = 10 5.5 and N = 10 6.0 , it is well fitted by a growth law ǫ ∝ t 1/2 . The bottom left panel exhibits emittance growth for the same orbit ensemble evolved in the smooth potential. Here the evolution is clearly linear, rather than square root, the expected behaviour for smooth orbits in generic integrable potentials where nearby initial conditions correspond to slightly different orbital frequencies. (Recall that, for the oscillator Model 1, there is zero emittance growth in the continuum limit.) That ǫ grows much faster for the frozen-N model with N = 10 6.0 than for orbits in the smooth potential indicates clearly that, for the thermal equilibrium model, N = 10 6.0 is not a sufficiently large particle number to justify a continuum approximation, even over an interval as short as t = 512.
The top four right hand panels demonstrate that discreteness effects can again be well-mimicked by energyconserving Gaussian white noise with Θ = E and N and η related as in Eq. (15) although, in this case, the best fit value p ≈ 1.5, rather than p ≈ 0.5. The bottom right panel exhibits the emittance growth for an ensemble evolved with η = 10 −6.5 , the largest noise amplitude that does not significantly alter emittance growth in the smooth potential. Presuming that the scaling (15) can again be extended to larger N and smaller η, one infers that, for the purpose of predicting emittance growth in this regular ensemble, the smallest value of N for which the continuum limit can be justified is N ∼ 10 8 . Figure 9 exhibits analogous data for a higher energy ensemble which, in the continuum limit, corresponds completely to chaotic orbits. It is evident that, as for the chaotic model in Section II, the evolution is exponential overall, rather than power law; and that discreteness effects again have an important effect. Also evident from a comparison of left and right hand panels is that, as for regular orbits, discreteness effects can be well-mimicked by noise with N ∝ 1/η and p ≈ 1.5. Most striking, however, is the fact that, in this case, even much weaker noise can accelerate emittance growth appreciably. For this chaotic ensemble, discreteness effects must correspond to a noise amplitude satisfying η < 10 −8.0 or so before a continuum limit can be justified. Chaotic orbits are far more susceptible to low amplitude noise than are regular orbits. Presuming again that the scaling relation (15) holds, this implies that for the case of chaotic orbit ensembles, the continuum limit cannot be justified for N < 10 9.5 .
This, coincidentally, is roughly the number of particles in the equilibrium proton bunch described in Sec. III A. Accordingly, in studying the dynamics of beams with moderate space charge, one may not be able to assume the validity of the continuum limit with complete confidence, even for a system in equilibrium. The situation may be even more problematic for a beam that is significantly out of equilibrium, since the resulting timedependent potential would be expected to generate a larger population of chaotic orbits [25].

C. Transitions between regular and chaotic behaviour
At very low energies, where the total potential is nearly harmonic, all smooth potential orbits are regular, so that discreteness effects can only act to deflect frozen-N orbits from one regular trajectory to another. However, at higher energies the smooth potential admits a complex coexistence of regular and chaotic orbits. This implies the possibility that discreteness effects can deflect frozen-N orbits from regular to chaotic characteristics and vice versa [12]. Of obvious interest then is how fast, as a function of N , such transitions occur.
For example, an accelerator designer relying on the Vlasov equation and an analysis based on a smooth, macroscopic potential would neglect these microscopic transitions and their corresponding impact on chaotic mixing. Thus the physics of collective relaxation and global emittance growth would be improperly modeled and the results, at least in principle, would be suspect. The consequences of such an omission depend on the problem at hand, but one might expect them to be especially severe for nonequilibrium beams where chaotic dynamics is likely to be more prevalent [25].
If one selects a localised ensemble of initial conditions corresponding to regular orbits in the smooth potential and integrates these initial conditions into the future, discreteness effects will, if sufficiently strong, eventually trigger transitions from regularity to chaos. That such transitions actually occur can be determined by a visual inspection of individual frozen-N orbits which can be observed to become abruptly 'more irregular' in appearance. If, moreover, large numbers of transitions occur over very short times, this can make the emittance associated with an initially localised ensemble, which ought to grow as a power law in time, exhibit instead a more rapid, roughly exponential, increase.
However, an accurate determination of the relative fraction of the orbits which are still regular requires a more careful analysis. This was done by recording the phase space coordinates of individual frozen-N orbits at various times t > 0, and evolving these into the future in the smooth potential to determine whether the resulting smooth potential characteristics were still regular or whether instead they had become chaotic. The most straightforward fashion in which to determine whether the orbits are chaotic would have been to compute an estimate of the largest (finite time) Lyapunov exponent. Given, however, that the potential cannot be expressed in terms of elementary functions, this would have proven extremely expensive computationally. For that reason, distinctions between regularity and chaos were based instead on the computed complexities of the characteristics. As discussed elsewhere (e.g., [11], [17], and references cited therein), such a criterion typically coincides almost exactly with more conventional criteria based on Lyapunov exponents. Presuming that the system is ergodic and that discreteness effects are sufficiently strong that they can in principle convert any orbit from regular to chaotic, and vice versa, it would seem clear what such an analysis ought to reveal. (1) At sufficiently late times, independent of N the relative fraction of chaotic orbits generated from any initial ensemble should (to within statistical uncertainties) coincide with the relative measure of chaotic orbits on the constant energy hypersurface, i.e., to the relative volume of the chaotic portions of the constant energy hypersurface. (2) Assuming, however, that discreteness effects are more important for smaller N , the time required to converge towards this asymptotic value should be an increasing function of N . As N increases, transitions should become more rare.
As illustrated in Fig. 10, which summarises computations with particle number between 10 4.5 and 10 6.0 , this expectation was in fact confirmed. For frozen-N systems with number as small as N = 10 4.5 , nearly 40% of the orbits in an initially regular ensemble had become chaotic within a time t = 64, a time corresponding to only ∼ 3t D , and the relative measure f of chaotic orbits appears to have asymptoted towards a time-independent value by t = 128. The relative measure of chaotic orbits for computations with N = 10 5.0 grew more slowly in time; but, by t = 512, the relative measure had again approached a value f ∼ 0.4. For larger values of N , f remains a monotonically increasing function of time, but transitions are sufficiently rare that one does not approach an equilibrium population within a time as short as t = 512 ∼ 25t D .

A. Ordinary Lyapunov exponents
As N increases, frozen-N orbits come to more closely resemble smooth potential characteristics generated from the same initial condition, both visually and in terms of their Fourier spectra. One might therefore expect that, at least for a regular, i.e., nonchaotic, smooth potential, the value of the largest Lyapunov exponent χ N should decrease with increasing N and converge towards zero for N → ∞. Such, however, is not the case. Rather, as is also true for gravitationally interacting systems of particles [22], for both regular and chaotic orbits the value of χ N is comparatively insensitive to N ; and there are even indications that χ N might increase with increasing N . Figure 11 exhibits the value of the largest Lyapunov exponent χ N as a function of N for a single initial condition evolved in frozen-N integrations with softening parameters varying between e = 10 −5 and e = 10 −1 . Figure 12 exhibits the same data, now plotting χ N as a function of e for different values of N . In both Figures, each point was generated by integrating the same initial condition used to generate Figures 5 and 6 for a total time t = 256 in 10 different frozen-N realisations of Model 1. In the continuum limit this initial condition corresponds to an integrable circular orbit with vanishing Lyapunov exponents; and, as was evident from Figure 5, the frozen-N orbits for larger N look much more regular in appearance than do the orbits with smaller N . Despite this, however, at least for the smallest values of e, χ N does not decrease with increasing N . As probed by χ N , orbits with N = 10 2.5 and 10 5 are comparably chaotic! However, for larger values of e, χ N does decrease with increasing N . That this should be the case is easily understood. Because the bulk potential is integrable, the chaos must at some level be associated with close encounters between nearby charges; but introducing a softening parameter de facto 'turns off' encounters on scales shorter than e. If the charge density is sufficiently large that encounters with separation < e become common, the simulation will have artificially, and erroneously, reduced this source of chaos, yielding a smaller χ N .
That the value of χ N for nearly unsoftened frozen-N integrations is insensitive to whether the smooth potential is regular or chaotic is illustrated in Figure 13, which exhibits χ N as a function of N for integrations with e = 10 −5 for the same initial condition integrated in both Models 1 and 2. That Figure also emphasises another important point, namely that χ N is typically larger than any Lyapunov exponent χ S associated with motion in the smooth potential by an order of magnitude or more. For this particular initial condition, χ N ≈ 0.82 whereas χ S ≈ 0.056.
That χ N should be insensitive to the choice of N , at least for unsoftened simulations, might seem rather curious. However, it is not hard to understand why this might be so for a 1/r 2 force. Given that the microchaos disappears completely in the continuum limit, it would seem clear that it must be associated with a sequence of 'random' interactions between a 'test' charge and a collection of 'field' charges. However, this suggests that the Lyapunov time t * ≡ χ −1 N associated with the growth of a small initial perturbation can be estimated by considering tidal effects associated with a pair of charges separated by a distance r comparable to (some fixed fraction of) the typical interparticle separation r sep . This tidal acceleration scales as with q the magnitude of an individual charge. Given, however, that r sep ∼ n −1/3 ∼ N −1/3 R sys , with R sys the size of the system and n a characteristic number density, and that q = Q/N , with Q the total charge, it follows that the time scale t * , and hence χ N , should be independent of N . As N increases, the sizes of the individual charges and the cube of the typical interparticle separation both decrease as N −1 , so that the ratio is independent of particle number. A more careful argument [26] actually allows one to prove analytically that, in the absence of softening, χ N cannot decrease towards zero with increasing N . For simple geometries, an analytic expression for the average value of the stability matrix entering into the definition of χ N can be formulated in terms of a 3N -dimensional integral. This integral cannot be evaluated analytically, but one can derive rigorous bounds which ensure that the largest eigenvalue remains positive even for N → ∞; and, even more strikingly, Monte Carlo evaluations of the integrals suggests that χ N should be a slowly increasing function of N . In other words, viewed in terms of χ N , orbits become more chaotic as N increases, even though, for the case of a regular potential, they become more regular in appearance! [27] The obvious inference is that N -body Lyapunov exponents χ N do not provide a useful characterisation of the degree of chaos associated with an orbit, at least when that orbit is viewed macroscopically.

B. Microchaos and macrochaos in the N -body problem
That frozen-N orbits have large positive Lyapunov exponents χ N , even for the case of an integrable potential, but that distinctions between regular and chaotic potentials are clearly manifested in the phase mixing of initially localised clumps would suggest that there are two different, and comparatively distinct, potential sources of chaos in the N -body problem.
On the one hand, one would expect microchaos, triggered by close encounters between individual charges, which is manifested only on very short scales, comparable to, or perhaps somewhat larger than, the typical interparticle spacing. This source of chaos should be generic to the N -body problem, irrespective of the form of the bulk potential, generating randomness qualitatively similar to what arises in pinball. On the other hand, there is also the possibility of larger scale macrochaos, which would be expected if and only if the bulk potential admits global stochasticity.
If these expectations are in fact correct, two nearby initial conditions evolved in a frozen-N realisation of any potential should diverge exponentially at a rate ∼ χ N until their separation becomes somewhat larger than a typical interparticle spacing, at which point the microchaos would 'turn off.' If the bulk potential is regular, no other source of chaos could act and the two orbits would continue to diverge as a more modest power law. If, however, the bulk potential is chaotic, macrochaos would still act, resulting in a continued exponential divergence, albeit at a rate ∼ χ S typically much smaller than χ N . In this case, exponential divergence should be replaced by a power law divergence only once the separation has become macroscopic.
This three-stage evolution for chaotic orbits is clearly illustrated in Figure 14, which exhibits data for frozen-N simulations with N = 10 5.0 , 10 5.5 , and 10 6.0 . The three curves in that figure were each generated by selecting 50 fiducial initial conditions in a phase space region of size 5 × 10 −4 and 50 perturbed initial conditions that were subjected to a displacement δx = 10 −5 , evolving each initial condition into the future and computing as a function of time the phase space separation [29] δZ = |δr| 2 + |δv| 2 . (32) That δZ experiences two distinct stages of exponential evolution, at very different rates is especially evident in the curve with N = 10 6.0 . The solid lines accompanying that curve have slopes 0.82 and 0.056, corresponding, respectively, to the mean values of the N -body χ N and the smooth potential χ S for those initial conditions.
As N increases, the initial exponential phase terminates for smaller values of δZ. The microchaos responsible for this first phase will 'turn off' when |δr| becomes large compared with a typical interparticle spacing, but that interparticle spacing scales as N −1/3 . Figure 15 shows the analogue of Figure 14, now generated for regular initial conditions in the thermal equilibrium model. Once again there is an initial exponential growth at a rate comparable to χ N , but in this case there is no evidence of a second exponential phase. Rather, the initial exponential phase is followed immediately by an interval of power law divergence.
One other point, not obvious from these Figures, is that the scaling with N observed for the final power law phase is different for the regular and chaotic systems. For regular orbits, the phase space separation δZ, like |δr| and |δv|, satisfies a linear growth law For chaotic orbits, one finds that δZ again grows linearly in time, but that the growth time

C. Alternative interpretations of Lyapunov exponents
The standard definition of Lyapunov exponents implies that they probe the average rate of divergence for two nearby chaotic orbits in a single system. However, for the N -body problem, Lyapunov exponents also quantify two other effects which, as a practical matter, are of equal importance: 1. Lyapunov exponents probe the rate at which orbits generated from the same initial condition but evolved in two different frozen-N systems diverge from one another. 2. Lyapunov exponents probe the rate of divergence associated with orbits evolved from the same initial condition in both a frozen-N system and in the smooth potential.
This means that Lyapunov exponents also provide information about the degree to which characteristics generated in the smooth potential can be interpreted as providing a pointwise approximation to real frozen-N orbits with the same initial condition, as well as the degree to which orbits in two different frozen-N systems remain close to one another in a pointwise sense.
This is important at a practical level. In a real experiment one may perhaps be able to ensure that a given N -body system constitutes (nearly) a fair sampling of some specified density distribution, but the details of the actual N -body distribution are likely inaccessible. Of obvious interest, therefore, are the questions: to what extent will orbits in two different N -body systems coincide? and to what extent do such orbits coincide with characteristics in the bulk potential associated with the smooth density distribution? Figure 16 exhibits the analogue of Figure 14, now generated by comparing the same 50 initial conditions evolved in two different frozen-N systems. Figure 17 compares orbits in a frozen-N system with orbits in the smooth potential. In each case, the duration of the initial interval of especially fast exponential divergence is significantly reduced, but the second interval with divergence at a rate ∼ χ S is still conspicuous.
It is not hard to understand why the smooth potential χ S provides information about orbits in different frozen-N simulations and/or their relation to orbits in the smooth potential. As noted already, discreteness effects can be extremely well-mimicked by noise, at least mesoscopically. However, after the rapid decay of any initial transients, multiple noisy realisations of the same initial condition corresponding to a chaotic orbit typically diverge exponentially in such a fashion that [24] δZ where the second proportionality follows from the observed scaling η ∝ 1/N .

V. CONCLUSIONS
Viewed macroscopically, there is a precise sense in which, as N increases, trajectories in frozen-N systems converge towards characteristics in the corresponding smooth potential. For very small particle number, N < 10 4 or so, the notion of an average bulk potential fails and orbits in frozen-N systems are very different from smooth potential characteristics. In particular, the usual distinctions between regularity and chaos that exist in a smooth potential seem completely lost. [28] However, for larger N one begins to observe clear distinctions between orbits evolved from initial conditions which, in the continuum limit, correspond to regular versus chaotic orbits.
In particular, although discreteness effects cannot be neglected, phase mixing of initially localised orbit ensembles in a frozen-N environment allows for clear distinctions between 'regular' and 'chaotic' clumps. Just as for clumps evolved in a smooth potential, emittance growth for a regular frozen-N clump proceeds as a power law in time, whereas it is roughly exponential for a chaotic clump. However, in both cases the growth is more rapid than in the smooth potential. Discreteness effects accelerate emittance growth for both regular and chaotic clumps.
In terms of both the statistics of orbit ensembles and the complexity of individual orbits, discreteness effects can be extremely well-mimicked by Gaussian white noise in the context of a Fokker-Planck, or Langevin, description, with a coefficient of dynamical friction η and a diffusion constant D consistent with the predicted scaling D ∝ η ∝ (ln Λ)/N , with ln Λ the Coulomb logarithm. A Fokker-Planck/Langevin description appears justified even when considering the short time behaviour of individual orbits. This suggests strongly that Langevin simulations can be used to assess the importance of discreteness effects in systems where N is too large to allow honest direct summation integrations.
To the extent that such an extrapolation is justified, one concludes that discreteness effects can remain important even for very large N , especially for the case of chaotic orbits. Consider, e.g., the role of discreteness effects in accelerating emittance growth for an initially localised clump. For the case of the thermal equilibrium model, a relatively benign system without particularly large density contrasts and without internal substructures, one needs N = 10 8 or more to justify the continuum approximation in tracking the evolution of a regular clump confined initially to a region ∼ 10 −3 the size of the accessible phase space. For the case of a chaotic clump of comparable size, one needs at least N = 10 9.5 . This has obvious implications for beams in that it affects macroscopic mixing and associated changes in the overall phase-space volume. Discreteness effects are also important because they can trigger transitions between regular and chaotic behaviour, a potentially serious problem for charged particle beams. One might, e.g.,, try to initialise a bunch in such a fashion that, although the bulk potential admits chaotic orbits, only regular regions are populated, thus aiming to facilitate emittance compensation. The problem, however, is that discreteness effects could transform significant numbers of orbits from regular to chaotic, thus making compensation far more difficult.
The time scale associated with transitions between regularity and chaos increases with increasing N , such transitions being impossible in the continuum limit; but for any finite N there is presumably a maximum time over which it is safe to ignore these transitions. The critical point, then, as is evident from Fig. 10, is that that time can be much shorter than the collisional relaxation time t R . To the extent that discreteness effects in the thermal equilibrium model can be mimicked by Gaussian white noise, particle number N = 10 6 corresponds to η ∼ 10 −4.5 , which in turn implies a relaxation time t R ∼ η −1 ∼ 10 4.5 . It is, however, evident that, for a N = 10 6 realisation of the thermal equilibrium model, transitions from regularity to chaos can be important already within a time t < 10 2.5 or so! By contrast, the dynals time scale t D ∼ 20.
It should also be stressed that, even if discreteness effects are too weak to facilitate frequent transitions between regularity and chaos, they could well play an important role in accelerating diffusion through a complex chaotic phase space. Generic smooth potentials admitting both regular and chaotic orbits have chaotic phase space regions partitioned by complex structures associated with cantori in two dimensions and the Arnold web in three which, albeit not acting as absolute obstructions, serve as 'entropy' barriers to slow phase space transport [30]. The important point, then, is that even very low amplitude Gaussian white noise can dramatically accelerate diffusion through such barriers [31]. To the extent that discreteness effects can be modeled as Gaussian white noise, they too should act as a significant source of accelerated phase space transport.
The meaning of 'chaos' in the N -body problem is necessarily somewhat subtle. In particular, it is important to recognise that two distinct sources of chaos can exist, associated with physics on different scales. Short range microchaos, associated with close encounters between individual charges, is a generic feature of the Nbody problem, independent of the form of the bulk potential. However, there is also the possibility of larger scale macrochaos which arises if and only if, in the continuum limit, the bulk potential admits chaos. The important point, then, is that these two distinct sources of chaos can be characterised separately by different sets of Lyapunov exponents. Close encounters trigger an exponential separation of nearby trajectories at a rate χ N . The bulk potential triggers an exponential separation at a rate χ S which is typically much smaller.
Standard numerical computations of Lyapunov exponents yield estimates of the much larger χ N , a quantity that does not decrease with increasing N . This leads to the seemingly oxymoronic conclusion that the N -body problem remains strongly chaotic for very large N , even if the potential is integrable in the N → ∞ limit and even if orbits 'look' nearly regular and have Fourier spectra that are nearly periodic. The crucial point, however, is that even though microchaos remains strong in the sense that χ N does not decrease with increasing N , it becomes progressively less important macroscopically. The scale on which the exponential divergence saturates is comparable to a typical interparticle separation r sep , a distance that decreases as N −1/3 with increasing N . By tracking the divergence of nearby orbits, starting with initial separations ≪ r sep and continuing until the separation becomes macroscopic, it is possible to extract estimates of both χ N and χ S .
Finally, it should be noted that, as applied to the Nbody problem, the smooth potential Lyapunov exponent χ S does not simply quantify the average divergence of two nearby trajectories in a single frozen-N simulation. It also quantifies the rate at which a frozen-N trajectory will diverge from a smooth potential characteristic with the same initial condition and the rate at which orbits with the same initial condition diverge in different frozen-N simulations, two quantities which, in some settings, could be even more important physically.   , defined as the mean number of frequencies required to capture 95% of the total power, computed for a collection of 100 initial conditions integrated in frozen-N realisations of the integrable Model 1 with variable N . Triangles show the same quantity for the same initial conditions integrated in the smooth potential but subjected to Gaussian white noise with variable η The solid line exhibits the mean complexity for orbits evolved in the unperturbed smooth potential. (b) The same as (a), generated from the same initial conditions but now computed for the strongly chaotic Model 2. FIG. 17: The mean phase space separation |δZ(t)| for 50 initial conditions evolved in Model 2, both in the smooth potential and in a frozen-N system with (from top to bottom) N = 10 5.0 , 10 5.5 , and 10 6.0 . The solid line again has a slope χS = 0.056, equal to the mean value of the largest smooth potential Lyapunov exponent.