text:
# Chaotic orbits in thermal-equilibrium beams: existence and dynamical implications

## 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 3$`N`$ 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-charge-depressed betatron periods reiser94 . Simulations of the experiment reveal a substantial fraction of globally chaotic orbits bohn02 , 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” kandrup94 ; merritt96 . 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 sagan94 , 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 davidson99 ; brown95 ; reiser93 . 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 reiser2 . 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 chao93 . 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 lee96 .

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 reiser95 . 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 davidson98 , 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. Additional details are available elsewhere bohn00 ; kandrup01 .

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. casetti96 : (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) long-time-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`$ kandrup01 .

Action principles in classical mechanics are tantamount to extremals of “arc lengths”; thus, one can infer a metric tensor from an action principle goldstein50 . 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 $`𝐪`$ and $`𝐪+\delta 𝐪`$ as governed by the equation of geodesic deviation:

$$\frac{D^2\delta q^\alpha }{ds^2}+R_{}^{\alpha }{}_{\beta \gamma \delta }{}^{}\frac{dq^\beta }{ds}\delta q^\gamma \frac{dq^\delta }{ds}=0,$$
(1)
in which $`D/ds`$ denotes covariant differentiation with respect to the “proper time” $`s`$, $`R_{}^{\alpha }{}_{\beta \gamma \delta }{}^{}`$ 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 $`\chi `$ as a measure of the system’s largest Lyapunov exponent, a quantity that reflects the long-time behavior of the separation vector:

$$\chi =\underset{t\mathrm{}}{lim}\frac{1}{t}\mathrm{ln}\frac{|\delta 𝐪(t)|}{|\delta 𝐪(0)|}.$$
(2)
Any number of action principles, and therefore any number of metric tensors, can be selected to proceed further. Eisenhart’s metric eisenhart29 , 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 szczesny99 . Eisenhart’s metric operates over an enlarged configuration space-time manifold in which the geodesics are parameterized by the real time $`t`$, i.e., $`ds^2=dt^2=2V(𝐪)(dq^0)^2+\delta _{ij}dq^idq^j+2dq^0dq^{3N+1}`$, in which $`V(𝐪)`$ is the potential energy per unit mass (hereafter called the “potential”); $`\delta _{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_0^t𝑑t^{}L(𝐪,\dot{𝐪})`$; 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 $`R_{0i0j}=_i_jV`$, in which $`_i=/q^i`$.

Using the aforementioned assumptions and approximations, Pettini and others casetti96 ; cipriani98  derive an expression for $`\chi `$ in terms of the curvature and its standard deviation averaged over the microcanonical ensemble. The idea is that, as $`t\mathrm{}`$, chaotic orbits of total energy $`E`$ mix through the configuration space toward an invariant measure, taken per assumption (4) to be the microcanonical ensemble $`\delta (HE)`$, over which time averages become equivalent to phase-space averages. Specifically, for an arbitrary function $`A(𝐪)`$, the averaging process is

$$A\underset{t\mathrm{}}{lim}A_t=\frac{𝑑𝐪𝑑\dot{𝐪}A(𝐪)\delta [H(𝐪,\dot{𝐪})E]}{𝑑𝐪𝑑\dot{𝐪}\delta [H(𝐪,\dot{𝐪})E]}.$$
(3)
Per Eisenhart’s metric, the average curvature $`\kappa `$ and the ratio $`\rho \sigma /\kappa `$, with $`\sigma `$ denoting the standard deviation of the curvature, are

$$\kappa =\frac{^2V}{3N1},\rho =\frac{1}{\kappa }\frac{\sqrt{(^2V)^2^2V^2}}{\sqrt{3N1}},$$
(4)
in which $`^2`$ denotes the Laplacian $`_i^i`$, and $`\rho `$ corresponds physically to the ratio of the average curvature radius to the length scale of fluctuations casetti95 . 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 $`\chi `$:

$`\chi (\rho )`$ $`=`$ $`{\displaystyle \frac{1}{\sqrt{3}}}{\displaystyle \frac{L^2(\rho )1}{L(\rho )}}\sqrt{\kappa };`$
$`L(\rho )`$ $`=`$ $`\left[T(\rho )+\sqrt{1+T^2(\rho )}\right]^{1/3},T(\rho )={\displaystyle \frac{3\pi \sqrt{3}}{8}}{\displaystyle \frac{\rho ^2}{2\sqrt{1+\rho }+\pi \rho }}.`$ (5)
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 kandrup01 . We take the external focusing potential $`V_f`$ to be quadratic in the coordinates $`𝐱`$ comoving with the bunch, i.e., $`V_f(𝐱)=(\omega 𝐱)^2/2`$, wherein $`\omega =(\omega _x,\omega _y,\omega _z)`$ corresponds to the focusing strength; the total potential is $`V=V_f+V_s`$. Per Eq. (4) and Poisson’s equation the quantities $`\kappa `$ and $`\sigma `$ are determined from $`^2V=\omega _f^2\omega _p^2(𝐱)`$, in which $`\omega _f^2=\omega _x^2+\omega _y^2+\omega _z^2`$, $`\omega _p^2(𝐱)=n(𝐱)q^2/(ϵ_om)`$, $`n(𝐱)`$ 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

$$\kappa =\frac{1}{2}\left[\omega _f^2\omega _p^2(𝐱)\right],\rho =\frac{1}{\kappa }\frac{\sqrt{n^2(𝐱)n(𝐱)^2}}{n(0)\sqrt{2}}.$$
(6)
Inserting these results into Eq. (5) gives the associated time scale for chaotic mixing, $`t_m1/\chi `$. When the standard deviation of the density distribution is large, as can be the case when substructure is present, $`\rho `$ will be appreciable, and in turn Eq. (5) makes clear that $`t_m`$ will be a few space-charge-depressed periods $`2\pi /\sqrt{\kappa }`$. Accordingly, the space-charge-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 rosenbluth . The result is $`t_R/t_D0.1N/\text{ln}N`$, wherein the Coulomb logarithm is conservatively taken to be $`\text{ln}N`$. If we substitute plausible parameter values for real high-brightness beams, we find $`t_Rt_D`$; for example, $`N=6.25\times 10^9`$ (1 nC) gives $`t_R10^7t_D`$; hence, $`t_mt_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:

$$E=\frac{1}{2}mv^2+\frac{1}{2}m(\omega 𝐱)^2+q\varphi (𝐱);$$
(7)
wherein $`\omega =(\omega _𝐱,\omega _𝐲,\omega _𝐳)`$ corresponds to the focusing strength; $`𝐱=(x,y,z)`$ denotes coordinates; $`m`$, $`v`$, and $`q`$ are the particle’s rest mass, speed, and charge, respectively; and $`\varphi (𝐱)=(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 struckmeier96 ; they cause nearby particle trajectories to separate exponentially. The rate of exponential separation, i.e., the Lyapunov exponent, is an increasing function of $`N`$ hemsendorf02 . 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 $`\mathrm{𝑚𝑖𝑐𝑟𝑜}`$chaos kansid01 ; sideris02 ; kansid02 ; Nbody . 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 $`\mathrm{𝑚𝑎𝑐𝑟𝑜}`$chaos. Thus, initially nearby chaotic orbits evolve in three stages kandrup02 : (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 $`(\mathrm{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`$ kansid02 .

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 is just the Maxwell-Boltzmann distribution, $`f(𝐱,𝐯)\mathrm{exp}(H/kT)`$, in which $`H=E`$ is 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:

$$n(𝐱)=n(0)\mathrm{exp}\left[\frac{\frac{1}{2}m(\omega 𝐱)^2q\varphi (𝐱)}{kT}\right];$$
(8)
$$^2\varphi (𝐱)=\frac{q}{ϵ_o}n(𝐱),\varphi (𝐱=0)=\varphi (𝐱=0)=0,$$
(9)
wherein $`ϵ_o`$ is the permittivity of free space.

A much more convenient formulation arises by using dimensionless variables. We introduce the Debye length $`\lambda _{D0}`$ and angular plasma frequency $`\omega _{p0}`$, both defined in terms of the centroid number density $`n(0)`$:

$$\lambda _{D0}^2\frac{ϵ_okT}{n(0)q^2};\omega _{p0}^2\frac{n(0)q^2}{ϵ_om}.$$
(10)
We then measure all lengths in the unit of $`\lambda _{D0}`$, i.e., $`𝐱𝐱/\lambda _{D0}`$, and all times in the unit of $`1/\omega _{p0}`$, i.e., $`t\omega _{p0}t`$. In addition, we introduce the dimensionless potential $`\mathrm{\Phi }(𝐱)q\varphi (𝐱)/(kT)`$, and we normalize $`n(𝐱)`$ to the centroid density $`n(0)`$, i.e., $`n(𝐱)n(𝐱)/n(0)`$. The number density and Poisson’s equation then reduce to their dimensionless forms:

$$n(𝐱)=\mathrm{exp}\left[\frac{1}{2}\mathrm{\Omega }^2R^2(𝐱)\mathrm{\Phi }(𝐱)\right];$$
(11)
$$^2\mathrm{\Phi }(𝐱)=n(𝐱),\mathrm{\Phi }(𝐱=0)=\mathrm{\Phi }(𝐱=0)=0;$$
(12)
wherein $`\mathrm{\Omega }^2(\omega _y/\omega _{p0})^2\omega _y^2`$, and $`R^2(𝐱)=(x/a)^2+y^2+(z/c)^2`$, with the “scale lengths” $`a`$ and $`c`$ defined as $`a\omega _y/\omega _x`$ and $`c\omega _y/\omega _z`$.

Equations (11) and (12) self-consistently provide the structure of the entire family of smoothed TE configurations. The parameter $`\mathrm{\Omega }`$ governs the strength of external focusing vis-$`\stackrel{`}{\text{a}}`$-vis space charge. The scale lengths $`a`$ and $`c`$ set the overall geometry: $`a=c=1`$ corresponds to spherical symmetry, $`a=c1`$ corresponds to cylindrical symmetry, and $`a1`$, $`c1`$, $`ac`$ establishes a triaxial configuration. The extreme case of maximum space charge corresponds to a density that is strictly uniform over the volume of the configuration, in which case $`\mathrm{\Phi }=\mathrm{\Omega }^2R^2/2`$ inside the beam. Upon substituting into Poisson’s equation, we see that the associated configuration carries the parameter $`\mathrm{\Omega }=\mathrm{\Omega }_u=1/\sqrt{(1/a^2)+1+(1/c^2)}`$. This is the minimum permissible focusing strength; the bunch is unconfined if $`\mathrm{\Omega }<\mathrm{\Omega }_u`$, and the corresponding constraint on the parameter space is

$$\frac{1}{a^2}+\frac{1}{c^2}\frac{1\mathrm{\Omega }^2}{\mathrm{\Omega }^2}.$$
(13)
Hence, the parameter set $`[a,c;\mathrm{\Omega }]`$ fully specifies a TE configuration.

Upon solving for the space-charge potential $`\mathrm{\Phi }(𝐱)`$, one can calculate orbits of test particles in the total potential. Their trajectories follow from the (dimensionless) equation of motion:

$$\frac{d^2𝐱}{dt^2}=\left[\frac{1}{2}\mathrm{\Omega }^2R^2(𝐱)+\mathrm{\Phi }(𝐱)\right].$$
(14)
One can, of course, introduce arbitrary initial conditions for the orbits. In our experiments, the initial condition on the velocity is $`𝐯(0)=0`$, and the total energy $`E`$ of a particle thereby corresponds to the potential energy associated with the initial position $`𝐱(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 bohn83 . 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.

### IV.1 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(𝐱)`$ takes a constant value. Then the spherically symmetric model corresponds to solving

$`{\displaystyle \frac{1}{R^2}}{\displaystyle \frac{d}{dR}}\left(R^2{\displaystyle \frac{d\mathrm{\Phi }_0(R)}{dR}}\right)`$ $`=`$ $`\mathrm{exp}\left[{\displaystyle \frac{1}{2}}\mathrm{\Omega }^2R^2\mathrm{\Phi }_0(R)\right];`$
$`\mathrm{\Phi }_0(0)`$ $`=`$ $`{\displaystyle \frac{d\mathrm{\Phi }_0}{dR}}|_{R=0}=0.`$ (15)
This model defines the “zeroth approximation” $`\mathrm{\Phi }_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 $`\mathrm{\Phi }_0(R)`$ is determined, the corresponding homeoidally striated density becomes the first approximation to the number density:

$$n_1(R)=\mathrm{exp}\left[\frac{1}{2}\mathrm{\Omega }^2R^2\mathrm{\Phi }_0(R)\right].$$
(16)
By inspection efe , 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:

$`\mathrm{\Phi }_1(𝐱)`$ $`=`$ $`{\displaystyle \frac{ac}{2}}{\displaystyle _0^{\mathrm{}}}{\displaystyle \frac{du}{\mathrm{\Delta }(u)}}{\displaystyle _0^{R(𝐱;u)}}𝑑rrn_1(r)`$ (17)
$`=`$ $`{\displaystyle \frac{ac}{2}}{\displaystyle _0^{\mathrm{}}}{\displaystyle \frac{du}{\mathrm{\Delta }(u)}}\left[r{\displaystyle \frac{d\mathrm{\Phi }_0(r)}{dr}}+\mathrm{\Phi }_0(r)\right]_{r=R(𝐱;u)},`$
wherein the second equality follows from an integration by parts, and the quantities $`\mathrm{\Delta }(u)`$ and $`R(𝐱;u)`$ are

$`\mathrm{\Delta }(u)`$ $`=`$ $`\sqrt{(a^2+u)(1+u)(c^2+u)},`$
$`R(𝐱;u)`$ $`=`$ $`\sqrt{{\displaystyle \frac{x^2}{a^2+u}}+{\displaystyle \frac{y^2}{1+u}}+{\displaystyle \frac{z^2}{c^2+u}}}.`$
Hence, in the first approximation the number density is homeoidally striated, but the space-charge potential is not.

### IV.2 Determination of the structure in the second and higher approximations

The number density in the second approximation, $`n_2(𝐱)`$, follows upon substituting $`\mathrm{\Phi }_1(𝐱)`$ 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 $`\mathrm{\Phi }_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 $`\mathrm{\Phi }(𝐱)`$ 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 $`\mathrm{\Phi }_0`$, derived from Eq. (15), is easy to evaluate, and it enables fast, high-precision orbital integration. However, $`\mathrm{\Phi }_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 $`\mathrm{\Phi }_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 $`\mathrm{\Phi }(𝐱)`$, then there is neither computational benefit nor motivation for using $`\mathrm{\Phi }_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, $`\mathrm{\Phi }_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, $`\mathrm{\Phi }(𝐱)`$.

### V.1 Solution for the zeroth-order potential

Per Eq. (13), the minimum possible focusing strength corresponding to a spherically symmetric system is $`\mathrm{\Omega }=\mathrm{\Omega }_u=1/\sqrt{3}`$. Accordingly, we choose focusing strengths in keeping with the following labeling convention:

$$\mathrm{\Omega }_i=\frac{1}{\sqrt{3}}\left[1+10^{1i}10^{1I_{max}}\right],$$
(18)
with $`i=1,2,\mathrm{},I_{max}`$ and $`I_{max}=10`$. Choices of parameters $`[a,c;\mathrm{\Omega }_i]`$ then must be selected based on the constraint of Eq. (13) which bounds the parameter space of exact potentials. Note, however, that we can examine any desired geometry: oblate axisymmetric (for which $`a=1,c<1`$), prolate axisymmetric (for which $`a=c<1`$), and the full range of oblate-through-prolate triaxial systems. Hereafter we refer to “Case 1, Case 2,…” according to “$`i=1,2,\mathrm{}`$” in Eq. (18), respectively.

Plots of $`\mathrm{\Phi }_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 quasi-linear; however, it is manifestly nonlinear in the “Debye fall-off region” (henceforth called the “Debye tail”). Fig. 1 shows that the choices of $`\mathrm{\Omega }_i`$ per Eq. (18) span a wide range of space charge. At the one extreme, zero space charge, the density profile is gaussian. Then, the range of $`i`$ spans from small space charge ($`i=1`$) for which the density profile is approximately gaussian, through the fully space-charge-dominated, uniform beam ($`i=I_{max}=10`$) for which $`\mathrm{\Omega }=1/\sqrt{3}`$. Note that Case 5 ($`\mathrm{\Omega }1.0001/\sqrt{3}`$) represents “intermediate space charge”; the density falls off over a length scale comparable to that of the core.

### V.2 Methodology for orbital analysis

#### V.2.1 Samples of orbits, power spectra, and complexity

After choosing a parameter set $`[a,c;\mathrm{\Omega }_i]`$ and solving for $`\mathrm{\Phi }_0[R(𝐱)]`$, we began by generating 2000 initial coordinates uniformly spanning the volume occupied by the core and Debye tail of the density $`n_1[R(𝐱)]`$. The initial velocities were all chosen to be zero. Next, starting with these initial conditions, for every orbit we integrated the equations of motion, cf. Eq. (14), using $`\mathrm{\Phi }=\mathrm{\Phi }_0`$ for at least 100 orbital periods and, in most cases, for $`200`$ orbital periods. The integrations were done using a fifth-order Runge-Kutta algorithm ptvf  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 bennetin . 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 ptvf . 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 complexkandrup .

#### V.2.2 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 $`\chi `$ and the complexity $`n`$ are well-established, conventional measures of chaos tabor . A first piece of information for defining the criterion comes from plotting $`n`$ versus $`\chi `$ for all of the orbits. Fig. 2 provides an example. It shows that (a) $`n`$ increases approximately linearly with $`\chi `$, as one might expect since both quantities are measures of chaos, and (b) the regular orbits occupy a sharply defined region close to the origin of the $`n`$-vs.-$`\chi `$ plot. The borders of this “region of regularity” thereby offer three possible criteria for chaos: one involving only the Lyapunov exponent, one involving only the complexity, and one involving both quantities.

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 orbitref ; 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.-$`\chi `$ 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. Fig. 2(c)\]. This choice was also checked and confirmed by carefully inspecting hundreds of plots of individual orbits in several geometries.

Plots of two representative orbits, one regular and one chaotic, as well as their power spectra in the $`x`$-direction and their surfaces of section in the $`dz/dt\text{-vs.-}z`$ phase space, appear in Fig. 3. Both orbits have similar total energies and evolve in $`\mathrm{\Phi }_0`$ corresponding to Case 5 and $`a^2=1.0,c^2=0.25`$, a strongly oblate spheroid. The first orbit has $`n=8`$ and $`\chi =0.028t_D^1`$; the second orbit has $`n=178`$ and $`\chi =0.640t_D^1`$. Their differences are striking. One major difference concerns their power spectra. The regular orbit exhibits two very distinct frequencies in its spectrum, whereas the chaotic orbit features a near-continuum spanning a large number of frequencies. The surfaces of section are computed by recording $`dz/dt`$ and $`z`$ when $`x=y=0`$. The regular orbit is seen to stay within a localized region of the $`dz/dt`$-vs.-$`z`$ phase space, whereas the chaotic orbit largely fills a global area of phase space commensurate to its total energy.

Results of a numerical experiment that highlights the largest Lyapunov exponent, i.e., the rate of global chaotic mixing, appears in Fig. 4. Here, four initially localized clumps comprising 500 chaotic orbits evolve in the zeroth-order potential $`\mathrm{\Phi }_0`$ corresponding to Case 5 and $`a^2=0.5,c^2=1.5`$, a triaxial configuration. Each clump initially occupies a cube of size $`0.05^3`$ in the configuration space. The figure reveals that each clump mixes through a global region of phase space with an e-folding time comparable to a dynamical time $`t_D`$, taken 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.

### V.3 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 $`\mathrm{\Phi }_0`$ found from Eq. (17). Recall that this potential depends on the focusing strength $`\mathrm{\Omega }`$ and, through $`R(𝐱)`$, the scale lengths ($`a,c`$); we thus choose a specific parameter set $`[a,c;\mathrm{\Omega }]`$. Then we integrate 2000 initial conditions generated by uniformly sampling (1) the volume spanned by the configuration, and afterward (2) the volume spanned only by the Debye tail. Having established the criterion for chaos (complexity $`n>20`$), we identify and count the chaotic orbits in the respective sample, and we express this number as a percentage of the sample. This percentage may be viewed as an indication of the extent to which a given parameter set supports globally chaotic orbits. However, it should not be taken too literally in that the initial conditions are distributed uniformly through a volume; they are not weighted by the actual density distribution.

Our numerical experiments fall largely into two categories. Category I pertains to keeping $`\mathrm{\Omega }`$ fixed to its Case 5 value (intermediate space charge), i.e., $`i=5`$ in Eq. (18), and then varying $`a`$ and $`c`$ within the constraint of Eq. (13). Category II pertains to keeping $`c^2=0.5`$ fixed, and then varying $`\mathrm{\Omega }`$ and $`a`$.

Figures 5 and 6 depict the percentage of chaotic orbits for a portion of the Category I experiments. Also shown in the bottom panels of these figures are the initial conditions, projected onto the $`(x,z)`$-plane, corresponding to the axisymmetric configurations for which $`a^2=0.5`$. In addition, results from analogous experiments in exact potentials are plotted for comparison. The figures exhibit a number of features. First and foremost is the indication that most of the configurations support a considerable population of chaotic orbits. This is true even for axisymmetric configurations, but of course it is not true for spherically symmetric configurations in that their potentials are integrable. Second, orbits for which the initial conditions all lie within the Debye tail reflect a higher percentage of chaos than were they distributed through the entire configuration space. Third, prolate axisymmetric configurations (for which only one datum is shown in each of these figures) support little chaos. Fourth, for reasons to be elaborated in Sec. VI below, the exact potentials support less chaos than their zeroth-order counterparts.

Further analysis reveals that, for configurations in which they are present, essentially all of the chaotic orbits originate in the Debye tail. Figure 7 dramatically illustrates this finding. To assemble this figure, the initial conditions are sorted and binned into increments of $`R`$ spanning 0.1 units of length. Then, in each increment, the complexity $`n`$ of every orbit is computed, and the complexities are averaged to obtain $`n`$. The process is repeated for all of the increments, and then for many different parameter sets. The results, $`n(R)`$ versus $`R`$ for Case 5 with $`a^2=0.5`$ and several choices of $`c^2`$, comprise Fig. 7.

Figure 8 depicts the percentage of orbits for the Category II experiments for which the initial conditions are uniformly distributed over the volume spanned by the configuration. Cases corresponding to intermediate space charge would seem to support more chaos. This finding makes sense when juxtaposed against the limiting cases of zero space charge at one extreme and zero Debye tail, i.e., the uniform beam, at the other extreme. With zero space charge, only the linear forces of the external potential influence the particles. For the uniform beam the external and space-charge potentials cancel one another so that the particles move freely, apart from reflections at the boundary of the configuration. In both extremes the orbits are all regular, excepting billiard effects, if any, associated with shapes of boundary surfaces of uniform bunches. The figure also suggests that prolate, axisymmetric TE configurations support little chaos.

## 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 kandrup94 . Our survey of the parameter space of TE configurations centered on the use of $`\mathrm{\Phi }_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 $`\mathrm{\Phi }_0`$. For example, in the limit of distances very far from the centroid, the exact space-charge potential will approach spherical symmetry, whereas $`\mathrm{\Phi }_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 $`\mathrm{\Phi }(𝐱)`$, which is a fully three-dimensional partial differential equation (PDE), we chose a multigrid algorithm brandt . 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 ptvf .

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.3, 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 $`\mathrm{\Phi }_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(𝐱)]`$, the density derived from the exact space-charge potential $`\mathrm{\Phi }(𝐱)`$ is quasi-uniform over a larger volume and falls to small values over a shorter scale length. Accordingly, the configuration-space 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 $`\kappa `$ and $`\rho `$ of Eq. (6) take the form

$$\kappa =\frac{1}{2}\left[\mathrm{\Omega }^2\left(\frac{1}{a^2}+1+\frac{1}{c^2}\right)n(𝐱)\right],\rho =\frac{\sqrt{n^2(𝐱)n(𝐱)^2}}{\kappa \sqrt{2}}.$$
(19)
and Eq. (5) then yields the mixing rate $`\chi `$. 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 $`\chi `$ 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 $`\mathrm{\Phi }(𝐱)=\mathrm{\Phi }_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 $`11R25`$, 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, time-independent, 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 bohn02 . 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 kansid02 . 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 time-independent potentials. Results will be described in a forthcoming paper Nbody . 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.

###### Acknowledgements.

We thank Henry Kandrup for helpful discussions and for financially supporting one of us (IVS) through NSF Grant AST-0070809 in the early stages of this work. We also benefited from useful conversations with Rami Kishek. Most of this work was supported through U.S. Department of Education Grant G1A62056.


format: text

text:
# Bright and dark breathers in Fermi-Pasta-Ulam lattices

## I Introduction and model

Discrete breathers, also called intrinsic localized modes, are classical exact spatially localized time-periodic solutions which can be sustained by many nonlinear lattices (see Willis ; SiPage  for a review). In 1994 MacKay and Aubry Mackay  rigorously proved their existence in Hamiltonian lattices with anharmonic on-site potentials and weak coupling. Breathers are obtained by continuation from the uncoupled case in which trivial breathers exist (e.g. breathers with only one oscillator excited, the others being at rest). With the same technique, the existence of breathers was also proved for diatomic Fermi-Pasta-Ulam (FPU) chains Livi . In this model, two different masses alternate on the chain and are nonlinearly coupled to their nearest-neighbours via an interaction potential $`V`$. This result is valid for a large mass ratio, since light masses (next-nearest neighbours) are weakly coupled due to the presence of a heavy mass in between.

Unfortunately this method is not applicable to homogeneous FPU lattices, which do not possess an uncoupled limit in which trivial breathers exist. For some years the only existence result concerned the particular potentials $`V(x)=x^{2m}`$ with $`m2`$ Flach . FPU lattices seemed elusive to a rigorous mathematical treatment in spite of numerous approximate ST88 ; Pag90 ; Huang  and numerical Bickham ; Sandusky  studies, which have indicated the existence of discrete breathers in these systems.

Nevertheless, recent papers have presented rigorous existence proofs for discrete breathers in infinite FPU lattices. On the one hand, Aubry et al AKK01  have proved the existence of breathers with frequencies above the phonon spectrum, when $`V`$ is a strictly convex polynomial of degree $`4`$. These results are obtained via a variational method and apply in fact to higher dimensional generalizations of FPU lattices. Without additional assumptions on $`V`$, these results give only a partial information on breathers amplitudes. Under the additional condition that $`V`$ is even, Aubry et al prove the existence of breathers of arbitrarily small amplitudes in one-dimensional FPU lattices.

On the other hand, the existence (resp. non existence) of small amplitude breathers with frequencies slightly above the phonon band has been established when $`V`$ satisfies (resp. violates) a local hardening condition J01 . The proof is based on a centre manifold technique, which shows more generally that all small amplitude time-periodic solutions of the FPU system (including breathers) are determined by a two-dimensional map, provided their frequency lies near the phonon band edge. The discrete centre manifold method has been put in a general framework J03  and applied to other systems, namely diatomic FPU chains JN  (far from the uncoupled regime) and spin lattices N03 .

The aim of this work is to test numerically the range of validity of the centre manifold method J01  and to explore new phenomena, far from the small amplitude regime.

The one-dimensional FPU system is given by the equations:

$$\ddot{x}_n=V^{}(x_{n+1}x_n)V^{}(x_nx_{n1}),n,$$
(1)
where $`x_n`$ represents mass displacements from their equilibrium positions and $`V`$ is a smooth interaction potential satisfying $`V^{}(0)=0`$, $`V^{\prime \prime }(0)=1`$. In ref. J01 , the existence of small amplitude breathers (SAB) with frequencies $`\omega _b`$ slightly above the phonon band ($`\omega _b>2`$) is obtained for $`B>0`$, where

$$B=\frac{1}{2}V^{(4)}(0)(V^{(3)}(0))^2,$$
(2)
and their non-existence is proved for $`B<0`$ (see T72 , F96  for related results on the modulational instability of nonlinear normal modes and the tangent bifurcation of standing waves respectively).

The parameter $`B`$ can be interpreted as a hardening coefficient, since breathers with amplitude $`A0`$ have frequency $`\omega _b2+\frac{B}{8}A^2`$. Note that $`B`$ is slightly different from the classical hardening coefficient of an anharmonic potential ($`V`$ is hard if $`\frac{3}{5}V^{(4)}(0)(V^{(3)}(0))^2>0`$ and soft for $`\frac{3}{5}V^{(4)}(0)(V^{(3)}(0))^2<0`$).

In this paper we consider anharmonic potentials

$$V(u)=\frac{u^2}{2}+\frac{K_3}{3}u^3+\frac{K_4}{4}u^4$$
(3)
(at least one of the coefficients $`K_3`$, $`K_4`$ is nonzero). One has $`B=3K_44K_3^2`$ and thus the parameter $`B`$ is positive if $`|K_3|<\sqrt{3K_4}/2`$. In the case $`K_3=0`$ we have an even potential and the sign of $`B`$ coincides with the sign of the quartic coefficient $`K_4`$.

In section II we compute numerically SAB and find that they are very well fitted by leading order analytical expressions deduced from ref. J01 . In section III we consider a class of potentials for which numerical computations yield large amplitude breathers (LAB) with frequencies near the top of the phonon band. In this case there is an energy threshold for the existence of exact discrete breathers. The existence of an energy threshold for breather creation in higher-dimensional lattices is well known FKM ; W99 , but only a few one-dimensional examples have been given where such a threshold exists MW ; AKK01 ; Bernardo . Section III explores in more details the example of ref. Bernardo . Section IV deals with dark breathers, i.e. spatially modulated standing waves whose amplitude is constant at infinity and vanishes at the center of the chain. Leading order analytical expressions for small amplitude dark breathers (with frequencies inside the phonon band) fit the numerically computed dark breathers very well in the case of even potentials. In section V we analyze bright and dark breathers superposed on uniformly stretched or compressed static states. In each section the linear stability of the above mentioned solutions is numerically investigated.

## II Small amplitude breathers

### II.1 Summary of the theory

We consider time-periodic solutions of (1) having a given frequency $`\omega _b`$. We introduce the interaction force

$$y_n=V^{}(u_n)=u_n+K_3u_n^2+K_4u_n^3,$$
(4)
where

$$u_n=x_nx_{n1}.$$
(5)
The time-average of $`y_n`$ is independent of $`n`$ (integrate (1)), and is fixed to $`0`$ in reference J01 , as it is the case for all spatially localized solutions. Then, problem (1) leads to

$$\frac{d^2}{dt^2}(V^{})^1(y_n)=y_{n+1}2y_n+y_{n1},n.$$
(6)
The main results of J01  can be formulated as follows. For $`\omega _b2`$ (maximal phonon frequency), all small amplitude solutions $`y_n`$ (even in $`t`$) have the form

$$y_n(t)=b_n\mathrm{cos}(\omega _bt)+\phi (t,b_{n1},b_n,\omega _b),$$
(7)
where $`\phi `$ is a smooth time-periodic function with frequency $`\omega _b`$ and $`\genfrac{}{}{0pt}{}{\text{Sup }}{t[0,2\pi /\omega _b]}|\phi (t,b_{n1},b_n,\omega _b)|=O((|b_{n1}|+|b_n|)(|b_{n1}|+|b_n|+|\omega _b2|^2))`$. This result can be seen as an exact version of the rotating-wave approximation method ST88 , in the small amplitude limit. Note that $`\phi `$ can be computed at an arbitrary order (see J03  p. 51). Moreover, $`b_n`$ satisfies the second order nonlinear recurrence relation

$$b_{n+1}+2b_n+b_{n1}=\mu b_n+Bb_n^3+\text{h.o.t}$$
(8)
where $`\mu =\omega _b^241`$ and the coefficient $`B`$ is given in (2). Higher order terms are $`O(|b_n|((|b_{n1}|+|b_n|)^2+|\mu |)^2)`$ (note that the right side of (8) can be computed at an arbitrary order). Consequently, the problem of finding small amplitude time-periodic solutions of (6) reduces to the problem (8), which can be viewed as a two-dimensional map $`(b_n,b_{n1})(b_{n+1},b_n)`$.

For $`B>0`$ fixed and $`\mu >0`$ sufficiently small, the recurrence relation (8) has homoclinic solutions to $`0`$ satisfying $`\underset{n\pm \mathrm{}}{lim}b_n=0`$ (a proof has been given in J01 , J03 , using the invariance $`nn`$ of (8)). More precisely, there exist homoclinic solutions to $`0`$ denoted as $`\pm b_n^1`$, $`\pm b_n^2`$, having different symmetries $`b_{n1}^1=b_n^1`$, $`b_n^2=b_n^2`$.

Beyond this particular result, one can expect the existence of infinitely many homoclinic solutions to $`0`$ (not necessarily symmetric). Indeed, for $`B>0`$ the intersections of the stable and unstable manifolds of $`b_n=0`$ are generically transverse, which yields the existence of “homoclinic tangles” guck .

A formal way to understand why homoclinic orbits exist is the following. Setting $`\xi _n=(1)^nb_n`$ yields the recurrence relation

$$\xi _{n+1}+\xi _{n1}2\xi _n=\mu \xi _nB\xi _n^3+\text{h.o.t.}$$
(9)
Since $`0<\mu 1`$, equation (9) can be approximated by an integrable differential equation

$$v^{\prime \prime }=vv^3$$
(10)
with

$$\xi _n=\sqrt{\frac{\mu }{B}}v(n\sqrt{\mu }).$$
(11)
Equation (10) has the homoclinic solutions $`v(x)=\pm \sqrt{2}/\mathrm{cosh}(x+c)`$.

Using the above analysis, it has been proved that for $`B>0`$ and $`\omega _b2>0`$ sufficiently small, there exist SAB solutions of (6) given by equation (7) and homoclinic solutions of (8) (see J01 , J03 ). SAB have the form

$$y_n=(1)^n\xi _n\mathrm{cos}(\omega _bt)+\text{h.o.t.},$$
(12)
where $`\xi _n=(1)^nb_n`$ satisfies the recurrence relation (9) and $`|\xi _n|0`$ as $`n\pm \mathrm{}`$. The above homoclinic solutions $`b_n^1`$, $`b_n^2`$ correspond via equation (7) to SAB $`y_n^1`$, $`y_n^2`$ having different symmetries $`y_{n1}^1(t)=y_n^1(t+T_b/2)`$, $`y_n^2=y_n^2`$, where $`T_b=2\pi /\omega _b`$ is the breather period (see J03  p.58,59). Note that $`y_n(t+T_b/2)y_n(t)`$ in general due to higher order terms in (12) (however the egality holds if $`V`$ is even). Homoclinic solutions $`b_n^i`$ of (8) simply correspond to $`y_n^i(t+T_b/2)`$.

The exact solutions $`y_n^1`$, $`y_n^2`$ can be approximated at leading order using (10)-(11) with $`v(x)=\sqrt{2}/\mathrm{cosh}x`$. This yields

$$y_n^2(t)(1)^n\sqrt{\frac{2\mu }{B}}\frac{\mathrm{cos}\omega _bt}{\mathrm{cosh}(n\sqrt{\mu })},$$
(13)
and by symmetry we construct the second approximation

$$y_n^1(t)(1)^n\sqrt{\frac{2\mu }{B}}\frac{\mathrm{cos}\omega _bt}{\mathrm{cosh}((|n+1/2|1/2)\sqrt{\mu })}.$$
(14)
Under these approximations one has $`y_n^i(t+T_b/2)y_n^i(t)`$ and $`y_0^1(t)y_0^2(t)`$. In fact, the exact solutions satisfy $`y_n^i(t+T_b/2)=y_n^i(t)+O(|\mu |)`$ and $`y_0^1(t)=y_0^2(t)+O(|\mu |^{3/2})`$ as $`\mu 0^+`$.

Notice that approximation (13) can also be derived using multiscale expansions T72 . Moreover, these expressions also approximate $`u_n`$ since $`u_n=y_n+O(y_n^2)`$.

These calculations show that the maximum amplitude $`A\sqrt{\frac{2\mu }{B}}`$ of the SAB is $`O(|\omega _b2|^{1/2})`$ as $`w_b2^+`$. More precisely, breathers with amplitude $`A0`$ have frequency $`\omega _b2+\frac{B}{8}A^2`$.

Although SAB decay exponentially as $`|n|+\mathrm{}`$, their width is $`O(|\omega _b2|^{1/2})`$ and diverges as $`w_b2^+`$. As a consequence, if $`B>0`$ breathers exist for any small value of energy in the FPU system (1).

### II.2 Numerical continuation

We have computed these solutions numerically in order to check the range of validity of approximations (13), (14). We have fixed $`K_4=1`$ in eq. (3) and thus the parameter $`B`$ is positive if $`|K_3|<K_3^{}=\sqrt{3}/20.86`$.

We have performed our computations using a numerical scheme based on the concept of anti-continuous limit and path continuation with the Newton method MA . With this technique it is more convenient to use the difference displacements variables $`u_n`$. Indeed, with these new variables the dynamical equations become

$$\ddot{u}_n+2V^{}(u_n)C\left[V^{}(u_{n+1})+V^{}(u_{n1})\right]=0,n,$$
(15)
where $`C=1`$ for our system, but it makes possible to consider the system (15) as a Klein–Gordon system, with an on–site potential $`V`$ and an anticontinuous limit at $`C=0`$. Note that $`u_n`$ is one-to-one related to the forces $`y_n`$ at small amplitudes, because $`V^{}`$ is locally invertible since $`V^{\prime \prime }(0)0`$.

We use periodic boundary conditions $`u_{n+2p}(t)=u_n(t)`$. The periodicity is considered for an even number of sites in order that the maximum frequency of the linear phonons is exactly 2 as in the infinite lattice (this frequency corresponds to the phonon with wave number $`\pi `$).

We have obtained breathers with symmetries $`u_n(t)=u_n(t)`$ (site-centred mode) and $`u_n(t)=u_{n1}(t+T_b/2)`$ (bond-centred mode). Note that the coordinates transformation (5) produces an exchange of the symmetry properties between both modes, i.e. in the difference displacements variables the site-centred mode is the Page mode Pag90  and the bond-centred mode is the Sievers-Takeno mode ST88 .

Fig.1 shows the excellent agreement between the approximation (13) (dashed line) and an exact site-centred mode obtained numerically (circles) for $`\omega _b=2.01`$ ($`\mu 0.04`$).

We have continued the SAB obtained numerically for $`B>0`$ as $`\omega _b`$ goes away from the phonon band. In Fig.2 we have plotted the maximum amplitude of the forces sup$`|y_n|`$ (dash lines) versus $`\mu ^{1/2}`$ for two different values of the cubic coefficient, $`K_3=0.6`$ ($`B=1.56`$) and $`K_3=0.3`$ ($`B=2.64`$) from top to bottom, respectively. The maximum amplitude of the forces is an approximately linear function of $`\mu ^{1/2}`$ up to $`\mu 0.6`$. The continuous lines represent the maxima given by approximation (13). This approximation is not satisfactory for very small values of $`B`$ since higher order terms in (8) cannot be neglected when $`B0`$. More precisely, the smaller $`B`$ is, the smaller $`\mu `$ must be chosen in (13) in order to have a good approximation.

The maximum amplitude of the relative displacements sup$`|u_n|`$ is also an approximately linear function of $`\mu ^{1/2}`$ (see squares in fig.2). This was expected for small $`\mu `$ since $`u_n=y_n+O(y_n^2)`$, but it occurs far from the phonon band, at least until values of $`\mu 2.25`$ . In fact, we have checked that, surprisingly, expressions (13) and (14) fit very well the profile of the relative displacements $`u_n`$ even far from the phonon band as Fig. 3 shows. Note that the vibration amplitudes of the breather in Fig. 3 are quite large.

In Fig 4 we have plotted again the maximum amplitude of the relative displacements (squares), but now for an even potential ($`K_3=0`$) and larger frequencies, in order to compare approximation (14) with an approximation derived by Sievers and Takeno in the large amplitude regime ST88 . The Sievers-Takeno approximation is obtained using the rotating wave approximation (only the first Fourier component of solutions is taken into account) and Green’s function techniques. For the Sievers-Takeno mode and $`K_4=1`$ ($`B=3`$) it predicts

$$sup|u_n|\frac{2}{3}\sqrt{w_b^23}.$$
(16)
This amplitude-frequency relation corresponds to the dashed line in Fig 4. On the other hand, from eq. (14) it follows (continuous line)

$$sup|u_n|\sqrt{\frac{2(w_b^24)}{3}}.$$
(17)
Clearly we see that approximation (16) fits better very localized solutions with large amplitudes, whereas approximation (17) works better up to moderate amplitudes ($`\omega _b2.4`$ in this case). In this sense both approximations can be considered complementary.

In addition we have performed a Floquet analysis in order to study the linear stability of the SAB Aubry . In the symmetric potential case the result is well known Bickham : the Page mode is stable while the Sievers-Takeno mode has a harmonic instability (a pair of real eigenvalues $`\sigma ,\sigma ^1`$ close to $`1`$) that increases with the breather frequency. When a cubic term in the potential $`V`$ is introduced $`(K_30)`$, the situation is more complex. The Sievers-Takeno mode shows again a harmonic instability (see Fig 5a) but it also appears oscillatory instabilities. As Fig 5b shows for the Page mode, these oscillatory instabilities increase with $`|K_3|`$. Although most of them are size dependent effects due to the discretization of continuous spectrum MA98 , for the infinite system some of them might remain out of the unit circle and in this case SAB would be unstable. We shall not examine this problem here and leave it for future works.

## III Large amplitude breathers with an energy threshold

In this section we study breathers in the parameter region $`B<0`$. Using the recurrence relation (8), one can show that SAB with frequencies slightly above the phonon band do not exist for $`B<0`$ (see J01 , J03 ). However, large amplitude breathers (LAB) do exist in this parameter region AKK01 . More precisely, for $`V`$ strictly convex ($`K_4>0`$, $`|K_3|<\sqrt{3K_4}`$) and $`B<0`$ ($`\sqrt{3K_4}/2<|K_3|`$), breathers whose amplitudes do not tend to zero as $`w_b2^+`$ have been numerically obtained AKK01 ; Bernardo ; Kastner . As a consequence there is an energy threshold for breathers creation in these FPU systems. The criterion $`\sqrt{3K_4}/2<|K_3|<\sqrt{3K_4}`$ for the existence of an energy threshold originates from reference Bernardo , and follows by combining analytical results in J03  (non existence of small amplitude breathers) and AKK01  (existence of large amplitude breathers).

The existence of an energy threshold for breather creation has been formally analyzed by Flach, Kladko and MacKay FKM  for general Hamiltonian systems, in a different case when small amplitude breathers bifurcate at an edge of the phonon spectrum. An energy threshold exists above a critical lattice dimension $`d`$ (typically $`d=2`$) which depends on the nonlinear terms of the Hamiltonian. This analysis has been extended to a class of non-analytic Hamiltonian systems in Kastner . The existence of excitation thresholds has been rigorously analyzed by Weinstein W99  for the $`d`$-dimensional discrete non linear Schrödinger equation, using variational methods (see also MW  for formal variational arguments in the one-dimensional case).

In Fig. 6 we compare a SAB and a LAB with the same frequency $`w_b=2.01`$. As it is shown in Fig 7a the energy of a LAB (dash line) does not go to zero when approaching the phonon band edge. It lies always above a certain positive lower bound. In the case of SAB (continuous lines) the breather energy can be lowered to arbitrarily small values.

The dependence of the energy threshold on the asymmetry parameter $`K_3`$ is also shown in Fig. 7b, where we have fixed $`K_4=1`$. One observes that the energy threshold goes to $`0`$ as $`B0^{}`$ (for $`K_3\sqrt{3}/20,86`$). Moreover, LAB families also exist for nonconvex potentials ($`K_3<\sqrt{3}1,73`$), and their energy threshold increases with $`|K_3|`$. This shows that the criterion $`\sqrt{3K_4}/2<|K_3|<\sqrt{3K_4}`$ for the existence of an energy threshold is not optimal. We conjecturate that an energy threshold appears under the more general condition $`K_4>0`$, $`\sqrt{3K_4}/2<|K_3|`$ ($`B<0`$) in the one-dimensional FPU system.

Figure 8 shows the energy of a LAB (upper line) and its amplitude (lower line) versus frequency, close to the phonon band edge (the upper line is a close-up of figure 7a near $`\omega _b=2`$). The reason for the non monotonous behaviour of energy is that there are two competing tendencies: on the one hand amplitudes are decreasing with frequency and going to a nonzero lower bound, but on the other hand the width of the breather is increasing and more and more particles begin to oscillate as $`\omega _b`$ decreases.

We have found LAB with the same symmetries as SAB (Page and Sievers-Takeno modes). LAB have an exponential decay for $`\omega _b>2`$. As Fig. 9 shows, we have been able to fit the LAB profiles to the following expression in the case of the Page-mode:

$$u_n(t)\alpha (\omega _b)\sigma (\omega _b)^{|n|}\mathrm{cos}(\omega _bt)$$
(18)
where

$$\sigma (\omega _b)=1\frac{\omega _b^2}{2}+\frac{\omega _b}{2}(\omega _b^24)^{1/2}(1,0),$$
(19)
and $`\alpha (\omega _b)`$ is a fitting parameter. A similar expression corresponds to the Sievers-Takeno mode.

We have computed the decay rate $`|\sigma (\omega _b)|`$ as follows. The FPU system (1) can be formulated as a mapping

$$\left(\genfrac{}{}{0pt}{}{y_n}{y_{n+1}}\right)=F\left(\genfrac{}{}{0pt}{}{y_{n1}}{y_n}\right)$$
(20)
where the map

$$F\left(\genfrac{}{}{0pt}{}{y_{n1}}{y_n}\right)=\left(\genfrac{}{}{0pt}{}{y_n}{\frac{d^2}{dt^2}(V^{})^1(y_n)+2y_ny_{n1}}\right)$$
(21)
acts in a space of smooth time-periodic functions with frequency $`\omega _b>2`$. The linearized operator

$$DF(0)\left(\genfrac{}{}{0pt}{}{y_{n1}}{y_n}\right)=\left(\genfrac{}{}{0pt}{}{y_n}{(\frac{d^2}{dt^2}+2)y_ny_{n1}}\right)$$
(22)
has a purely hyperbolic spectrum and $`\sigma (\omega _b)`$, $`\sigma (\omega _b)^1`$ are the closest eigenvalues to $`1`$ (with $`\sigma (2)=1`$).

We have observed that LAB remain localized for $`\omega _b=2`$ (see pluses in figure 9) and thus approximation (18) does not work in this limit.

The stability properties of LAB are the same as SAB with non-even potentials, i.e. both modes show oscillatory instabilities and also a harmonic instability in the case of the Sievers-Takeno mode. As we mention in section II, most of the oscillatory instabilities are size dependent effects due to the discretization of continuous spectrum MA98 . However, a band analysis AKK01  suggests that, at least for some parameter values (close to $`\omega _b=3`$, $`K_3=1`$, $`K_4=1`$), some eigenvalues remain out of the unit circle for the infinite system and LAB are unstable. We have studied the evolution of Floquet eigenvalues of the Page mode for $`\omega _b=3`$, $`K_3=1`$ and $`K_4=1`$ when the system size becomes larger. In figure 10 one can see that one pair of eigenvalues seems to converge to values out of the unit circle for large system sizes, as stated in AKK01 .

## IV Dark Breathers

Using the recurrence relation (8), it has been shown that for $`B<0`$ and $`\mu <0`$ small enough ($`\omega _b`$ lies in the phonon band) there exist small amplitude heteroclinic solutions of (6) connecting two nonlinear normal modes $`y_n^\pm `$ at infinity, i.e. $`\underset{n\pm \mathrm{}}{lim}y_n(t)y_n^\pm (t)=0`$ (see J01 , J03 ).

We now briefly describe the nonlinear normal modes $`y_n^\pm `$. These solutions are spatially periodic ($`y_{n+2}^\pm (t)=y_n^\pm (t)`$), have $`0`$ time-average (the mean interaction force is $`0`$) and neighbours oscillate out of phase ($`y_{n+1}^\pm (t)=y_n^\pm (t+T_b/2)`$). Moreover, solutions $`y_n^\pm `$ are equal up to a half period phase shift ($`y_n^+(t)=y_n^{}(t+T_b/2)`$). Note that the corresponding displacement patterns $`x_n`$ are in general different from the zone-boundary modes $`x_n(t)=(1)^nf(t)`$ (also called binary oscillations). Indeed, in this latter case the time-average of relative displacements $`u_n`$ is $`0`$, whereas interaction forces have $`0`$ time-average in our case. However, both types of solutions coincide for even potentials $`V`$.

To obtain an approximate expression of the heteroclinic solutions, we can rewrite the recurrence relation (9) in the following form

$$\frac{\rho _{n+1}+\rho _{n1}2\rho _n}{\mu }=\rho _nB\rho _n^3+\text{h.o.t.},$$
(23)
using the variable change $`\xi _n=\sqrt{\mu }\rho _n`$. Equation (23) can be approximated by the differential equation:

$$v^{\prime \prime }=vBv^3$$
(24)
with $`\rho _n=v(n\sqrt{\mu })`$. Consequently, the system (6) has a family of small amplitude solutions which can be approximated by

$$y_n(t)(1)^n\sqrt{\frac{\mu }{B}}\mathrm{tanh}\left(\frac{n\sqrt{\mu }}{\sqrt{2}}\right)\mathrm{cos}(\omega _bt).$$
(25)
This type of solutions are called dark breathers dark  or antisymmetric kinks BFW97 . Typically the term breathers refers to bright breathers, that is a few excited units in a lattice. Here we have the opposite kind of localization: most oscillators are excited except one or a few of them.

Expression (25) approximates exact site-centred solutions of (6) satisfying $`y_n(t)=y_n(t+T_b/2)`$. Under our approximation one has $`y_n(t+T_b/2)y_n(t)`$, and in fact the exact solutions satisfy $`y_n(t+T_b/2)=y_n(t)+O(|\mu |)`$ (uniformly in $`n`$) as $`\mu 0^+`$.

Equation (6) has also exact bond-centred dark breather solutions satisfying $`y_n(t)=y_{n+1}(t)`$. We shall approximate them by

$$y_n(t)(1)^n\sqrt{\frac{\mu }{B}}\mathrm{tanh}\left(\frac{(n+1/2)\sqrt{\mu }}{\sqrt{2}}\right)\mathrm{cos}(\omega _bt).$$
(26)
Note that expressions (25), (26) also approximate $`u_n`$ since $`u_n=y_n+O(y_n^2)`$.

For computing dark breathers numerically we consider a lattice with an odd number of sites ($`51`$ lattice sites). We first choose an even potential ($`K_3=0`$) with $`K_4<0`$ (so that $`B<0`$). We solve equation (15) in Fourier space by the Newton method, fixing the time-average of $`u_0`$ to $`0`$ and using (25) as an initial condition. Dark breathers are computed for non-even potentials by continuation with respect to $`K_3`$, starting from $`K_3=0`$ (during the continuation the time-average of $`y_n`$ remains very small, due to the constraint on $`u_0`$).

Fig. 11 shows the excellent agreement between an exact site-centred solution obtained numerically (circles) and the approximation (25) (dashed line), for an even potential with $`K_4=1`$ ($`B=3`$) and a frequency $`\omega _b=1.9`$. For non even potentials (full circles) the agreement is not so good. The reason is clearly that the approximation of $`u_n`$ by $`y_n`$ with equation (25) has the symmetry property $`u_n(t+T_b/2)=u_n(t)`$ which is not true any more when $`K_30`$.

Similar results are obtained for bond-centred solutions (see Fig. 12).

For a fixed potential we can continue numerically these solutions decreasing $`\omega _b`$. As Fig. 13 shows, then the oscillation amplitudes increase and for even potentials expression (25) continues matching very well the relative displacements $`u_n`$ (circles), even far from the top of the phonon band. For non even potentials (squares) the agreement is also excellent close enough to the edge of the phonon band. For lower frequencies we observe discrepancies becoming larger.

We have also performed a Floquet analysis to investigate the linear stability of these family of solutions when we decrease the frequency.

On the one hand, the site-centered mode is stable close to the edge of the phonon band for even potentials as figure 14a shows. However, a cascade of harmonic instabilities appears for low enough frequencies (larger amplitudes). These instabilities are visible in Fig. 14a at $`\omega _b1.43`$. Before coming out the harmonic instabilities one observes an instability bubble at $`\omega _b1.53`$. Nevertheless, its magnitude decreases as the system size increases and thus it should not be relevant for the infinite system.

On the other hand, the bond-centered mode shows a harmonic instability (see figure 15a) that increases when the frequency decreases. One also observes another cascade of harmonic instabilities at $`\omega _b1.43`$.

Non-even potentials induce oscillatory instabilities for low enough frequencies (larger amplitudes), both for the site-centered and the bond-centered modes. Fig. 14b and fig. 15b show that this occurs for frequencies lower than $`1.85`$, for $`K_3=0.6`$ and $`K_4=1`$. Numerically we haven’t found these oscillatory instabilities near the phonon band edge, but we do not exclude that they might exist and be harder to detect.

The harmonic instability before mentioned for the bond-centered mode is also present for non even potentials, but is not visible in fig. 15b due to oscillatory instabilities. Note that frequency scales in figures a) and b) are different. For this reason in figures b) we do not see either the cascade of harmonic instabilities around $`\omega _b1.43`$. We have checked that they are also present with non-even potentials.

## V  Bright and dark breathers with an uniform stress

Breathers considered in section III for $`B<0`$ can be numerically continued at fixed frequency up to values $`B>0`$, bringing about a new type of solutions which satisfy $`lim_{n\pm \mathrm{}}u_n=c`$, where $`c`$ is a nonzero constant. As Fig. 16 shows, these new type of solutions can coexist with SAB for the same values of the parameters and can be seen as time-periodic oscillations around static solutions of (1) given by $`x_n=cn`$. The lattice is uniformly stretched at infinity for $`c>0`$, and uniformly compressed for $`c<0`$. Such breather solutions with a constant static strain $`c<0`$ have been also observed by Sandusky and Page in the case $`B<0`$ (see SanduskyP  p.876).

This section provides analytical and numerical results concerning breather solutions with a constant static strain. The case of dark breathers with an uniform stress is also examined.

In the sequel we assume $`V^{\prime \prime }(c)>0`$. We make the change of variable $`x_n(t)=cn+\stackrel{~}{x}_n(t\sqrt{V^{\prime \prime }(c)})`$ where $`\stackrel{~}{x}_n1`$ and introduce the modified potential

$$\stackrel{~}{V}(u)=(V(c+u)V^{}(c)u)/V^{\prime \prime }(c).$$
(27)
Equation (1) leads to

$$\frac{d^2}{dt^2}\stackrel{~}{x}_n=\stackrel{~}{V}^{}(\stackrel{~}{x}_{n+1}\stackrel{~}{x}_n)\stackrel{~}{V}^{}(\stackrel{~}{x}_n\stackrel{~}{x}_{n1}),n,$$
(28)
with $`\stackrel{~}{V}^{}(0)=0`$ and $`\stackrel{~}{V}^{\prime \prime }(0)=1`$. One can see that equation (28) is exactly equation (1) where $`V`$ has been changed to $`\stackrel{~}{V}`$. Consequently, the results of the above sections (both analytical and numerical) readily apply to equation (28).

Returning to the original variables, the following analytical results for the FPU system (1) follow from the analysis of section II-A. Let us fix $`c`$ and define

$$B(c)=\frac{1}{2}V^{\prime \prime }(c)V^{(4)}(c)(V^{(3)}(c))^2.$$
(29)
If $`B(c)>0`$ and $`\omega _b2\sqrt{V^{\prime \prime }(c)}>0`$ is sufficiently small, the FPU system has a family of exact solutions which can be approximated by

$$u_n(t)c+(1)^n\sqrt{\frac{2\mu }{B(c)}}\frac{\mathrm{cos}(\omega _bt)}{\mathrm{cosh}(n\sqrt{\mu /V^{\prime \prime }(c)})},$$
(30)
where $`\mu =\omega _b^24V^{\prime \prime }(c)`$, $`0<\mu 1`$. These solutions are time-periodic (with frequency $`\omega _b`$) and their oscillating part is spatially localized. As $`n\pm \mathrm{}`$ they converge towards a uniformly stressed static state (stretched or compressed depending on $`c`$). Note that their lower frequency lies inside the phonon band if $`V^{\prime \prime }(c)<1`$, and above for $`V^{\prime \prime }(c)>1`$. The exact solutions approximated by (30) are site-centred, i.e. they satisfy $`u_n(t)=u_n(t)`$. There exists also exact bond-centred solutions satisfying $`u_n(t)=u_{n1}(t+\pi /\omega _b)`$. In the sequel we shall denote these solutions as breathers with an uniform stress.

As we previously mentioned, the numerical study of these solutions is already contained in section II-B, which applies to the renormalized displacements $`\stackrel{~}{x}_n`$ satisfying (28). For the polynomial potential (3), breathers with an uniform stress have the same stability properties as the SAB corresponding to the modified potential (27). In particular, oscillatory instabilities are observed for $`K_3=0`$ ($`V`$ is even), since the modified potential $`\stackrel{~}{V}`$ is non-even in this case. On the contrary, if $`K_30`$, $`K_4>K_3^2/6`$ then the modified potential $`\stackrel{~}{V}`$ is even for $`c=\frac{K_3}{3K_4}`$ (and the conditions $`V^{\prime \prime }(c)>0`$, $`B(c)>0`$ are realized). In this case oscillatory instabilities are not observed for the SAB with specific static strain $`c`$.

If $`B(c)<0`$ and $`K_4>0`$, we numerically obtain large amplitude solutions of the same type, whose oscillating parts do not vanish as $`w_b2\sqrt{V^{\prime \prime }(c)}`$ (with $`w_b>2\sqrt{V^{\prime \prime }(c)}`$ ). This result follows directly from the numerical study of section III, as well as the stability properties of these solutions. These LAB with an uniform stress have been previously observed by Sandusky and Page (see SanduskyP  p.874-876), for the interaction potential (3) with $`B<0`$. In their examples $`V^{\prime \prime }(c)>1`$ (since $`c<0`$, $`K_3<0`$, $`K_4>0`$) and thus the lower frequencies of these LAB families lie strictly above the phonon band. This phenomenon can be observed in fig. 8 of reference SanduskyP , which provides frequency vs amplitude plots of such solutions (note that $`c`$ may not be constant on each curve).

If $`B(c)<0`$ and $`\omega _b2\sqrt{V^{\prime \prime }(c)}<0`$ is sufficiently small, the FPU system has a family of exact solutions which can be approximated by

$$u_n(t)c+(1)^n\sqrt{\frac{\mu }{B(c)}}\mathrm{tanh}\left(\frac{n\sqrt{\mu }}{\sqrt{2V^{\prime \prime }(c)}}\right)\mathrm{cos}(\omega _bt),$$
(31)
where $`\mu =\omega _b^24V^{\prime \prime }(c)<0`$, $`|\mu |1`$. These solutions are time-periodic (with frequency $`\omega _b`$) and consists in a dark breather superposed with a uniformly stressed static state. Note that their higher frequency lies above the phonon band if $`V^{\prime \prime }(c)>1`$. These solutions are site-centred ($`u_n(t)=u_n(t+\pi /\omega _b)`$), but their bond-centred analogues also exist. They have the same stability properties as dark breathers in section IV, for the modified potential (27).

## VI Conclusions

In this work we have analyzed both numerically and analytically a rich variety of nonlinear solutions in FPU lattices. There are four main results.

On the one hand, we have numerically explored the range of validity of a recent breathers (bright or dark) existence proof J01  for one-dimensional FPU lattices. Our computations have been carried out with polynomial interaction potentials of degree $`4`$. For $`B>0`$, we have found that the maximum amplitude of bright breathers in the relative displacement variables is an approximately linear function of $`\mu ^{1/2}`$ even far from the phonon band and, in fact, the approximate expression (14) describes very well the profile of the breathers when their amplitudes are large. For $`B<0`$ and $`K_4<0`$, we find in the same way that dark breathers can be very well described by approximate analytical expressions.

On the other hand, we have found numerically for $`B<0`$ large amplitude bright breathers having frequencies arbitrary close to the top of the phonon band, which are out of range of the local analysis J01 . This breather family exhibits an energy threshold, which is a rarely observed phenomenon in one-dimensional lattices.

Moreover, with respect to stability properties we have checked that non-even potentials induce oscillatory instabilities in both bright and dark cases. In the case of even potentials, the Page mode is linearly stable while the Sievers-Takeno mode is unstable. We have also found that bond-centred dark breathers are unstable. The site-center dark breather is unstable for low enough frequencies and linearly stable with an even potential near the edge of the phonon band.

At last, we have analyzed small amplitude bright and dark breathers superposed to uniformly stressed static states $`x_n=cn`$. In particular, the lowest frequency of bright breathers is inside the phonon band if $`0<V^{\prime \prime }(c)<1`$, and the highest frequency of dark breathers is above the phonon band if $`V^{\prime \prime }(c)>1`$. We have obtained local conditions on $`V`$ for their existence as well as approximate analytical expressions.

###### Acknowledgements.

This work has been supported by the European Union under the RTN project LOCNET, HPRN–CT–1999–00163. J. Cuevas acknowledges a FPDI grant from la Junta de Andalucía. G. James acknowledges stimulating discussions with R. MacKay.