Chaotic Orbits in Thermal-Equilibrium Beams: Existence and Dynamical Implications

Phase mixing of chaotic orbits exponentially distributes these orbits through their accessible phase space. This phenomenon, commonly called ``chaotic mixing'', stands in marked contrast to phase mixing of regular orbits which proceeds as a power law in time. It is operationally irreversible; hence, its associated e-folding time scale sets a condition on any process envisioned for emittance compensation. A key question is whether beams can support chaotic orbits, and if so, under what conditions? We numerically investigate the parameter space of three-dimensional thermal-equilibrium beams with space charge, confined by linear external focusing forces, to determine whether the associated potentials support chaotic orbits. We find that a large subset of the parameter space does support chaos and, in turn, chaotic mixing. Details and implications are enumerated.


I. INTRODUCTION
Rapid, inherently irreversible dynamics is a practical concern in producing high-brightness charged-particle beams. Time scales of irreversible processes place constraints on methods for compensating against degradation of beam quality caused by, for example, space charge. This is a very important practical matter because compensation must be fast compared to these processes, and this affects the choice and configuration of the associated hardware.
A beam bunch with space charge comprises an N -body system with typically 3N degrees of freedom. Upon coarse-graining, i.e., "smoothing" the system to remove granularity, the collective space-charge force remains. One might conjecture that this force, when nonlinear, may support chaotic orbits. One example is the University of Maryland five-beamlet experiment that shows presumably irreversible dissipation of the beamlets after a few space-chargedepressed betatron periods [1]. Simulations of the experiment reveal a substantial fraction of globally chaotic orbits [2], and phase mixing of these orbits thereby presents itself as a contributing evolutionary mechanism. This example pertains to a strongly time-dependent nonequilibrium system, yet one might conjecture that nonlinear space-charge forces in a static system could support chaotic orbits as well. We shall explore this conjecture.
An initially localized clump of chaotic orbits will, via phase mixing, grow exponentially and eventually reach an invariant distribution. This is "chaotic mixing" [3,4]. Strictly speaking, the process is reversible in that it is collisionless and its dynamics is included in, e.g., Vlasov's equation. Nonetheless, when the invariant distribution spans a global region of the system's phase space, chaotic mixing is a legitimate relaxation mechanism in that it drastically smears correlations. Moreover, from a practical perspective, the process is strictly irreversible because infinitesimal fine-tuning is needed to reassemble the initial conditions. It is also distinctly different from phase mixing of regular orbits, i.e., linear Landau damping [5], a process that winds an initially localized clump into a filament over a comparatively narrow region of phase space. Whereas chaotic mixing proceeds exponentially over a well-defined time scale and can cause global, macroscopic changes in the system, phase mixing of regular orbits carries a power-law time dependence, proceeds on a time scale depending on the distribution of orbital frequencies across the clump, and acts only over a portion of the phase space. Accordingly, ascertaining conditions for, and time scales of, chaotic mixing in beams is an undertaking of practical importance.
In this paper we consider a family of thermal-equilibrium (TE) configurations of beam bunches with space charge, i.e., nonneutral plasmas, confined by linear external forces [6,7,8]. For simplicity, we treat the dynamics in a reference frame that comoves with the bunch and has its origin affixed to the bunch centroid. Particle motion in this reference frame is taken to be nonrelativistic; transforming from the bunch frame to the laboratory frame is straightforward [9]. In the laboratory frame the space-charge force decreases inversely with the square of the beam energy. For the transverse component, this arises from the partial cancellation between the self-magnetic and self-electrostatic forces; while for the longitudinal component, it is due to Lorentz contraction [10]. Nonetheless, there are many situations involving high-brightness beams wherein space charge is important. Contemporary examples include low-to-medium-energy hadron accelerators such as those that drive spallation-neutron sources or serve as boosters for high-energy machines, heavy-ion accelerators, and low-energy electron accelerators such as photoinjectors [11].
Thermal-equilibrium beams are of practical interest in connection with, e.g., high-current radiofrequency linear accelerators. While conventional designs of such machines lead to bunches that are out of equilibrium, a design strategy that keeps the beam at or near thermal equilibrium has been formulated [12]. The principal motivation for this alternative strategy is to circumvent equipartitioning processes that cause emittance growth and halo formation.
Because a TE configuration is a maximum-entropy configuration, is static, and is manifestly stable [13], one might expect its intrinsic dynamics to be entirely benign. The expectation is questionable. The density distribution of such a configuration is uniform in its interior and falls to zero over a distance commensurate to the Debye length. Thus, large-amplitude orbits will explore this "Debye tail," during which time they experience a nonlinear force.
The question we seek to answer is whether the nonlinear force in the Debye tail can cause a significant number of orbits to be chaotic. The answer is unequivocally "no" for spherically symmetric or infinitely long cylindrically symmetric configurations because their potentials are integrable and thereby support only regular orbits. However, breaking the symmetry can generate chaotic orbits, as will become apparent in the analysis to follow.
Our study involves a comprehensive suite of numerical experiments concerning orbital dynamics in smooth (coarse-grained) TE configurations. We establish a quantitative measure of chaos in orbits and use this measure to distinguish between regular and chaotic orbits. We then evolve initially localized clumps of particles in the smooth potentials. The experiments are fast if the potentials are analytic, but they are much slower if the potentials must first be tabulated numerically over a grid. As part of the preliminaries, Sec. II presents a semianalytic theory for estimating the time scale for chaotic mixing. In general the TE configurations, specified in Sec. III, must be found numerically. Section IV presents a means for rapidly constructing approximate, semianalytic models of their potentials. With these models we are able to survey the parameter space and obtain a zeroth-order assessment of the prevalence and degree of chaos; this is done in Sec. V. Section VI concerns examples for which the potential is accurately determined via a numerical solution of Poisson's equation on a grid. For these examples the experiments of Sec. V are repeated, and the results are compared to those derived from the semianalytic approximation. Section VII summarizes the findings, discusses their implications while providing a comparison with the theory of Sec. II, and presents a path for follow-on work.

II. ESTIMATED TIME SCALE FOR CHAOTIC MIXING
Before embarking on numerical studies, it is wise to ascertain whether chaotic mixing can indeed proceed rapidly. One can construct an analytic tool to estimate the chaotic-mixing rate, although its application involves the tacit assumption, or initial knowledge, that chaotic orbits are present. In this section we sketch the methodology leading to analytic predictions.
The past few years have seen development of a geometric method proposed by M. Pettini to quantify chaotic instability in Hamiltonian systems with many degrees of freedom. The central idea is to describe the dynamics in terms of average curvature properties of the manifold in which the particle orbits are geodesics. The method hinges on the following assumptions and approximations; they are discussed thoroughly in Ref. [16]: (1) a generic geodesic is chaotic; (2) the manifold's effective curvature is locally deformed but otherwise constant; (3) the effective curvature reflects a gaussian stochastic process; and (4) longtime-averaged properties of the curvature are calculable as phase-space averages over an invariant measure, specifically, the microcanonical ensemble. The gaussian process is the zeroth-order term in a cumulant expansion of the actual stochastic process; assumption (3) is that the zeroth-order term suffices. The end result relates chaotic instability to the geometric properties of the manifold defined by the long-time-averaged orbits. In short, the theory is based on (often questionable) assumptions that chaos exists and is characterized by ergodicity and a microcanonical ensemble, and it treats chaotic orbits as arising from a parametric instability that can be modeled by a stochastic-oscillator equation. It has recently been adapted for application to low-dimensional, autonomous (time-independent) Hamiltonian systems and, in tests against a wide variety of such systems, it was found commonly to yield estimates of mixing rates that are good to within a factor ∼ 2 [15].
Action principles in classical mechanics are tantamount to extremals of "arc lengths"; thus, one can infer a metric tensor from an action principle [17]. The metric tensor manifests all of the properties of the manifold over which the system evolves, with these properties being calculable following standard methods of differential geometry. Of special interest is the divergence of two initially nearby 3N-dimensional geodesics q and q + δq as governed by the equation of geodesic deviation: in which D/ds denotes covariant differentiation with respect to the "proper time" s, R α βγδ is the Riemann tensor derivable from the metric tensor, and summation over repeated indices is implied with each index spanning the 3N degrees of freedom. Equation (1) is the basis for determining the mixing rate χ as a measure of the system's largest Lyapunov exponent, a quantity that reflects the long-time behavior of the separation vector: Any number of action principles, and therefore any number of metric tensors, can be selected to proceed further. Eisenhart's metric [18], which is consistent with Hamilton's least-action principle, is probably the most convenient choice. It offers easy calculation of the Riemann tensor, and it avoids spurious results traceable to the singular boundary of the perhaps better-known Jacobi metric that is derivable from Maupertius' least-action principle [19]. Eisenhart's metric operates over an enlarged configuration space-time manifold in which the geodesics are parameterized by the real time t, i.e., is the potential energy per unit mass (hereafter called the "potential"); δ ij (with the indices i, j running from 1 to 3N) is the unit tensor corresponding (without loss of generality) to a cartesian spatial coordinate system, q 0 = t; q 3N +1 = t/2 − t 0 dt ′ L(q,q); and L is the Lagrangian. The resulting geodesic equations for the spatial coordinates q i are Newton's equations of motion, so the particle trajectories correspond to a canonical projection of the Eisenhart geodesics onto the configuration space-time manifold. A convenient byproduct of the Eisenhart metric is that the only nonzero components of the Riemann tensor are Using the aforementioned assumptions and approximations, Pettini and others [16,20] derive an expression for χ in terms of the curvature and its standard deviation averaged over the microcanonical ensemble. The idea is that, as t → ∞, chaotic orbits of total energy E mix through the configuration space toward an invariant measure, taken per assumption (4) to be the microcanonical ensemble δ(H − E), over which time averages become equivalent to phase-space averages. Specifically, for an arbitrary function A(q), the averaging process Per Eisenhart's metric, the average curvature κ and the ratio ρ ≡ σ/κ, with σ denoting the standard deviation of the curvature, are in which ∇ 2 denotes the Laplacian ∂ i ∂ i , and ρ corresponds physically to the ratio of the average curvature radius to the length scale of fluctuations [21]. By taking the curvature to vary randomly along a chaotic orbit, one can reduce Eq. (1) to a stochastic-oscillator equation that can be solved analytically. The solution yields an estimate of the largest Lyapunov exponent χ: The geometric quantities derive from the 6N-dimensional microcanonical distribution.
Anticipating that granularity takes a long time to affect mixing, and wishing to identify conditions for rapid mixing, we now consider the influence of the 3-dimensional coarse-grained space-charge potential V s on a generic chaotic orbit. The largest Lyapunov exponent for the coarse-grained system equates to the chaotic-mixing rate. We presume the assumptions and approximations stated at the outset carry over to the coarse-grained system; the main justification is that the aforementioned previous work concerning low-dimensional autonomous Hamiltonians has shown the mixing rate in such systems usually depends only weakly on the dynamical details [15]. We take the external focusing potential V f to be quadratic in the coordinates x comoving with the bunch, i.e., V f (x) = (ω · x) 2 /2, wherein ω = (ω x , ω y , ω z ) corresponds to the focusing strength; the total potential is V = V f + V s . Per Eq. (4) and Poisson's equation the quantities κ and σ are determined from is the (smoothed) particle density, q and m are the single-particle charge and rest mass, respectively, and ǫ o is the permittivity of free space. We then have Inserting these results into Eq. (5) gives the associated time scale for chaotic mixing, t m ≡ 1/χ. When the standard deviation of the density distribution is large, as can be the case when substructure is present, ρ will be appreciable, and in turn Eq. (5) makes clear that t m will be a few space-charge-depressed periods 2π/ √ κ. Accordingly, the spacecharge-depressed period, a quantity commensurate to the orbital period of a typical particle, constitutes a "dynamical time" t D for charged-particle beams.
To underscore the potential impact of collisionless relaxation via chaotic mixing, it is of interest to compare t m to the collisional relaxation time t R . Perhaps the simplest way to develop an order-of-magnitude estimate of t R in a charged-particle bunch (a nonneutral plasma) is to calculate the time required for a typical particle velocity to change by of order itself presuming collisions comprise a sum of incoherent binary interactions [22]. The result is t R /t D ∼ 0.1N/lnN, wherein the Coulomb logarithm is conservatively taken to be lnN. If we substitute plausible parameter values for real high-brightness beams, we find t R ≫ t D ; for example, N = 6.25 × 10 9 (1 nC) gives t R ∼ 10 7 t D ; hence, t m ≪ t R when chaotic mixing is prominent. The remaining question is whether there is a significant population of globally chaotic orbits to mix, a question to which we now turn our attention.

III. THE EQUATIONS OF THERMAL EQUILIBRIUM
Consider a system, i.e., a bunch, of N identical charged particles, e.g., electrons or protons. For simplicity, invoke a Cartesian coordinate system whose origin lies at the bunch centroid. Assume all particle velocities in this coordinate system are nonrelativistic. The particles mutually interact via the Coulomb force and are confined by a static, externally applied, linear focusing force. The focusing force may have different strengths along the three Cartesian axes. Assume, apart from this focusing force, that the system is isolated and is in thermal equilibrium. Accordingly, the total energy E of each particle is conserved: wherein ω = (ω x , ω y , ω z ) corresponds to the focusing strength; x = (x, y, z) denotes coordinates; m, v, and q are the particle's rest mass, speed, and charge, respectively; and φ(x) = (m/q)V s is the space-charge potential arising from the collective Coulomb force.
To proceed, one would in principle work with the 6N-dimensional microcanonical distribution of particles. This distribution includes interactions at all scales, ranging from particle-on-particle to a single particle interacting with the bulk, smooth potential from all other particles. Discreteness effects from 1/r 2 particle collisions generate chaos [23]; they cause nearby particle trajectories to separate exponentially. The rate of exponential separation, i.e., the Lyapunov exponent, is an increasing function of N [24]. In this sense, larger N gives rise to more chaos. However, the scale at which the separation saturates is a decreasing function of N. Accordingly, in large-N, high-charge-density systems such as beams with space charge, discreteness establishes microchaos [25,26,27,28]. At the other extreme, that of a single particle interacting with the bulk, smooth potential, exponential separation of nearby chaotic particles (if any are present) saturates at a global scale, corresponding to a state of macrochaos. Thus, initially nearby chaotic orbits evolve in three stages [29]: (1) very rapid exponential divergence that saturates at a scale large compared to the initial interparticle spacing but small compared to the system size; followed by (2) rapid exponential divergence that persists until the particles are globally dispersed; followed by (3) less rapid power-law divergence on a time scale ∝ (ln N)t D , in which t D is a dynamical time commensurate to the orbital period. If, in the smooth potential, the initially nearby particles execute regular motion rather than chaotic, then stage (2) is absent, and stage (3) proceeds on the much longer time scale ∝ (N 1/2 )t D [27].
Our interest here is in stage (2). Specifically, we are concerned about the existence of, and time scale for, macroscopic chaos, i.e., chaotic mixing into the global region of phase space that is energetically accessible to the individual particles. Accordingly, we specialize to the smooth 6-dimensional distribution function of a single particle, recognizing that discreteness effects vanish on macroscopic scales as the number density grows. For the TE beam, this the Hamiltonian, k is Boltzmann's constant, and T is the beam temperature. The number density follows upon integrating over velocity space, and the space-charge potential follows upon solving Poisson's equation: wherein ǫ o is the permittivity of free space.
A much more convenient formulation arises by using dimensionless variables. We introduce the Debye length λ D0 and angular plasma frequency ω p0 , both defined in terms of the centroid number density n(0): We then measure all lengths in the unit of λ D0 , i.e., x ↔ x/λ D0 , and all times in the unit of 1/ω p0 , i.e., t ↔ ω p0 t. In addition, we introduce the dimensionless potential Φ(x) ≡ qφ(x)/(kT ), and we normalize n(x) to the centroid density n(0), i.e., n(x) ↔ n(x)/n(0).
Equations (11) and (12)  This is the minimum permissible focusing strength; the bunch is unconfined if Ω < Ω u , and the corresponding constraint on the parameter space is Hence, the parameter set [a, c; Ω] fully specifies a TE configuration.
Upon solving for the space-charge potential Φ(x), one can calculate orbits of test particles in the total potential. Their trajectories follow from the (dimensionless) equation of motion: One can, of course, introduce arbitrary initial conditions for the orbits. In our experiments, the initial condition on the velocity is v(0) = 0, and the total energy E of a particle thereby corresponds to the potential energy associated with the initial position x(0).
A key challenge in exploring orbital dynamics throughout the parameter space is to integrate large numbers of orbits rapidly for sufficiently long evolutionary times. Ideally, one would have analytic solutions for the density-potential pairs, from which the force on a particle at each time step can be quickly evaluated. Unfortunately, the equations of equilibrium generally do not submit to analytic techniques. Thus, in principle, one must solve these equations numerically, e.g., over a grid. However, as delineated in the following section, it is possible to formulate approximate, semianalytic solutions, and these solutions enable a search of a broad range of the parameter space for regions that support chaotic orbits. We now turn to that exploration. Subsequently, for select cases, we compare these results against those derived from fully self-consistent numerical solutions.

IV. APPROXIMATE SOLUTIONS TO THE EQUATIONS OF EQUILIBRIUM
A method to solve the equations of equilibrium is through a sequence of successive approximations [30]. A way to begin such a sequence is as follows: (1) As a first approximation, represent the system as a configuration stratified on similar and similarly situated concentric ellipsoids. A "homeoid" is a shell that is bounded by two similar and similarly situated concentric ellipsoids, and in which the surfaces of constant density are ellipsoids that are similar to, concentric with, and similarly situated with respect to the bounding ellipsoid.
Thus, the charge density is "homeoidally striated" in the first approximation (as is later illustrated in Fig. 13). Determine the stratification by solving a spherically symmetric model of the equations of equilibrium. (2) In the second approximation, derive the space-charge field corresponding to the homeoidally striated charge density, and then solve exactly the equations of equilibrium in this field. (3 and up) Repeat the process until the density and potential converge. In practice, one can carry out steps (1) and (2) of this recipe using semianalytic methods; to go further requires numerical techniques.
A. Determination of the structure in the first approximation To invoke a spherically symmetric model of Eq. (12), we take the potential to be stratified over ellipsoids on which R(x) takes a constant value. Then the spherically symmetric model corresponds to solving This model defines the "zeroth approximation" Φ 0 (R) to the potential. In general Eq. (15) must be solved numerically; however, the solution is rapidly and easily accomplished with the aid of, e.g., a Runge-Kutta algorithm.
Once Φ 0 (R) is determined, the corresponding homeoidally striated density becomes the first approximation to the number density: By inspection [31], one can write down the space-charge potential corresponding to the number density n 1 (R), and this becomes the first approximation to the potential: wherein the second equality follows from an integration by parts, and the quantities ∆(u) and R(x; u) are ∆(u) = (a 2 + u)(1 + u)(c 2 + u) , Hence, in the first approximation the number density is homeoidally striated, but the spacecharge potential is not.
B. Determination of the structure in the second and higher approximations The number density in the second approximation, n 2 (x), follows upon substituting Φ 1 (x) calculated from Eq. (17) into Eq. (11). For the special case of spherical symmetry, all orders of approximation agree with one another, but this is of course not true for a general triaxial geometry. To go further requires numerical methods, e.g., solving Poisson's equation for the potential Φ 2 corresponding to the density n 2 , substituting the result into Eq. (11) to obtain n 3 , and successively repeating the process until convergence is achieved. As discussed in Sec. VI below, we use a different method, a multigrid algorithm, for solving Eqs. (11) and (12) numerically.

V. SURVEY OF THE PARAMETER SPACE
Gathering sufficient data to support precise, statistically based conclusions concerning orbital behavior in a given potential requires integrating thousands of orbits in that potential.
And before these orbits can be tracked, the potential needs to be ascertained to sufficient accuracy. In principle, and for each choice of parameters, one must construct the "exact" potential Φ(x) by numerically solving the corresponding Poisson equation. This can be a computationally tedious process, and the solution is by necessity defined over a grid. Next, orbit integration through the grid requires accurate interpolation to evaluate the potential and corresponding particle acceleration between grid points. For sufficient resolution, the time steps need to be appropriately small; accordingly, many interpolations are required, and integrating many orbits is computationally time-consuming. This process is feasible for studying a few choices of parameter sets, and it underlies the results of Sec. VI below.
However, to survey the entire parameter space, i.e., to investigate many choices of parameter sets, the process becomes prohibitive. For this purpose one must resort to using approximate potentials.
Sec. IV above details a sequence of approximations, the first elements of which are semianalytic. The zeroth-order potential Φ 0 , derived from Eq. (15), is easy to evaluate, and it enables fast, high-precision orbital integration. However, Φ 0 itself may be a crude approximation to the exact potential; the approximation gets progressively worse as the parameter sets deviate further from spherical symmetry. One might expect the potential Φ 1 of the first approximation to provide a better model. However, its underlying integral, given in Eq. (17), adds additional complexity and time to the orbit integrations. We tried evaluating this integral at each time (thus position) step along the orbit, but doing so made the orbit computations prohibitively long. The alternative is to evaluate the integral over a grid and then do orbit integrations through the grid. As previously mentioned, integrations through a grid are too computationally expensive to enable a parameter survey. Moreover, if one is able to solve Poisson's equation for the exact potential Φ(x), then there is neither computational benefit nor motivation for using Φ 1 . Our strategy is to explore a few choices of parameter sets in the exact potential to strengthen conclusions from our survey, and this necessitated developing the Poisson solver described in Sec. VI. For these reasons, we use the potential of the zeroth approximation, Φ 0 , to survey the parameter space. For a few specific parameter sets for which the results of the zeroth approximation look especially interesting, we then check the results using the numerically evaluated exact potential, Φ(x).
Plots of Φ 0 (R) versus R derived from Eq. (15) appear in Fig. 1. Also shown are the corresponding profiles of the number densities n 1 (R) constructed in the first approximation.
For larger "case numbers" i, the density contains larger quasi-uniform central regions. In the outer regions the density decreases, over a length commensurate to the Debye length, to a low-density tail. The space-charge force in the quasi-uniform "core" is correspondingly quasilinear; however, it is manifestly nonlinear in the "Debye fall-off region" (henceforth called the "Debye tail"). Fig. 1 shows that the choices of Ω i per Eq. (18)  The integrations were done using a fifth-order Runge-Kutta algorithm [32] with variable time step. As the integration proceeded, we computed the largest short-time Lyapunov exponent of each orbit using a well-established algorithm in the field of chaotic dynamics [33]. The idea is to evolve two initial conditions that start from a very close distance for about one dynamical time, then renormalize to bring the two particles close together again, and repeat the process until the average exponent associated with the orbital separation converges to an almost stable value. Typically convergence was achieved within ∼ 100 orbital periods.
After computing the orbits, we extracted the power spectrum for each orbit using a fast-Fourier-transform algorithm [32]. In doing so, we recorded each orbit at a rate ∼ 40 times per orbital period. From the spectrum we computed the total power. Then we sorted the spectral frequencies in descending order, and starting from the highest frequency we added as many frequencies as were needed to reach 90% of the total power. The required number of frequencies is defined to be the "complexity" n of the orbit [34].

Criterion for chaos
Our first and foremost interest is to determine how many of the 2000 orbits in our sample are chaotic in a given TE configuration. Accordingly an objective, quantitative criterion for chaos is needed. There is no universally accepted criterion; hence, we developed our own using the following rationale. Both the largest short-time Lyapunov exponent χ and the complexity n are well-established, conventional measures of chaos [35]. A first piece of information for defining the criterion comes from plotting n versus χ for all of the orbits. We chose to base our criterion only on the complexity for the following reason. Though no problem arises in computing Lyapunov exponents in the zeroth approximation, such is not the case in the exact potential, wherein we found that longer integration times are required to achieve adequate convergence. Recall that the exact potential must be specified over a three-dimensional grid. Accordingly, interpolation errors and discontinuities between the cells can affect computations of Lyapunov exponents because they involve the distance between two initially nearby orbits, which is a local property that is sensitive to the grid size and the order of interpolation. By contrast, a computation of complexity, in that it involves the Fourier spectrum of an individual orbit, avoids reference to nearby orbits and is thus a global property influenced little by the grid size, a notion that we corroborated during the course of our numerical studies. For simulations in exact potentials, we chose a grid size and interpolation algorithm (cf. Sec. VI below) such that numerical errors had negligible effect on the computation of individual orbits. A standard measure of the "goodness" of an orbital integration is the degree to which total energy is conserved [36]; for every orbit we achieved conservation of total energy with relative error ≤ 10 −6 . This is some two orders of magnitude better than contemporary standard practice.
Investigations of n-vs.-χ plots for many zeroth-order TE potentials led us to a specific quantitative criterion, namely, an orbit is categorized as chaotic if its complexity n > 20 [cf. here to be the orbital period corresponding to the total energy of the individual particles comprising the clump. After some tens of t D each clump has spread through a volume commensurate to the total particle energy.

C. Survey results
Our strategy for surveying the parameter space of the TE configurations is as follows. We conduct the survey using the zeroth-order potential Φ 0 found from Eq. (17). Recall that this potential depends on the focusing strength Ω and, through R(x), the scale lengths (a, c);

VI. NUMERICAL EXPERIMENTS IN EXACT POTENTIALS
As a matter of principle, one must be concerned about the extent to which subtle structure in the potential can influence the qualitative behavior, and in particular the chaoticity, of an orbit. One well-known example is that of the Toda potential; the full Toda potential is integrable and supports only regular orbits, but generally a truncated Toda potential is not integrable and supports a population of chaotic orbits [3]. Our survey of the parameter space of TE configurations centered on the use of Φ 0 , a generally crude approximation to the true potential. The survey suggests a large region of the parameter space supports sizeable populations of chaotic orbits, wherein all of these orbits reach into the Debye tail.
We may expect in general that the density profile, particularly that of the Debye tail, corresponding to the exact potential is considerably different from that corresponding to Φ 0 . For example, in the limit of distances very far from the centroid, the exact space-charge potential will approach spherical symmetry, whereas Φ 0 is everywhere homeoidally striated. Accordingly, to check and have confidence in the qualitative results of Sec. V, we must repeat the numerical experiments in a suitably broad collection of exact potentials. As mentioned earlier, the reason we did not base the survey on exact potentials is that the respective numerical experiments are computationally expensive.
To integrate Eq. (12) governing the exact potential Φ(x), which is a fully threedimensional partial differential equation (PDE), we chose a multigrid algorithm [37]. The algorithm requires boundary conditions be specified over the surface of the volume occupied by the grid. We chose a cubic grid volume greatly exceeding the volume of interest, i.e., that spanning the Debye fall-off of the density. Then we calculated the boundary conditions over the surface of this volume using the formalism of the first approximation, specifically, Eq. (17). Because the resulting boundary conditions are only first approximations to the true boundary conditions, we checked our numerical solutions by varying the positions of the bounding surfaces of the grid by factors of five, and we found negligible change in the results over the volume of interest. Applying the multigrid algorithm to three dimensions involves nontrivial manipulations of the inherent restriction, interpolation, and relaxation routines. In the process, a nonlinear algebraic equation emerges due to the nonlinearity of the PDE. It was solved using an iterative method that combines Newton-Raphson and bisection techniques [32].
After tabulating the exact potential in three dimensions, we used a fifth-order Runge-Kutta algorithm with variable time step to evolve the individual orbits in time, within which the force at every time (thus position) step was computed using a three-dimensional interpolation scheme. The resulting orbits conserved total energy with relative error better than 10 −6 and sometimes as low as 10 −7 . Figure 9 exemplifies the difference between the zeroth-order and exact potentials. The top two panels present isopotential contours in the (x, z)-plane for Case 5 with a 2 = 0.5, c 2 = 1.5, a triaxial configuration. The bottom panels show how the two potentials compare along each of the (x, y, z)-axes. Obviously there are, and there should be, differences, but the important question is to what extent these differences alter the qualitative evolution of the orbits and, in turn, the complexities that characterize them? Analogous graphs for Case 5 and a 2 = 1.0, c 2 = 0.25, a strongly oblate spheroid, appear in Fig. 10, from which the respective differences are seen to be much more pronounced. Figure 11 provides a visual comparison of a chaotic orbit starting from the same initial condition and evolving in the zeroth-order and exact potentials of Fig. 9. Although orbits in the two potentials differ quantitatively, in many cases they are qualitatively similar in that they explore a similar volume of phase space and have similar morphology. This pertains to the example of Fig. 11; however, this one example does not in any way guarantee every orbit that is chaotic in the potential of the zeroth approximation is also chaotic in the exact potential. Statistical comparisons of orbital complexities respective to the zeroth-order and exact potentials for several Case 5 configurations appear in Fig. 12, and these show that the complexities in the two potentials can differ considerably depending on the specific parameter set under study.
Following the procedure delineated in Sec. V C, we also computed the percentage of chaotic orbits in a broad range of exact Case 5 potentials. The results, juxtaposed against their counterparts computed using Φ 0 , appear in Table I and in Figs. 5 and 6. For these examples, there is generally a smaller percentage of chaotic orbits in the exact potential than in the zeroth approximation. The explanation is simple: compared to the density n 1 [R(x)], the density derived from the exact space-charge potential Φ(x) is quasi-uniform over a larger volume and falls to small values over a shorter scale length. Accordingly, the configurationspace volume over which the space-charge force is markedly nonlinear, i.e., the Debye tail, is smaller. Figure 13 illustrates the difference in the density profiles corresponding to the approximate and exact solutions. In the limit of spherical symmetry the profiles are identical, and they disagree more strongly as they become less spherically symmetric. Most notable is the comparison between Figs. 13(e) and (f) concerning a strongly oblate spheroid, where we see that the corresponding exact density distribution is much more uniform than that of the zeroth approximation. This accounts for the strong discrepancy revealed in Fig. 12 (d) concerning orbital chaoticity. It is also consistent with expectations based on first principles: the closer a system is to being one-dimensional (e.g., sphere, cylindrically symmetric disc, infinite symmetric cylinder, in which particle motion is integrable), the less is the population of chaotic orbits. However, the essential observation is that, in most cases, the exact TE configurations do indeed support substantial populations of chaotic orbits in keeping with expectations that surfaced from the survey based on approximate solutions.

VII. SUMMARY, IMPLICATIONS, AND FUTURE WORK
We have explored orbital dynamics and phase mixing in thermal-equilibrium beams for which the potential is the superposition of an external potential quadratic in the coordinates and the self potential arising from space charge. The associated parameter space spans the full range of symmetries, i.e., spherical, cylindrical, and triaxial, and the full range of density profiles, i.e., gaussian (corresponding to negligible space charge) through uniform (corresponding to maximal space charge). To reiterate, the main findings concerning chaos in these systems, "discovered" in the context of zeroth approximations to the space-charge potentials and affirmed with the respective exact potentials, are: (1) configurations corresponding to a large portion of the parameter space support considerable populations of chaotic orbits, (2) essentially all of the orbits that are chaotic reach into the Debye tail where the collective space-charge force is manifestly nonlinear, (3) prolate axisymmetric configurations support little chaos, but prolate triaxial configurations can support considerable chaos, and (4) strongly oblate spheroids support little chaos, but moderately oblate spheroids can support considerable chaos.
It is of interest to compare theoretical predictions concerning TE configurations, for which we herein have established the existence of chaotic orbits, with results of our numerical experiments. In terms of the dimensionless quantities introduced in Sec. III, the parameters κ and ρ of Eq. (6) take the form and Eq. (5) then yields the mixing rate χ. A comparison between theory and numerical experiments appears in Fig. 14 wherein the simulation results reflect statistics from initially localized clumps of 2000 particles that were started at zero velocity at various points in configuration space corresponding to various total particle energies E. The figure presents a plot of the mixing rate χ versus |E| in the Case 5 configuration with a 2 = 4/5, c 2 = 4/3, a slightly triaxial system. This configuration is "not too far away" from spherical symmetry, which means the zeroth approximation Φ(x) = Φ 0 is correspondingly reasonable. It also means only a modest population (∼ 5% for this parameter set) of chaotic orbits is supported.
The figure was derived within the framework of the zeroth approximation because therein the Lyapunov exponents, i.e., mixing rates, can be accurately computed from the simulations and the microcanonical averages required for the theory likewise can be easily and accurately evaluated. The numerical experiments span a range 0.5 ≤ |E| ≤ 60, corresponding to 11 ≤ R ≤ 25, i.e., extending from within to well beyond the Debye drop-off in the density profile. The agreement between theory and numerical experiments is remarkably close.
One can see from the numerical experiments described herein that chaotic mixing takes place on an e-folding time scale comparable to a dynamical time (an orbital period). This is very fast compared to, e.g., collisional relaxation; hence, one must account for this collisionless process when designing an accelerator for the production of high-peak-current, high-brightness beams. For example, particles comprising a beam out of equilibrium will, if globally chaotic, redistribute themselves globally and irreversibly on a dynamical time scale.
Because perturbations induced, e.g., by transitions in the beamline will drive a beam away from equilibrium, chaotic mixing can be a dynamic of practical importance. Consider the case of a TE configuration: a small perturbation from image charges passing through an external irregularity in the beamline will distort the Debye tail. If a substantial fraction of particles in the Debye tail are chaotic, which is the case for a wide range of bunch geometries, a corresponding fraction of the orbits comprising the distortion will quickly mix throughout the volume of the configuration. The work done by the external perturbation in setting up the distortion will thereby appear in the form of a larger configuration-space volume. If the perturbation is strong enough so that mixing in momentum space associated with consequent time-dependence in the potential is also substantial, then some of the work done will also appear in the form of a larger momentum space. The net effect is a larger emittance.
If there are many such perturbations along the beamline, the cumulative emittance growth may be troublesome.
The present investigation and its associated implications concern only very specific, timeindependent, single-species systems, i.e., beams (or nonneutral plasmas) in thermal equilibrium. These are the most benign systems imaginable, yet we found even they can support chaotic orbits. Any perturbation will create a nonequilibrium, time-dependent system that will subsequently evolve self-consistently. Accordingly, the space-charge potential can be complicated, particularly if the perturbation is strong. The only sensible conjecture under such conditions is that the corresponding population of chaotic orbits will be larger, and in turn chaotic mixing will be more prevalent. Exploratory numerical simulations of an equipartitioning system and of merging beamlets have supported this notion [2]. Further exploration of time-dependent beams is warranted and will likely prove illuminating, particularly in regard to deciphering time scales for emittance growth, halo formation, etc.
By using only smooth potentials we have restricted our analysis to the six-dimensional phase space of a single particle. Accordingly we have suppressed dissipative effects of collisions in particular, and force fluctuations in general. Such effects can only enhance chaos, as has been demonstrated, e.g., in numerical experiments concerning self-gravitating systems [27]. As the next step, we have constructed frozen N-body representations of the charge densities of the TE configurations and with these representations are repeating the numerical experiments described herein. One of our objectives is to determine the minimum number of particles needed to reproduce the dynamics associated with smooth timeindependent potentials. Results will be described in a forthcoming paper [28]. In the future it will be of interest to do likewise for time-dependent systems and ultimately ascertain, e.g., conditions under which the Vlasov equation governing the six-dimensional phase space of a single particle can be applied with confidence.