# Chaotic Dynamics Driven by Particle-Core Interactions

Konstantin Batygin<EMAIL_ADDRESS>author to whom correspondence should
be addressed. Division of Geological and Planetary Sciences California
Institute of Technology, Pasadena, CA 91125, USA  Yuri K. Batygin
<EMAIL_ADDRESS>Los Alamos National Laboratory, Los Alamos, NM 87545, USA

###### Abstract

High-intensity beams in modern linacs are frequently encircled by diffuse
halos, which drive sustained particle losses and result in gradual degradation
of accelerating structures. In large part, the growth of halos is facilitated
by internal space-charge forces within the beams, and detailed
characterization of this process constitutes an active area of ongoing
research. A partial understanding of dynamics that ensue within space-charge
dominated beams is presented by the particle-core interaction paradigm – a
mathematical model wherein single particle dynamics, subject to the collective
potential of the core, are treated as a proxy for the broader behavior of the
beam. In this work, we investigate the conditions for the onset of large-scale
chaos within the framework of this model, and demonstrate that the propensity
towards stochastic evolution is strongly dependent upon the charge
distribution of the beam. In particular, we show that while particle motion
within a uniformly charged beam is dominantly regular, rapid deterministic
chaos readily arises within space-charge dominated Gaussian beams.
Importantly, we find that for sufficiently high values of the beam’s space
charge and beam pulsation amplitude, enhanced chaotic mixing between the core
and the halo can lead to an enhanced radial diffusion of charged particles. We
explain our results from analytic grounds by demonstrating that chaotic motion
is driven by the intersection of two principal resonances of the system, and
derive the relevant overlap conditions. Additionally, our analysis illuminates
a close connection between the mathematical formulation of the particle-core
interaction model and the Andoyer family of integrable Hamiltonians.

##  I Introduction

An important and long-standing challenge to the performance of high-intensity
particle accelerators is presented by the emergence of halos around space-
charge dominated beams. Simply put, a beam halo is comprised of a tenuous
population of charged particles that fill a phase-space volume that is more
extended than that occupied by the beam’s core. Within a halo, particles can
execute long-range radial excursions, leading to their eventual absorption by
the accelerating structure. Correspondingly, this process limits the intensity
of the beam and results in radio-activation of the linac[1].

Although the aforementioned practical issues are broadly appreciated [6, 7,
8], efforts to suppress halo formation through beam collimation cannot resolve
the problem completely, as halos have a distinct tendency to regenerate over
relatively short timescales. Fundamentally, this propensity towards
reemergence indicates that intrinsic self-excitation of the beam through non-
linear space-charge forces plays an appreciable, if not dominant, role in the
process of halo formation [9, 10]. In this vein, one potential suppression
mechanism for beam halo formation employs nonlinear focusing, as originally
proposed in [2, 3] and further developed in [4, 5]. Nevertheless, a detailed
understanding of this halo-formation mechanism remains incomplete, and
constitutes an active area of investigation [11, 12, 13, 14].

One illustrative approach towards characterizing phase-space growth of space-
charge-dominated beams rests on the particle-core interaction model of halo
formation [9, 7]. Within the context of this picture, the evolution of a
single test-particle – subject to the collective potential of a breathing
axisymmetric beam – is treated as a proxy for typical dynamics exhibited by
charged particles. More specifically, the radial oscillations of the beam are
taken to stem from an integrable Hamiltonian, while the motion of the test-
particle is governed by a non-autonomous Hamiltonian with only a single degree
of freedom [15]. For typical parameters, this model yields a pulsating core
that is encircled by a separatrix, which can transport particles from the
outer edge of the core to a considerably larger radial distance, thereby
creating a diffuse halo.

Despite the particle-core interaction model’s apparent lack of complexity, it
has been shown to readily yield non-trivial behavior [7, 16], and has
displayed a noteworthy degree of agreement with experimental data [17].
Moreover, the model’s capacity to capture the transition from ordered, quasi-
periodic motion to deterministic chaos renders it an ideal paradigm for
exploring halo formation through stochastic transport of particles. We note,
however, that even though the appearance of nonlinear resonances and the
materialization of chaotic bands in phase-space have been routinely observed
in numerical simulations of the particle-core interaction model, the phase-
space attributes of these resonances, the specific conditions for the onset of
large-scale chaos, as well as the dependence of the beam’s stochastic behavior
on its radial charge distribution remain to be quantified from analytic
grounds.

Carrying out this analysis is the primary purpose of this paper. More
specifically, a key aim of this paper is to employ canonical perturbation
theory to systematically explore two (uniform and Gaussian) variants of the
particle-core interaction model, with an eye towards identifying the principal
resonant dynamics that facilitate halo formation. The remainder of the
manuscript is organized as follows. We begin our analysis in section II by
reviewing particle motion within a uniformly-charged beam. We then explicitly
identify the primary harmonic, responsible for large-amplitude radial
excursions, through time-averaging of the Hamiltonian. With a simplified model
in hand, we highlight a hitherto-overlooked connection between the particle-
core interaction model and the so-called “second fundamental model for
resonance” used extensively in celestial mechanics [18]. Subsequently, we
derive the associated resonance proximity parameter in terms of the beam’s
current and pulsation amplitude.

In section III, we consider the particle-core interaction model with a
Gaussian charge distribution, and repeat our perturbative analysis, this time
carrying out the time-averaging procedure semi-analytically. By comparing our
results with stroboscopic surfaces of section, we demonstrate that the beam
core of a Gaussian variant of the model is markedly more susceptible to
chaotic mixing, and outline a criterion for the onset of chaos through an
analogy with a modulated pendulum. We further demonstrate that resonant
excursion of particles to large radii can be sourced from deep within the
core. Finally, we confirm our results with a suite of numerical simulations,
and identify a parameter range where halo formation is strongly enhanced by
long-range stochastic diffusion. We briefly summarize our findings in section
IV.

##  II A Uniformly Charged Beam

The mathematical formulation of the particle-core interaction model amounts to
two second order (equivalently, two pairs of first order), coupled non-linear
ODEs: one for the radial pulsation of the beam core and the second for the
radial motion of the test-particle. As already alluded to in the previous
section, however, the coupling between these equations is uni-directional,
such that oscillations of the core around its equilibrium radius, $R$,
modulate the particle’s phase-space trajectory, but not vice-versa.
Accordingly, let us begin by quantifying the core’s radial breathing as a
function of time.

###  II.1 Pulsation of the Core

The dynamics of the beam core’s radius is governed by the Hamiltonian [9]:

$\displaystyle\mathcal{H}_{\rm{c}}=\frac{p_{r}^{2}}{2}+\frac{r^{2}}{2}+\frac{1}{2\,(1+b)}\frac{1}{r^{2}}-\frac{b\,\ln(r)}{1+b},$
(1)

where $r=x_{\rm{c}}/R$ is a dimensionless core radius and $p_{r}=\dot{r}$ is
the conjugate momentum. Strictly speaking, this mono-dimensional description
of the beam’s evolution entails axial symmetry, and is valid for an axially-
symmetric beam in a continuous solenoidal focusing channel. In a quadrupole
channel, on the other hand, equation (1) describes the evolution of a smoothed
envelope, averaged in both $\hat{x}$\- and $\hat{y}$-directions, wherein
small-amplitude fast oscillations are filtered out. We note that although the
inclusion of 2D or 3D pulsation (see e.g., [19]) can result in additional
higher-order dynamics, such effects are by no means central to the problem at
hand.

The Hamiltonian is parameterized the dimensionless space-charge [16]

$\displaystyle
b=\frac{2\,m\,c}{p}\frac{I}{I_{\rm{c}}}\frac{R^{2}}{\epsilon^{2}}$ (2)

where $m$ is the particle mass, $c$ is the speed of light, $p$ is the particle
momentum, $I$ is the beam current, $I_{\rm{c}}=4\pi\,\epsilon_{0}\,m\,c^{3}/q$
is the characteristic beam current, $\epsilon_{0}$ is the permittivity of free
space, and $\epsilon$ is the normalized beam emittance. Note that instead of
$b$, Wangler et al. [7] parameterize $\mathcal{H}_{\rm{c}}$ in terms of the
in-core space-charge tune-depression ratio

$\displaystyle\frac{\sigma}{\sigma_{0}}=\frac{1}{\sqrt{1+b}},$ (3)

where $\sigma$ and $\sigma_{0}$ are the phase advances of transverse particle
oscillations per focusing period, with and without space-charge forces,
respectively. Moreover, notice that because $\mathcal{H}_{\rm{c}}$ is
expressed in terms of dimensionless variables, the unit of time must also be
appropriately non-dimensional, and here, $t=\tau\,\Omega$, where $\Omega$ is
the frequency of particle oscillation in a uniform focusing channel in absence
of space-charge, and $\tau$ is the dimensional time.

The Hamiltonian (1) possesses a global minimum – which corresponds to an
elliptic equilibrium – at $(r,p_{r})=(1,0)$. Since we are primarily interested
in variations of ($r,p_{r})$ around this fixed point, it is sensible to change
variables to $\rho=r-1,p_{\rho}=p_{r}$ (addition or subtraction of a constant
to either the coordinate or the momentum always corresponds to a canonical
transformation [20]), and expand the Hamiltonian around the origin. Carrying
out the power-series expansion to cubic order in $\rho$, we have:

$\displaystyle\mathcal{H}_{\rm{c}}$
$\displaystyle\approx\frac{p_{\rho}^{2}}{2}+\frac{2\,(2+b)}{1+b}\,\frac{\rho^{2}}{2}-\frac{(6+b)}{1+b}\,\frac{\rho^{3}}{3}$
$\displaystyle=\frac{p_{\rho}^{2}}{2}+\omega^{2}\,\frac{\rho^{2}}{2}-\zeta^{2}\,\frac{\rho^{3}}{3},$
(4)

where the characteristic frequency of core oscillation is
$\omega=\sqrt{2\,(2+b)/(1+b)}$, $\zeta=\sqrt{(6+b)/(1+b)}$, and we have
dropped a dynamically irrelevant constant term from the expression.

Figure 1: Anharmonic approximation to pulsation of the core. The background
contour plot corresponds to the level contours of Hamiltonian (1). Yellow
curves are given by the approximate anharmoinc solution (5), for oscillation
amplitudes of $\Delta=0.1,0.2,0.3,$ and $0.4$. In the considered range of
parameter space, equation (5) provides an adequate approximation for the
phase-space portrait of Hamiltonian (1).

Hamiltonian (4) corresponds to an anharmonic oscillator. Although a closed-
form solution to the resulting equations of motion can be expressed in terms
of elliptic integrals [21], for our purposes, an approximate solution for
$\rho$ expressed as a rudimentary function of $t$ will be sufficient. Such a
solution is readily obtained from classical perturbation theory [22, 23] and
has the form (see appendix for details):

$\displaystyle\rho$
$\displaystyle\approx\Delta\cos(\omega\,t+\varphi)+\frac{\Delta^{2}\,\zeta^{2}}{2\,\omega^{2}}$
$\displaystyle-\frac{\Delta^{2}\,\zeta^{2}}{6\,\omega^{2}}\,\cos(2\,(\omega\,t+\varphi)),$
(5)

where $\Delta$ is the amplitude of oscillation, and $\varphi$ is an arbitrary
phase. For consistency with previous studies [11, 7, 16], we set $\varphi=\pi$
such that $r$ is minimal and $p_{r}$ is null at $t=0$.

The solution (5) is shown on a $(\rho,p_{\rho})$ plane in Figure (1) for a
series of values of $\Delta\in(0,0.4)$ and $b\in(0.1,10)$ – a physically
relevant parameter range that we will adopt for the remainder of the paper. We
note that these bounds on $\Delta$ and $b$ are motivated by experiment: even
the most high-intensity machines in operation are characterized by $b\lesssim
10$, while beams with oscillation amplitude in excess of $\Delta\gtrsim 0.4$
typically suffer immense losses, rendering them uninteresting from a practical
standpoint. On the same figure, level curves of Hamiltonian (1) are shown as
the background contour map. The agreement between the approximate series
solution (5) and the evolution entailed by $\mathcal{H}_{\rm{c}}$ is
satisfactory for $\Delta\lesssim 0.4$ and improves for higher values of $b$.
In fact, this dependence can be expected, since $(\zeta/\omega)^{2}$ ranges
from 3/2 for $b\rightarrow 0$ to $1/2$ for $b\gg 1$, implying that higher
order terms are slightly suppressed for space-charge dominated beams.
Accordingly, for the remainder of the manuscript, we will retain equation (5)
as an acceptable approximation, but note that while the physical meaning of
$r$ denotes the outer edge of the beam in case of a uniformly charged core, it
corresponds to twice the rms radius of the beam in the case of a Gaussian
charge distribution.

###  II.2 Particle Dynamics

With the time-dependence of $r$ specified, let us now review the dynamics of
an embedded particle from a perturbative perspective. In essence, the goals of
the following calculation are two-fold. First, we aim to outline the relevant
approximation scheme that we will revisit in the next section. Second, by
developing an approximate model, we seek to contextualize the phase-space
portrait of the particle-core interaction model from purely analytic grounds.

####  II.2.1 An Integrable Approximation

The motion of a particle subject to the space-charge forces of a uniformly
charged beam core is governed by the following Hamiltonian [9, 7, 16]:

$\displaystyle\mathcal{H}_{\rm{p}}=\frac{p_{u}^{2}}{2}+\frac{u^{2}}{2}-\frac{b}{1+b}\,\times$
$\displaystyle\begin{dcases*}\frac{u^{2}}{2\,r^{2}}&$|u|<r$\\\
\bigg{(}\frac{1}{2}+\ln(|u|(1+b))-\ln(r(1+b))\bigg{)}&$|u|\geqslant
r,$\end{dcases*}$ (6)

where $u$ and $p_{u}$ refer to the particle’s position and momentum,
respectively. Because $\mathcal{H}_{\rm{p}}$ exhibits explicit dependence on
time through $r$, it is not integrable, and analytic solutions for $(u,p_{u})$
cannot exist. Conversely, numerical experiments of dynamics jointly
facilitated by Hamiltonians (1) and (6) are plentiful in the literature [6,
10, 11, 7, 15], and need not be rehashed at this stage. As a result, let us
instead derive an integrable approximation to Hamiltonian (6), leaving
simulations for later. In spirit, our analysis will follow the same general
theme as the work of Gluckstern [9] and Batygin [16], although our approach to
perturbation theory will be different, and will illuminate the criteria under
which the averaging process is warranted.

As a starting point, we follow Gluckstern [9], and replace the piecewise-
continuous function in equation (6) with a smooth quartic polynomial:

$\displaystyle\mathcal{H}_{\rm{p}}\approx\frac{p_{u}^{2}}{2}+\frac{u^{2}}{2}-\frac{b}{1+b}\bigg{(}\frac{u^{2}}{2\,r^{2}}-\frac{u^{4}}{16}\bigg{)}.$
(7)

Before proceeding further, we caution that although this polynomial fit has
been suitably employed within the literature, it constitutes a relatively
crude approximation for the potential function outlined in equation (6). In a
purely quantitative sense, the fractional error entailed by equation (7) is on
the order of $10-20\%$ for $\Delta=0$, and can even exceed unity for
$\Delta=0.4$ and $|u|\gtrsim 1+\Delta$. From a qualitative perspective, it is
further important to note that a strictly linear field interior to $|u|<r$ in
equation (6) is replaced by a weakly non-linear one in expression (7). Despite
these shortcoming, estimates obtained within the context of this analytic
model provide a satisfactory match to the results of numerical experiments
(e.g., the radial extent of the resonant separatrix is reproduced to within
$\sim 10\%$, etc; [16]). Nevertheless, these approximations are sufficiently
pronounced for us to emphasize that any results derived from Hamiltonian (7) –
including the reduction of the problem to the Andoyer model outlined in
section II.2.2 below – should not be mistaken for a quantitatively accurate
description of particle dynamics. Rather, the key goal of this analysis is to
provide an intelligible “blueprint” for particle dynamics, within the context
of which the topological character of the phase-space portrait can be
interpreted.

Returning to the approximate form of the Hamiltonian, we substitute equation
(5) for $r=1+\rho$ and expand the resulting expression as a power-series in
$\Delta$. For the purposes of our analysis, it suffices to truncate this
expansion at linear order in $\Delta$ for two reasons: first, retention of
higher-order terms does not generate any slowly-varying harmonics and
therefore does not contribute appreciably to the description of the system’s
principal resonance; second, the polynomial approximation for the piecewise-
continuous potential employed in equation (7) is only valid for small values
of $\Delta$, meaning that a more precise representation of $(1/r^{2})$ is
unlikely to yield more accurate results.

Figure 2: Phase-space portraits of the Andoyer model, for a series of values
of $\delta$. Equilibria of the Hamiltonian are shown with yellow dots, while
the separatrix is shown with a black curve, where present. The topological
character of the Andoyer Hamiltonian is fully determined by the resonance
proximity parameter, $\delta$. Notably, in presence of a homoclinic curve,
elliptic stability of the origin – which translates to the existence of an
inner circulation region of the phase-space – necessitates that $\delta>0$.

With these approximations in place, the Hamiltonian takes the form:

$\displaystyle\mathcal{H}_{\rm{p}}$
$\displaystyle\approx\frac{p_{u}^{2}}{2}+\bigg{(}\frac{1}{1+b}\bigg{)}\frac{u^{2}}{2}+\bigg{(}\frac{b}{1+b}\bigg{)}\frac{u^{4}}{16}$
$\displaystyle-\bigg{(}\frac{b\,\Delta}{1+b}\bigg{)}\,u^{2}\,\cos(\omega\,t).$
(8)

To proceed further, we introduce action-angle coordinates $(\psi,\Psi)$,
implicitly defined in terms of $(u,p_{u})$:

$\displaystyle u=\sqrt{2\,\Psi}\,\sin(\psi)$ $\displaystyle
p_{u}=\sqrt{2\,\Psi}\,\cos(\psi).$ (9)

To be precise, this change of variables simply amounts to a transformation
from cartesian to canonical polar coordinates, and stems from a type-1
generating function of the form
$\mathcal{F}_{1}=-(1/2)\,u^{2}\,\tan^{-1}(\psi)$. After some rearrangement,
equation (8) becomes:

$\displaystyle\mathcal{H}_{\rm{p}}$
$\displaystyle\approx\Psi-\bigg{(}\frac{b}{1+b}\bigg{)}\,\Psi\,\sin^{2}(\psi)+\bigg{(}\frac{b}{4\,(1+b)}\bigg{)}\,\Psi^{2}\,\sin^{4}(\psi)$
$\displaystyle-\frac{b\,\Delta}{1+b}\bigg{(}\Psi\,\cos(\omega\,t)-\frac{1}{2}\Psi\,\cos(2\,\psi-\omega\,t)$
$\displaystyle-\frac{1}{2}\Psi\,\cos(2\,\psi+\omega\,t)\bigg{)}.$ (10)

In order to identify the slowly-varying harmonic in equation (10), let us
consider $\mathcal{H}_{\rm{p}}$ in the limit of vanishing $b$ and $\Delta$. In
this regime, it is trivial to see that
$\dot{\psi}=\partial\mathcal{H}_{\rm{p}}/\partial\Psi\rightarrow 1$.
Simultaneously, equation (4) dictates that $\omega=2$ for $b=0$. This means
that the harmonic $\cos(2\,\psi-\omega\,t)$ will be slowly varying, and the
associated angle must correspond to the principal resonance of the system.
With this notion in mind, let us convert to a new set of action-angle
variables $(\phi,\Phi)$ that reflect this property of the dynamics.

To accomplish this, we extend the phase-space of the system, by adding a dummy
action, $\mathcal{T}$, conjugate to time [24], such that the new Hamiltonian
becomes:

$\displaystyle\tilde{\mathcal{H}}_{\rm{p}}=\mathcal{H}_{\rm{p}}+\mathcal{T}.$
(11)

Next, we employ a canonical contact transformation of variables, stemming from
the type-2 generating function:

$\displaystyle\mathcal{F}_{2}=(\psi-\omega\,t/2)\,\Phi+(t)\,\Xi.$ (12)

The transformation equations then trivially yield:

$\displaystyle\phi=\frac{\partial\,\mathcal{F}_{2}}{\partial\,\Phi}=\psi-\frac{\omega}{2}\,t$
$\displaystyle\Psi=\frac{\partial\,\mathcal{F}_{2}}{\partial\,\phi}=\Phi$
$\displaystyle\xi=\frac{\partial\,\mathcal{F}_{2}}{\partial\,\Xi}=t$
$\displaystyle\mathcal{T}=\frac{\partial\,\mathcal{F}_{2}}{\partial\,t}=\Xi-\frac{\omega}{2}\,\Phi.$
(13)

In preparation for canonical averaging, let us write the transformed
Hamiltonian $\tilde{\mathcal{H}}_{\rm{p}}$ as
$\tilde{\mathcal{H}}_{\rm{p}}=\check{\mathcal{H}}_{\rm{p}}+\hat{\mathcal{H}}_{\rm{p}}$,
where

$\displaystyle\check{\mathcal{H}}_{\rm{p}}=\bigg{(}\frac{1}{2}\bigg{(}\frac{2+b}{1+b}\bigg{)}-\frac{\omega}{2}\bigg{)}\,\Phi+\frac{3\,b}{32\,(1+b)}\,\Phi^{2}+\Xi$
(14)

represents a trivially integrable kernel of the dynamics, and

$\displaystyle\hat{\mathcal{H}}_{\rm{p}}$
$\displaystyle=b\,\Phi\,\big{(}4\,(4-\Phi)\cos(2\,\phi+\omega\,\xi)+\Phi\,\cos(4\,\phi+2\,\omega\,\xi)$
$\displaystyle+16\,\Delta\,(\cos(2\,(\phi+\omega\,\xi))-2\,\cos(\omega\,\xi)$
$\displaystyle+\cos(2\,\phi))\big{)}/(32\,(1+b))$ (15)

is a harmonic perturbation [23].

The terms in $\hat{\mathcal{H}}_{\rm{p}}$ that exhibit an explicit dependence
on $\xi$ can be removed through a near-identity canonical transformation of
variables. Adopting the Lie-series approach to perturbation theory [25, 26] we
take this transformation to emanate from a generating Hamiltonian, $\chi$ (in
this case, the un-averaged and averaged variables are related to one another
through $\Phi=\mathcal{S}_{\chi}^{t}\bar{\Phi}$,
$\phi=\mathcal{S}_{\chi}^{t}\bar{\phi}$, where $\mathcal{S}_{\chi}^{t}$ is the
Lie series operator [24]). Accordingly, carrying out the averaging procedure
to first order, the homologic equation reads:

$\displaystyle\hat{\mathcal{H}}_{\rm{p}}+\\{\check{\mathcal{H}}_{\rm{p}},\chi\\}=\frac{b\,\Delta}{2\,((1+b))}\,\Phi\,\cos(2\,\phi),$
(16)

where $\\{,\\}$ denote the Poisson bracket. For the problem at hand, equation
(16) is satisfied by:

$\displaystyle\chi$
$\displaystyle=-\frac{b\,(\Phi-4)\,\Phi}{16+b\,(3\,\Phi+8)}\sin(2\,\phi+\omega\,\xi)$
$\displaystyle+\frac{b\,\Phi^{2}}{8\,(16+b\,(3\,\Phi+8))}\sin(4\,\phi+2\,\omega\,\xi)$
$\displaystyle+\frac{4\,b\,\Delta\,\Phi}{b\,(3\,\Phi+8\,\omega+8)+8\,(\omega+2)}\sin(2(\phi+\omega\,\xi))$
$\displaystyle-\frac{b\,\Delta\,\Phi}{(b+1)\,\omega}\sin(\omega\,\xi).$ (17)

After this reduction, particle dynamics are governed by the integrable
Hamiltonian

$\displaystyle\bar{\mathcal{H}}_{\rm{p}}$
$\displaystyle=\bigg{(}\frac{1}{2}\bigg{(}\frac{2+b}{1+b}\bigg{)}-\frac{\omega}{2}\bigg{)}\,\Phi+\frac{3\,b}{32\,(1+b)}\,\Phi^{2}$
$\displaystyle+\frac{b\,\Delta}{2\,((1+b))}\,\Phi\,\cos(2\,\phi).$ (18)

where we have dropped the now-constant action $\Xi$, and the averaged nature
of the coordinates ($\Phi,\phi$) is implied. This expression agrees with
equation (18) of ref. 16, although our derivation is different.

The functional form of $\chi$ illuminates an important feature of the resonant
structure of the dynamics encapsulated by Hamiltonian (7) – namely that
conditions for $\chi$ to become singular are unphysical. In other words, the
system at hand appears to naturally circumvent the small divisor problem, at
least to leading order (see Ch. 2.3 of Morbidelli [24]). Qualitatively, this
implies that as long as Hamiltonian (7) constitutes a suitable approximation
to the dynamical evolution jointly facilitated by Hamiltonians (1) and (6),
the process of reducing $\mathcal{H}_{\rm{p}}$ to an integrable form through
time-averaging is justified under all physically relevant circumstances. In
turn, this means that the emergence of broad chaotic layer is not expected
within the parameter range of interest.

This interpretation is consistent with the numerical results in refs. 6, 7
(see also Figure 4 below). Nevertheless, we remark that the averaging
procedure – when not truncated at leading order (as done in equation 16) –
generates additional, higher-order resonant harmonics through ancillary terms
in the averaged Hamiltonian of the form
$\\{\\{\check{\mathcal{H}}_{\rm{p}},\chi\\},\chi\\}$,
$\\{\\{\\{\check{\mathcal{H}}_{\rm{p}},\chi\\},\chi\\},\chi\\}$, etc [23, 24].
If the strength of such terms increases with $b$, chaotic motion can be driven
by the modulation of the principal resonance’s separatrix by such high-order
resonances. We hypothesize that this is the reason why [7] observe chaotic
motion in numerical simulations with $b=100$ but not for lower values of the
space-charge parameter.

####  II.2.2 Particle-Core Interaction as an Andoyer Model

Beyond providing a purely analytical explanation for numerically computed
results, averaged Hamiltonian (18) is reminiscent of a mathematical pendulum.
Unlike a conventional pendulum, however, the harmonic term in expression (18)
possesses the D’Almbert characteristic of being multiplied by the action. In
fact, equation (18) belongs to a family of so-called Andoyer Hamiltonians
[27], which play a central role in the perturbative treatment of a broad array
of astrophysical phenomena, including the evolution of stellar and planetary
obliquities [28, 29, 30, 31], secular interactions [32, 33], as well as
orbital resonances [18, 34].

In general, Andoyer Hamiltonians have the form [35]:

$\displaystyle\mathcal{K}$
$\displaystyle=(2\,k-5)(1+k\,\delta)\,\Gamma+k\,\Gamma^{2}$
$\displaystyle+(3-k)^{1+k/2}\,\Gamma^{k/2}\,\cos(k\,\gamma),$ (19)

where $k$ is an integer greater than zero, $(\gamma,\Gamma)$ are action-angle
variables, and $\delta$ is the resonance proximity parameter – a number that
quantifies how close the system is to exact commensurability. Importantly, $k$
dictates the possible number of stable resonant libration regions in phase-
space, and for the problem at hand, $k=2$. Moreover, the topological character
of the phase-space portrait of $\mathcal{K}$ is regulated exclusively by
$\delta$, and it is straightforward to show that depending on the value of
$\delta$, the qualitative picture lies in one of three regimes: (i) for
$\delta<-1$, $\mathcal{K}$ only possesses a single elliptic equilibrium point
at the origin, and the dynamics are essentially equivalent to that of a
harmonic oscillator; (ii) for $-1\leqslant\delta<0$, the equilibrium point at
$\Gamma=0$ is rendered hyperbolic, and a figure-8 resonant separatrix emerges;
(iii) for $\delta\geqslant 0$, an inner circulation region develops, and the
separatrix takes on the shape of vesica piscis. To illustrate this, we show
the phase-space portraits of Hamiltonian (19) for $\delta=-3/2,-1/2,1/2,3/2$
in Figure (2).

Let us rewrite the particle-core interaction model in the form of equation
(19), with an eye towards expressing $\delta$ in terms of $b$ and $\Delta$.
The procedure of converting equation (18) into expression (19) consists fo two
consecutive steps: scaling of the actions, followed by a change of unit of
time. First, we take $\Gamma=\Phi/\eta$ and
$\mathcal{K}^{\prime}=\bar{\mathcal{H}}_{\rm{p}}/\eta$. Upon direct
substitution, we require the coefficient upfront the $\Gamma^{2}$ term to
differ from the harmoinc term by a factor of $2$, as dictated by equation (19)
for $k=2$. This immediately yields an equation for $\eta$:

$\displaystyle\frac{3\,b\,\eta}{32\,(1+b)}=\frac{b\,\Delta}{1+b},$ (20)

which reduces to $\eta=32\,\Delta/3$.

Figure 3: Resonance proximity parameter, $\delta$, as a function of the space
charge metric, $b$. The four depicted curves correspond to core oscillation
amplitudes of $\Delta=0.1$ (Blue), $\Delta=0.2$ (Red), $\Delta=0.3$ (Green),
and $\Delta=0.4$ (Gray). Within the parameter range of interest, $\delta$ is
always positive, implying that qualitatively, the phase-space portrait of the
uniform particle-core interaction model will resemble the two RHS panels of
Figure (2). Figure 4: Stroboscopic surfaces of section of particle motion,
with a uniformly charged core. The two LHS panels of the Figure have $b=1/3$
and lie in the emittance dominated regime, while the panels on the RHS have
$b=3$ and reside in the space charge dominated regime. Trajectories that
occupy the inner and outer circulation regions of the phase-space diagram are
shown in blue and gray respectively. Resonant orbits are shown in red. In all
cases, the qualitative character of the surface of section is identical to
that of the Andoyer model with $\Delta>0$. Moreover, even for large core
pulsation amplitude ($\Delta=0.4$), particle motion remains quasi-periodic, as
indicated by the fact that all trajectories fall on invariant tori.

Finally, we set $\mathcal{K}=\mathcal{K}^{\prime}\,(2(1+b))/(b\,\Delta)$,
which reduces the Hamiltonian to the form given by equation (19), with the
resonance proximity parameter given by:

$\displaystyle\delta$
$\displaystyle=\frac{b^{2}\left(\Delta^{2}-\Delta+\omega-1\right)+b\,(3\,\omega-2\,\Delta-4)+2\,(\omega-2)}{2\,b\,(b+2)\,\Delta}$
$\displaystyle\approx\frac{(\omega-2)+b\,(\omega-1)}{2\,b\,\Delta}-\frac{1}{2}+\mathcal{O}(\Delta).$
(21)

It is illuminating to examine the limiting regimes of equation (21). In the
space-charge dominated case of $b\gg 1$, the resonance proximity parameter
$\delta\approx(\sqrt{2}-1+\Delta^{2}-\Delta)/(2\,\Delta)$ is positive
definite. In the emittance dominated regime of $b\ll 1$,
$\delta\approx(1-2\,\Delta)/(4\,\Delta)$, meaning that the existence of a
stable inner circulation region of phase-space occupied by the core – granted
by the condition that $\delta>0$ – is also guaranteed as long as $\Delta<1/2$.

Although the possibility that $\delta$ can become negative for sufficiently
large values of the beam pulsation amplitude – causing the $\Phi=0$ fixed to
become hyperbolically unstable – is intriguing, we remind the reader that our
entire analysis is predicated on the assumption that $\Delta$ is small. To
this end, simultaneous numerical integration of equations of motion stemming
from Hamiltonians (1) and (6) demonstrates that for low values of $b$, the
$(u,p_{u})=(0,0)$ initial condition remains stable, even for large values of
$\Delta$. This suggests that the potential for $\delta$ to become smaller than
zero insinuated by equation (21) is an artifact of the $\Delta\ll 1$
assumption.

Figure (3) shows the resonance proximity parameter for $\Delta=0.1-0.4$ as a
function of $b$. Notably within our adopted range of beam parameters, $\delta$
is bracketed by $0\leqslant\delta\leqslant 2$. Importantly, the fact that
$\delta$ is of order unity despite large variations in underlying physical
parameters of the beam, signifies a pronounced insensitivity of the phase-
space portrait to the assumed values of the core pulsation amplitude,
$\Delta$, and the beam space-charge parameter, $b$. In other words, equation
(21) suggests that under all physically plausible conditions, the phase-space
portrait of the particle-core interaction model will always resemble the two
right-hand-side panels of Figure (2). To further illustrate this point, in
Figure (4) we show a series of stroboscopic surfaces of section, generated
from direct integration of equations of motion stemming from Hamiltonian (6).
Indeed, the qualitative agreement between the integrable approximation and
numerical experiments is satisfactory.

##  III A Gaussian Beam

Although the results presented in the proceeding section may illuminate the
underlying structure of the phase-space portrait of a particle embedded within
a uniformly charged beam, they do not yield any qualitative physical insights
that were not already evident within the established understanding of the
standard particle-core interaction model. That is, the notion of the existence
of resonant separatrixes that can transport particles from the outer edges of
the beam core to greater radial separations was already evident in the
pioneering work of Gluckstern [9] and was further demonstrated from numerical
grounds by Lagniel [11] and Wangler et al. [7]. In fact, the only genuinely
novel result of the analytic calculations presented in the previous section is
an explanation for why the basic aspects of the particle’s phase-space
portrait are persistently insensitive to the assumed values of beam
parameters.

Actual beams, however, are not uniformly charged, as assumed in the preceding
section. Rather, a broad collection of experimental measurements (see ref. 36
and the references therein) indicate that within real high-intensity linacs,
beam cores take on a profile that is closely approximated by a Gaussian
distribution. One example of such an observation of a Gaussian profile of a
space-charge-dominated beam stems from the SNS linac [37], where tune
depression achieves the value of $0.6$ [38], corresponding to space charge
parameter $b\approx 1.8$. Another example is the observation of a starkly non-
uniform beam profile in the J-PARC linac [39, 40], where tune depression
reaches the value of $0.5$ [41] corresponding to space charge parameter $b=3$.

It is important to understand that in real machines, the shape of the beam
profile is determined by a multitude of competing effects. Specifically,
space-charge forces themselves act to flatten the beam distribution [42],
while other phenomena, including misalignments of accelerating and focusing
elements, transverse-longitudinal coupling in the RF field, irregularities and
instabilities of RF and focusing fields, as well as excitation of higher-order
modes in the accelerating structures, act to widen the distribution.
Cumulatively, experimental results indicate that these phenomena result in
diffusive, Gaussian-type profile. Correspondingly, a Gaussian variant of the
particle-core interaction model represents a more self-consistent description
of particle dynamics.

As we show below, the more realistic case of a Gaussian beam exhibits greater
qualitative complexity. Particularly, within the context of a Gaussian beam,
the onset of large-scale stochastic motion can occur for physically plausible
parameters, and the resulting phase-space mixing can even connect the radial
excursions along the separatrixes to a larger section of the beam core.
Deriving specific conditions for the appearance of chaos and the delineation
of the associated irregular dynamics is the primary goal of this section.
However, in order to set the stage for our analysis of chaotic dynamics, it is
instructive to first derive an integrable model for a Gaussian beam, along the
same lines as what was done above.

###  III.1 A Semi-Analytic Model

The electromagnetic field of a beam with a Gaussian charge distribution has
the form:

$\displaystyle\Lambda=\frac{1}{u}\,\bigg{(}\frac{b}{1+b}\bigg{)}\,\bigg{(}1-\exp\bigg{(}-2\,\frac{u^{2}}{r^{2}}\bigg{)}\bigg{)}.$
(22)

We assume that $r$ still obeys equation (5) but is now interpreted as the
2-rms radius of the core. In doing so, we are assuming that beam oscillations
are dominated by space-charge forces, and are ignoring all auxiliary effects
that contribute to shaping the Gaussian profile of the beam. Nevertheless,
this interpretation is consistent with experimental data reported by Plum
[36], where the beam radius is shown to pulsate, while the beam profile
remains consistently Gaussian. Unlike the field of a uniform beam, the
function (22) smoothly varies from $\Lambda\propto u$ for $u\ll r$ to
$\Lambda\propto 1/u$ for $u\gg r$. Nevertheless, a power-series expansion
$\Lambda\approx(b/(1+b))(2\,u/r^{2}-2\,u^{3}/r^{4}+...)$ does not yield an
adequate approximation for equation (22) for $u\gtrsim r$, even if the degree
of the polynomial is high (Figure 5). Consequently, for the sake of
quantitative accuracy, here we proceed with the averaging method maintaining
the expression for $\Lambda$ in closed form.

Figure 5: The field of a Gaussian beam, $\Lambda$, evaluated setting the core
radius to its equilibrium value of $r=1$, is shown as a black curve.
Polynomial approximations to $\Lambda$ of various degrees are shown with
multi-colored curves. Clearly, power-series expansions fail to capture the
form of $\Lambda$ for $u\gtrsim 1$, even at a crude level. Figure 6:
Numerically averaged phase-space portraits of $\bar{\mathcal{H}}_{\rm{p}}$,
depicting the resonance in $(\Phi,\phi)$, given a Gaussian charge distribution
of the core. As in Figure (4), the two LHS and RHS panels lie in the emittance
dominated and space charge dominated regimes, respectively. Qualitatively, the
outlines of this resonance are indistinguishable from the Andoyer model with
$\delta>0$.

The Hamiltonian corresponding to particle motion within a Gaussian beam is:

$\displaystyle\mathcal{H}_{\rm{p}}$
$\displaystyle=\frac{p_{u}^{2}}{2}+\frac{u^{2}}{2}+\mathcal{T}$
$\displaystyle+\frac{b}{1+b}\bigg{(}\frac{\mathrm{Ei}(-2\,u^{2}/r^{2})}{2}-\ln(u)\bigg{)},$
(23)

where $\mathcal{T}$ is a dummy action conjugate to time and $\mathrm{Ei}$ is
the exponential integral. Changing variables in accord with canonical
transformations outlined in equations (9) and (13), the above expression
becomes:

$\displaystyle\mathcal{H}_{\rm{p}}$
$\displaystyle=\bigg{(}1-\frac{\omega}{2}\bigg{)}\,\Phi+\Xi-\frac{b\,\ln\big{(}\sqrt{2\Phi}\sin(\phi+\omega\,\xi/2)\big{)}}{1+b}$
$\displaystyle+\frac{b}{2\,(1+b)}\,\mathrm{Ei}\bigg{(}-\frac{4\,\Phi\,\sin^{2}(\phi+\omega\,\xi/2)}{r^{2}}\bigg{)},$
(24)

and we remind the reader that $\xi=t$.

Rather than attempting to reduce $\mathcal{H}_{\rm{p}}$ to a one-degree-of-
freedom Hamiltonian through the Lie series approach (which entails solving the
homologic equation 16 as in section II.2.1), here we instead take advantage of
separation of timescale inherent to the ($\Phi,\phi$) and ($\Xi,\xi$) degrees
of freedom. In particular, because $\dot{\phi}\sim 1-\omega/2\ll\dot{\xi}=1$,
the evolution of ($\Phi,\phi$) is adiabatic with respect to radial pulsation
of the core [43, 44]. In this regime, we are justified in directly averaging
the Hamiltonian over a single oscillation cycle in $\xi$:

$\displaystyle\bar{\mathcal{H}}_{\rm{p}}$
$\displaystyle=\frac{\omega}{2\,\pi}\oint\big{(}\mathcal{H}_{\rm{p}}-\Xi\big{)}\,d\xi.$
(25)

Because this procedure reduces $\bar{\mathcal{H}}_{\rm{p}}$ to an integrable
system semi-analytically, in order to illuminate the associated dynamics, we
can simply project the level curves of $\bar{\mathcal{H}}_{\rm{p}}$ onto the
($\Phi,\phi$) plane. Figure (6) shows a sequence of phase-space portraits of
$\bar{\mathcal{H}}_{\rm{p}}$ for $b=1/3$ and $b=3$ (we have also computed
analogous Figures for $b=0.1,1,10$, and found them to be very similar to those
depicted in Figure 6). Crudely speaking, this exercise suggests that the
qualitative character of averaged dynamics within a Gaussian and a uniform
beam are identical: in both cases, an inner circulation region is encircled by
a vesica piscis-shaped separatrix, with two resonant equilibria residing on
the $u$-axis. Moreover, the pronounced insensitivity of the phase-space
portraits to the beam parameters $b$ and $\Delta$ appears to be in line with
the functional form for the resonance proximity parameter suggested by
equation (21).

A more cursory inspection of the Hamiltonian, however, reveals that even at a
basic level, the picture cannot be so simple, since a prominent additional
resonance – previously removed from $\mathcal{H}_{\rm{p}}$ via the averaging
process (25) – emerges as an important driver of particle dynamics for
sufficiently large values of $b$. To demonstrate this, let us return to
equation (23), set $\Delta=0$ (such that $r=1$, which renders
$\mathcal{H}_{\rm{p}}$ autonomous), and examine its equilibrium. In
particular, it is instructive to consider the Hessian matrix of
$\mathcal{H}_{\rm{p}}$, evaluated at the $(u,p_{u})=(0,0)$ fixed point. Since
there is no $u-p_{u}$ coupling in the Hamiltonian,
$\partial^{2}\mathcal{H}_{\rm{p}}/\partial u\,\partial p_{u}=0$. It is further
trivial to see that $\partial^{2}\mathcal{H}_{\rm{p}}/\partial p_{u}^{2}=1$.
The determinant of the Hessian matrix is then simply

$\displaystyle\mathcal{D}$
$\displaystyle=\frac{\partial^{2}\mathcal{H}_{\rm{p}}}{\partial
u^{2}}=\lim_{u\rightarrow
0}\bigg{[}1+\frac{b}{1+b}\frac{e^{-2\,u^{2}}(e^{2\,u^{2}}-4\,u^{2}-1))}{u^{2}}\bigg{]}$
$\displaystyle=1-\frac{2\,b}{1+b}<0\ \ \ \ \ \mathrm{if}\ b>1.$ (26)

This rudimentary analysis shows that as $b$ crosses the threshold of unity
(which coincides with the transition from the emittance dominated regime to
the space-charge dominated regime), the origin of the $(u,p_{u})$ phase-space
transitions form being an elliptic fixed point (corresponding to a minimum of
$\mathcal{H}_{\rm{p}}$) to a hyperbolic fixed point (corresponding to a saddle
point of $\mathcal{H}_{\rm{p}}$), signaling the appearance of a resonant
separatrix for $b>1$. Importantly, these dynamics arise completely
independently of the breathing mode of the core. The emergence of a figure-
eight separatrix, along with new elliptic equilibria at $|u|>0$, $p_{u}=0$ is
shown in Figure (7), where we plot the phase-space portrait of
$\mathcal{H}_{\rm{p}}$ for $b=1/3,1,3,10$, setting $\Delta=0$. An important
consequence of the appearance of this resonance is that when $\Delta$ is made
sufficiently large, this resonance can overlap the resonance shown in Figure
(6), driving prominent chaotic motion within the core.

We note that herein lies an important distinction between the Gaussian and
uniformly-charged variants of the particle-core interaction model. That is,
contrary to the dynamics encapsulated by the Hamiltonian (23), linear
stability of the origin in a $\Delta=0$ uniformly charged beam is trivially
ascertainable. In particular, from equation (6), it is easy to see that near
the $(u,p_{u})=(0,0)$ point, the Hamiltonian has the form of a harmonic
oscillator: $\mathcal{H}_{\rm{p}}=p_{u}^{2}/2+(1/(1+b))\,u^{2}/2$, implying a
determinant of the Hessian matrix that is positive definite. Consequently, a
nonlinear resonance of the kind depicted in Figure (7) can never arise within
the uniformly-charged model, meaning that chaotic motion within the core is
only expected to arise within the context of a Gaussian beam.

Numerically computed surfaces of section shown in Figure (8) largely confirm
this expectation. In particular, unlike the case of a uniformly-charged beam,
where the trajectories that comprise the core remain largely quasi-periodic,
here a conspicuous stochastic layer appears close the origin for $b\gtrsim 1$.
These stroboscopic plots can be compared to equivalent surfaces of section of
the uniformly charged beam depicted in Figure (4). Importantly, in both cases,
the surfaces of section were obtained through direct integration of the
equations of motion derived from Hamiltonians (1, 6, 23), and are not
influenced by any of the approximations (e.g., polynomial fit to the field,
etc) invoked above.

A crucial consequence of the appearance of a stochastic layer within the beam
core is that it allows for chaotic transport to ensue over a broader phase-
space domain. That is, for sufficiently large values of $\Delta$, this layer
connects to the principal resonance’s separatrix, effectively linking the beam
core and the halo through chaotic mixing. Let us now quantify the onset of
chaos within a Gaussian particle-core interaction model, and derive the
relevant criteria from analytic grounds.

Figure 7: Phase-space portraits of
$u=\sqrt{2\,\Psi}\,\sin(\psi),p_{u}=\sqrt{2\,\Psi}\,\cos(\psi)$ dynamics
within a Gaussian beam for $\Delta=0$. A nonlinear resonance appears for
$b>1$, rendering the origin of the phase-space diagram hyperbolically
unstable. Importantly, this resonance is distinct from that generated by the
breathing mode of the core (when $\Delta>0$), and does not emerge in the
uniformly charged variant of the particle-core interaction model. The phase-
space topology of this resonance is equivalent to that of the Andoyer model
with $-1<\delta<0$.

###  III.2 Onset of Chaos

A quantitative understanding for conditions under which dynamical systems
begin to exhibit stochastic behavior dates back to the pioneering work of
Chirikov [45, 46], and by now, the essential statement of the Chirikov
criterion is widely known: in order for global chaos to ensue, the sum of the
half-widths of a pair of resonances in action space must exceed the distance
between the resonances’ equilibria. One approach towards calculating this
criterion in the context of charged particle motion within a Gaussian beam, is
to compute the locations and widths of individual resonant separatrixes semi-
analytically, as shown in Figures (6) and (7). A considerably less accurate
(see Figure 5), but more qualitatively informative approach is to adopt a
series-expansion of the governing Hamiltonian, with an eye towards reducing
$\mathcal{H}_{\rm{p}}$ to the simplest possible model that still captures
irregular dynamics, and estimating the Chirikov criterion in a purely
analytical manner. This is the strategy we adopt in this section.

As a starting point, let us write down an approximate variant of Hamiltonian
(23), truncating the polynomial expansion of $\Lambda$ at cubic order:

$\displaystyle\mathcal{H}_{\rm{p}}$
$\displaystyle=\frac{p_{u}^{2}}{2}+\frac{u^{2}}{2}-\frac{b}{1+b}\bigg{(}\frac{u^{2}}{r^{2}}-\frac{u^{4}}{2\,r^{4}}\bigg{)}.$
(27)

This Hamiltonian is clearly reminiscent of equation (7), although the quartic
term in the above expression is notably stronger, and is dependent on $r$. We
now follow a procedure similar to that outlined in section II.2.1. We
substitute equation (5) for $r$, but to keep the model as compact as possible,
we expand the resulting expression to leading order in $\Delta$. Moreover, in
spirit of a series-expansion approach, upon transforming to action-angle
variables ($\Psi,\psi$) via equation (9), we neglect all harmonics
proportional to $\Psi^{2}$, retaining only those that are linear in $\Psi$.
This reduces the Hamiltonian to the following simplified form:

$\displaystyle\mathcal{H}_{\rm{p}}$
$\displaystyle\approx\frac{\Psi}{1+b}+\frac{3\,b\,\Psi^{2}}{4\,(1+b)}$
$\displaystyle+\frac{b\,\Psi}{1+b}\big{(}1+2\,\Delta\,\cos(\omega\,t)\big{)}\cos(2\,\psi).$
(28)

Importantly, Hamiltonian (28) has a similar mathematical structure to the
Hamiltonian of a modulated pendulum. As a mechanical system, a modulated
pendulum is nothing more than a mass suspended via a rigid rod, whose base
hangs on a spring that oscillates with frequency $\omega$. In turn,
oscillations of the base result in a periodic variation in the effective
acceleration due to gravity [23]. Although qualitatively simple, depending on
system parameters and initial conditions, the modulated pendulum can display a
broad range of behavior, including strictly periodic motion, quasi-periodic
resonant trajectories, motion that lies on KAM tori, as well as full-fledged
deterministic chaos. For this reason, the modulated pendulum is routinely
adopted as a paradigm for stochastic dynamics [47, 48, 49, 24].

The application of the Chirikov resonance overlap criterion to the modulated
pendulum is straightforward. If we momentarily ignore the dependence of the
cosine term in $\mathcal{H}_{\rm{p}}$ upon $\Psi$, it is easy to see how
sinusoidal variation its strength would give rise to a total of three
resonances. In particular, because
$(1+2\,\Delta\,\cos(\omega\,t))\,\cos(2\,\psi)=\cos(2\,\psi)+\Delta\,\cos(2(\psi-\omega\,t/2))+\Delta\,\cos(2(\psi+\omega\,t/2))$,
the fixed points of these resonances would be separated in action space
(equivalently, frequency space) by $\omega/2$. Accordingly, global chaos would
ensue when the half-widths of these resonant multiplets become comparable to
$\omega/2$ such that the separatrixes cross.

Figure 8: A gallery of stroboscopic surfaces of section of particle dynamics
within a Gaussian beam. The same color scheme as that employed in Figure (4)
is also adopted here, with the added element of chaotic trajectories being
depicted in purple. In the emittance dominated ($b=1/3$) and transitionary
($b=1$) regimes, dynamical evolution within a Gaussian beam is largely
equivalent to that of the uniformly charged beam, and follows the well-known
narrative of the particle-core interaction model. On the other hand, in the
space charge-dominated regime (corresponding to panels with $b=3$ and $b=10$),
a discernible chaotic layer emerges close to the origin, facilitating
unpredictable particle evolution. Importantly, for sufficiently large core
pulsation amplitudes, separatrixes of the “intrinsic” Gaussian resonance (in
$\psi$) and the “breathing” resonance (in $\phi$) overlap strongly, connecting
the core and the halo through rapid stochastic transport. RHS panels where
chaotic core-halo mixing ensues are labeled with a checkmark, whereas those
where the chaotic layer is confined to the core are marked with a diagnoal
cross.

All of the dynamical features of a conventional modulated pendulum also
naturally arise within the context of Hamiltonian (28). The analogy, however,
is inexact in one important respect: while a standard modulated pendulum
always encompasses three resonant harmonics with well-defined equilibrium
conditions, the system at hand possesses the D’Almbert characteristic, which
limits the physically meaningful domain of the action $\Psi$ to positive
values (such that $u$ and $p_{u}$ are real). Additionally, the D’Almbert
characteristic modifies the distance between the resonant equilibrium points.
This means that in order for chaos to ensue, not only do the resonances have
to overlap, but at least two of the corresponding equilibria have to reside at
$\Psi>0$. As we already demonstrated in the previous sub-section via equation
(26), this condition translates to a requirement that $b>1$. We remind the
reader that this constraint can be written as $\sigma/\sigma_{0}<1/\sqrt{2}$,
corresponding to the transition from emittance-dominated regime to the space-
charge dominated regime. Assuming that this constraint is satisfied, let us
now write down the approximate criterion for the onset of chaos.

Unlike the width of the unperturbed resonance in $\psi$ (which only exists for
$b>1$ but is largely independent of $\Delta$), the width of the resonance
corresponding to $\phi=\psi-\omega/2$ depends on $\Delta$. Consequently, a
reasonable way to formulate the Chirikov criterion is to equate the distance
between the stable equilibrium points of the $\psi$ and $\phi$ resonances to
the half-width of the resonance in $\phi$. Let us now compute the relevant
quantities. The equilibrium action of the unperturbed $(\Psi,\psi)$ dynamics,
$[\Psi]$, follows from Hamilton’s equations in the $\Delta\rightarrow 0$
limit:

$\displaystyle\frac{d\,\psi}{dt}=\frac{\partial\,\mathcal{H}_{\rm{p}}}{\partial\,\Psi}=\frac{1}{1+b}+\frac{3\,b\,[\Psi]}{2\,(1+b)}+\frac{b\,\cos(2\,[\psi])}{1+b}=0$
$\displaystyle\frac{d\,\Psi}{dt}=-\frac{\partial\,\mathcal{H}_{\rm{p}}}{\partial\,\psi}=-\frac{2\,b\,[\Psi]\,\sin(2\,[\psi])}{1+b}=0.$
(29)

From inspection, it is evident that the equilibrium value of the angle lies at
$[\psi]=\pi/2,3\pi/2$ such that $\cos(2[\psi])=-1$. The value of the
unperturbed equilibrium action is then

$\displaystyle[\Psi]=\frac{2\,(b-1)}{3\,b}.$ (30)

The functional form of $[\Psi]$ is yet another reminder of the fact that
chaotic dynamics can only appear for $b\geqslant 1$: for values of $b$ smaller
than unity, resonance overlap is precluded by $[\Psi]$ taking on a negative
value, such that that the resonant equilibrium resides on the imaginary
$(u,p_{u})$ plane.

The value of the equilibrium action $[\Phi]$ can be computed in a similar
manner. Upon transforming variables according to equations (13) and dropping
(i.e., “averaging away”) $\xi$-dependent terms, Hamilton’s equations give:

$\displaystyle\frac{d\,\phi}{dt}=\frac{\partial\,\mathcal{H}_{\rm{p}}}{\partial\,\Phi}=\frac{1-\Delta\,b}{1+b}-\frac{\omega}{2}+\frac{3\,b\,[\Phi]}{2\,(1+b)}=0,$
(31)

where we have set $\cos(2\,[\phi])=-1$ by analogy with equation (29). After
mild rearrangement, the distance between resonant equilibria in action-space
reads:

$\displaystyle[\Phi]-[\Psi]=\frac{\omega+b\,(\omega+2\,\Delta-2)}{3\,b}.$ (32)

The remaining step is to compute the resonance half-width $\delta\Phi$. In the
pendulum approximation – which is justifiable for the $\phi$-resonance as long
as $\delta\Phi\lesssim[\Phi]$ [50] – we can substitute $[\Phi]$ for the value
of the action upfront the cosine term in $\mathcal{H}_{\rm{p}}$. This gives
[24]:

$\displaystyle\delta\Phi\approx
2\sqrt{\bigg{(}\frac{b\,\Delta\,[\Phi]}{1+b}\bigg{)}\bigg{(}\frac{2\,(1+b)}{3\,b}\bigg{)}}=\sqrt{\frac{8\,\Delta\,[\Phi]}{3}}.$
(33)

Recalling that the existence of the resonance in the ($\psi,\Psi$) degree of
freedom requires that $b>1$, the Chirikov criterion
$([\Phi]-[\Psi])/\delta\Phi\lesssim 1$ can now be written as

$\displaystyle\frac{\omega+b\,(\omega+2\,\Delta-2)}{2\,\sqrt{2}\sqrt{b\,\Delta\,(\omega+b\,(\omega+2\,\Delta)-2)}}\lesssim
1\ \ \mathrm{and}\ \ b>1.$ (34)

Although it is possible to express $b$ as a function of $\Delta$ (or
alternatively $\Delta$ as a function of $b$) in closed form, the resulting
expressions are cumbersome, and do not yield any additional insight that is
not already evident in the above expression.

Figure (9) depicts the Chirikov criterion (34) as a thick cyan curve on the
$(\Delta,b)$ plane. Upon inspection, two physically important characteristics
of this threshold are immediately evident. First, for diminishing values of
$\Delta$, progressively higher values of $b$ are necessary to drive chaotic
motion. This relationship stems from the $\sim 1/\sqrt{\Delta}$ dependence in
equation (34), and simply insinuates that in the limit of vanishing $\Delta$,
$\mathcal{H}_{\rm{p}}$ reduces to an integrable form, implying perfectly
regular dynamics.

Second, for $\Delta>(\sqrt{3}-1)/3\approx 0.24$, the expression (34) suggests
that a stochastic layer can emerge for a value of $b$ below unity. This is,
however, an artifact of our pendulum-like formulation of the Chirikov
criterion, and as we have already discussed above, $b>1$ constitutes a
necessary but insufficient condition for resonance overlap. Stated
differently, our analysis suggests that for core pulsation amplitude in excess
of $\sim 25\%$, particle dynamics within a Guassian beam become chaotic
essentially as soon as the space-charge dominated regime is reached.

Before leaving this sub-section, let us briefly comment on the suitability of
the assumed profile of the beam core. Aforementioned experimental support for
a Gaussian beam profile aside, the emergence of a prominent nonlinear
resonance within a stationary beam – and bonafide stochasticity within a
pulsating one – begs the question of self-consistency. That is, if initially
stationary particles positioned in the center of the beam begin to execute
significant radial excursions as soon as the beam becomes space charge-
dominated, surely the particle-core interaction model implies that a Gaussian
distribution cannot be maintained indefinitely. In fact, this is precisely
what is observed in idealized multi-particle numerical simulations: on rather
short timescales, an initially Gaussian beam exhibits the propensity to self-
regulate, rapidly approaching a uniform charge distribution [51]. It is
important to understand however, that this tendency is suppressed in real
machines, underscoring the limitations of numerical experiments that only
include space charge as a means of beam profile modulation, and fail to
account for other important physical effects. In particular, lattice
misalignments and beam mis-steering, together with nonlineaities arising from
the RF field and focusing elements, act to restore the (approximately)
Gaussian nature of the beam, suggesting that while chaotic diffusion within
the particle core interaction model is driven by nonlinear space-charge
forces, the very nonlinearity of these effects are maintained through
extrinsic forcing.

###  III.3 Numerical Experiments

The analytic conclusions of the preceding discussion can be tested and
extended with the aid of numerical simulations – an exercise we now carry out.
The series of questions we seek to answer in the remainder of this section are
three-fold. First, we wish to quantify the validity of equation (34) within
the context of our system. Second, we aim to compute the characteristic
Lyapunov exponent pertinent to charged particle dynamics. Finally, we wish to
address the practically important issue of mapping out conditions under which
chaotic diffusion within the core ceases to be bounded, and links to more
extended radial excursions, facilitating an accelerated formation of a beam
halo. To address these questions, we performed the following suite of
numerical experiments.

On a grid of beam parameters spanning $b\in(0,10),\Delta\in(0.01,0.4)$, we
computed $100$ stroboscopic surfaces of section, akin to those shown in Figure
(8). To populate each surface of section, we chose 30 initial conditions lying
on the $u-$axis, spanning the range $u_{0}\in(0,1-\Delta)$. We note that this
choice is made strictly for definitiveness, and we have confirmed that our
results do not exhibit any significant dependence on the distribution of
initial conditions. Recalling that our surfaces of section correspond to a
phase of core pulsation when $r$ is minimized, these initial conditions
correspond to ions that originate inside the core with zero radial velocity.
Every trajectory was evolved for 1000 core pulsation periods.

In conjunction with each integration of the equations of motion, we
additionally evolved a second, shadow trajectory (initialized in close
proximity to the main trajectory at
$\acute{p_{u}}=\epsilon,\acute{u}=u_{0}+\epsilon$, where $\epsilon=10^{-6}$)
through integration of linearized equations of motion (see appendix B).
Keeping track of the distance between the primary and shadow trajectories in
phase-space, $\nu=\sqrt{(p_{u}-\acute{p_{u}})^{2}+(u-\acute{u})^{2}}$, we
compute the Lyapunov coefficient, $\lambda$, of the trajectory. To be more
specific, we computed the so-called mean exponential growth factor of nearby
orbits [52]:

$\displaystyle
Y=\frac{2}{t}\int_{0}^{t}\frac{\dot{\nu}}{\nu}t^{\prime}\,dt^{\prime},$ (35)

which known to converge more rapidly than the standard renormalization method
for computing $\lambda$ [53]. To this end, we note that the running time-
average of $Y$ tends towards $\bar{Y}\sim 2$ if the trajectory is regular, and
towards $\bar{Y}\sim\lambda\,t/2$ if the trajectory is chaotic.

If all initial conditions on the surface of section yield quasi-period
dynamics, we record $\lambda=0$ for the given combination of $\Delta$ and $b$.
If, on the other hand, some trajectories result in chaotic motion, we record
the mean value of $\lambda$, obtained by averaging the Lyapunov coefficients
of chaotic initial conditions (i.e., initial conditions that do not give
stochastic behavior are excluded). The resulting distribution of $\lambda$ is
shown as a background contour-plot on Figure (9). Although imperfect, the
analytic criterion for the onset of chaos (equation 34) clearly tracks the
boundary between regular and stochastic dynamics. Thus, in spite of a series
of relatively crude approximations that were made in reducing the Gaussian
particle core interaction model to a simplified modulated pendulum-like
Hamiltonian (28), the resulting criterion for the onset of chaos displays
satisfactory agreement with numerical simulations.

Figure 9: Regular and chaotic dynamics within the Gaussian particle-core
interaction model. The background contour plot informs the topological
character of the stroboscopic surface of section by reporting the Lyapunov
coefficient, $\lambda$, where chaotic motion is detected. The analytical
Chirikov criterion for the onset of stochasticity (equation 34) is marked with
a cyan curve. Note that while the Chirikov criterion yields $b\approx 2.3$ for
the lowest numerically tested beam pulsation amplitude of $\Delta=0.01$, this
curve turns sharply towards $b\rightarrow\infty$ as $\Delta\rightarrow 0$,
where the dynamics becomes integrable. Finally, the threshold for chaotic
mixing between the core and the halo is outlined with a black curve.

Beyond the binary question of whether or not a given combination of beam
parameters yields chaotic motion or not, the value of $\lambda$ itself tells
an intriguing story. That is, the distribution of $\lambda$ is clearly
bimodal, with moderate values of $\Delta$ and $b$ corresponding to
$\lambda\sim 0.1-0.3$ and large values of $\Delta$ and $b$ corresponding to
$\lambda\gtrsim 0.5$. Crucially, the boundary between these two domains
appears rather sharp, and corresponds to the onset of chaotic linkage of the
separatrixes of $\phi$ and $\psi$ resonances. To check this, we kept track of
whether or not trajectories that originate within the core attain a value of
$u$ that is greater than that of the stable equilibrium point of the
$\phi$-resonance (computed using semi-analytical averaging as in section
III.1) on the surface of section, and derived an approximate threshold for
such mixing. This threshold is shown on Figure (9) as a black curve, and
crudely tracks the boundary between slow and rapid diffusion. Thus, the region
of parameter space characterized by the larger value of the Lyapunov
coefficient also corresponds to the regime where halo formation is facilitated
by enhanced stochastic mixing of trajectories, which joins the beam core and
halo through rapid chaotic transport, along the outlines of resonant
separatrixes in phase-space.

##  IV Summary

In this work, we have carried out a perturbative, as well as numerical
investigation of the particle-core interaction model of halo formation, with
an eye towards understanding the role played by large-scale chaotic diffusion,
and have quantified the conditions for its onset. A key result of our analysis
is the demonstration that the emergence of globally stochastic dynamics within
the framework of the particle-core interaction paradigm hinges upon the
assumed charge distribution within the core. That is, while particle motion
within a uniform beam is largely regular, an appreciable fraction of phase-
space occupied by the core can become chaotic given physically plausible
parameters if the charge distribution of the beam is Gaussian.

The analytic calculations presented in sections II.2 and III.1 explain this
difference from qualitative grounds. In particular, for the case of a uniform
beam, we show that the particle-core interaction model can be reduced to an
integrable Andoyer-type Hamiltonian through canonical perturbation theory, and
that the associated generating Hamiltonian is strictly non-singular, meaning
that the averaging procedure is justified – at least to first order in the
resulting series [24] – for all physically plausible beam parameters. This
elucidates why chaotic motion is confined to a thin band around the separatrix
when the charge-distribution of the core is uniform [7, 16].

The case of a Gaussian beam is considerably more complex. While a uniform
beam-like qualitative picture holds in the emittance-dominated regime, our
analysis reveals that as soon as the system crosses over into the space-charge
dominated regime, a second principal resonance emerges in phase space in
addition to the usual breathing mode of the particle-core interaction model.
This resonance renders the center of the beam hyperbolically unstable.
Modulation of the corresponding figure-$\infty$ separatrix by the core’s
breathing mode gives rise to chaotic dynamics within the core. Moreover, for
sufficiently high values of the beam’s space charge parameter $b$ and beam
pulsation amplitude $\Delta$, the two principal resonances can overlap
strongly, leading to enhanced chaotic mixing between the core and the halo
(Figure 9). Indeed, the boundaries of model parameter space within which
deterministic stochasticity reigns are predictable, and are given by a simple
analytic expression (34). Accordingly, the analysis presented herein outlines
a theoretical framework within the context of which more detailed numerical
and experimental results can be interpreted.

###### Acknowledgements.

K.B. is grateful to Caltech as well as the David and Lucile Packard Foundation
for their generous support. Y.B. is grateful to the Los Alamos National
Laboratory and the U.S. Department of Energy for their generous support, under
US DOE contract 89233218CNA000001.

## Data Availability Statement

The data that support the findings of this study are available from the
corresponding author upon reasonable request.

##  Appendix A An approximate solution for the pulsation of the core

In the vicinity of the $(\rho,p_{\rho})=(0,0)$ equilibrium, Hamilton’s
equations of motion – applied to expression (1) – give:

$\displaystyle\frac{d\,\rho}{d\,t}=\frac{\partial\,\mathcal{H}_{\rm{c}}}{\partial\,p_{\rho}}=p_{\rho}$
$\displaystyle\frac{d\,p_{\rho}}{d\,t}=-\frac{\partial\,\mathcal{H}_{\rm{c}}}{\partial\,\rho}=-\omega^{2}\,\rho+\zeta^{2}\rho^{2}.$
(36)

Following standard practice [22, 23], we note that for low-amplitude
oscillations the linear term in $\rho$ will dominate (yielding a harmonic
oscillator), and consider the trial solution

$\displaystyle\rho\approx\Delta\,\cos(\omega\,t+\varphi)+\zeta^{2}\,f,$ (37)

where $f$ is a yet-to-be-determined function of time.

Writing equations (36) as a single second order differential equation for
$\rho$, and substituting the trial solution, we have:

$\displaystyle\zeta^{2}\,\frac{d^{2}\,f}{d\,t^{2}}-\Delta\,\omega^{2}\,\cos(\omega\,t+\varphi)=-\Delta\,\omega^{2}\,\cos(\omega\,t+\varphi)$
$\displaystyle+\Delta^{2}\,\zeta^{2}\,\cos^{2}(\omega\,t+\varphi)-\zeta^{2}\,\omega^{2}\,f+\zeta^{6}\,f^{2}$
$\displaystyle+2\,\Delta\,\zeta^{4}\cos(\omega\,t+\varpi).$ (38)

To second order in $\zeta$, the above equation of motion simplifies to:

$\displaystyle\frac{d^{2}\,f}{d\,t^{2}}=\Delta^{2}\,\cos^{2}(\omega\,t+\varphi)-\omega^{2}\,f.$
(39)

This expression is satisfied by the solution:

$\displaystyle
f=\frac{\Delta^{2}\,\zeta^{2}}{2\,\omega^{2}}-\frac{\Delta^{2}\,\zeta^{2}}{6\,\omega^{2}}\,\cos(2\,(\omega\,t+\varphi)).$
(40)

Equation (5) is thus obtained. Although it is straightforward to extend this
procedure to derive higher-order corrections to this solution, Figure (1)
indicates that for our purposes a second-order approximation is sufficient.

##  Appendix B Linearized Equations of Motion

A well-known hallmark of chaotic evolution is the sensitive dependence on
initial conditions, which manifests through exponential divergence of nearby
trajectories. Since this divergence only emerges locally, it is standard
practice to linearize the equations of motion and compute the phase-space
distance between the primary and shadow trajectories $\delta u=\acute{u}-u$,
$\delta p_{u}=\acute{p_{u}}-p_{u}$ according to:

$\displaystyle\frac{d\,\delta
u}{dt}=\frac{\partial^{2}\,\mathcal{H}_{\rm{p}}}{\partial p_{u}^{2}}\,\delta
p_{u}+\frac{\partial^{2}\,\mathcal{H}_{\rm{p}}}{\partial u\,\partial
p_{u}}\,\delta u$ $\displaystyle\frac{d\,\delta
p_{u}}{dt}=-\frac{\partial^{2}\,\mathcal{H}_{\rm{p}}}{\partial u\,\partial
p_{u}}\,\delta p_{u}-\frac{\partial^{2}\,\mathcal{H}_{\rm{p}}}{\partial
u^{2}}\,\delta u$ (41)

For the problem at hand, these equations of motion take on a particularly
simple form. As we already saw in section III.1, the mixed derivatives vanish
because there is no $u-p_{u}$ coupling in the system. Moreover, the second
derivative of $\mathcal{H}_{\rm{p}}$ with respect to $p_{u}$ is just unity,
and $\partial^{2}\mathcal{H}_{\rm{p}}/\partial u^{2}$ is given by equation
(26). As a result, we have:

$\displaystyle\frac{d\,\delta
p_{u}}{dt}=-\bigg{(}1+\frac{b}{1+b}\frac{e^{-2\,u^{2}}(e^{2\,u^{2}}-4\,u^{2}-1))}{u^{2}}\bigg{)}\delta
u$ $\displaystyle\frac{d\,\delta u}{dt}=\delta p_{u}$ (42)

## References

  * Wangler [1998] Wangler, T. 1998. RF Linear Accelerators. RF Linear Accelerators, by Thomas Wangler, pp. 397. ISBN 0-471-16814-9. Wiley-VCH , May 1998. 397.
  * Batygin [1996] Batygin, Y. K. 1996. Conservation of space-charge-dominated beam emittance in a strong nonlinear focusing field. Physical Review E 53, 5358–5364. doi:10.1103/PhysRevE.53.5358
  * Batygin [1998] Batygin, Y. K. 1998. Self-consistent distribution of a high brightness beam in a continuous focusing channel and application to halo-free beam transport. Physical Review E 57, 6020–6029. doi:10.1103/PhysRevE.57.6020
  * Sonnad and Cary [2005] Sonnad, K. G., Cary, J. R. 2005. Control of beam halo formation through nonlinear damping and collimation. Physical Review Accelerators and Beams 8. doi:10.1103/PhysRevSTAB.8.064202
  * Pathak and Krishnagopal [2017] Pathak, A., Krishnagopal, S. 2017. Higher order mode beams mitigate halos in high intensity proton linacs. Physical Review Accelerators and Beams 20. doi:10.1103/PhysRevAccelBeams.20.014201
  * Jameson [1993] Jameson, R.A. 1993. Beam-Halo From Collective Core/Single-Particle Interactions. Los Alamos Report No. LA-UR-93-1209
  * Wangler et al. [1998] Wangler, T. P., Crandall, K. R., Ryne, R., Wang, T. S. 1998. Particle-core model for transverse dynamics of beam halo. Physical Review Accelerators and Beams 1.
  * Qiang et al. [2004] Qiang, J., Ryne, R. D., Hofmann, I. 2004. Space-Charge Driven Emittance Growth in a 3D Mismatched Anisotropic Beam. Physical Review Letters 92.
  * Gluckstern [1994] Gluckstern, R. L. 1994. Analytic model for halo formation in high current ion linacs. Physical Review Letters 73, 1247–1250.
  * Lagniel [1994] Lagniel, J.-M. 1994. On halo formation from space-charge dominated beams. Nuclear Instruments and Methods in Physics Research A 345, 46–53.
  * Lagniel [1994] Lagniel, J.-M. 1994. Chaotic behaviour and halo formation from 2D space-charge dominated beams. Nuclear Instruments and Methods in Physics Research A 345, 405–410.
  * Allen and Wangler [2001] Allen, C., Wangler, T. 2001. Parameters for Quantifying Beam Halo. 19th Particle Accelerator Conference (PAC 2001), 1732–1734.
  * Gerigk [2004] Gerigk, F. 2004. Beam halo in high-intensity hadron accelerators caused by statistical gradient errors. Physical Review Accelerators and Beams 7.
  * Lund et al. [2007] Lund, S. M., Barnard, J. J., Bukh, B., Chawla, S. R., Chilton, S. H. 2007. A core-particle model for periodically focused ion beams with intense space-charge. Nuclear Instruments and Methods in Physics Research A 577, 173–185.
  * Ikegami [1999] Ikegami, M. 1999. Particle-core analysis of beam halo formation in anisotropic beams. Nuclear Instruments and Methods in Physics Research A 435, 284–296.
  * Batygin [2010] Batygin, Y. K. 2010, Analytical treatment of particle-core interaction, Nuclear Instruments and Methods in Physics Research A, 618, 37
  * Allen et al. [2002] Allen, C. K. and 12 colleagues 2002. Beam-Halo Measurements in High-Current Proton Beams. Physical Review Letters 89.
  * Henrard and Lemaitre [1983] Henrard, J., Lemaitre, A. 1983. A Second Fundamental Model for Resonance. Celestial Mechanics 30, 197–218.
  * Qiang and Ryne [2000] Qiang, J., Ryne, R. D. 2000. Beam halo studies using a three-dimensional particle-core model. Physical Review Accelerators and Beams 3. doi:10.1103/PhysRevSTAB.3.064201
  * Goldstein et al. [2002] Goldstein, H., Poole, C., Safko, J. 2002. Classical mechanics. Classical mechanics (3rd ed.) by H. Goldstein, C. Poolo, and J. Safko. San Francisco: Addison-Wesley, 2002..
  * Dixon et al. [1991] Dixon, J. M., Tuszyński, J. A., Otwinowski, M. 1991. Special analytical solutions of the damped-anharmonic-oscillator equation. Physical Review A 44, 3484–3491.
  * Landau and Lifshitz [1969] Landau, L. D., Lifshitz, E. M. 1969. Mechanics. Course of Theoretical Physics, Oxford: Pergamon Press, 1969, 2nd ed..
  * Lichtenberg and Lieberman [1983] Lichtenberg, A. J., Lieberman, M. A. 1983. Regular and stochastic motion. Applied Mathematical Sciences, New York: Springer, 1983.
  * Morbidelli [2002] Morbidelli, A. 2002. Modern celestial mechanics : aspects of solar system dynamics. London: Taylor & Francis, 2002
  * Hori [1966] Hori, G. 1966. Theory of General Perturbation with Unspecified Canonical Variable. Publications of the Astronomical Society of Japan 18, 287.
  * Deprit [1969] Deprit, A. 1969. Canonical transformations depending on a small parameter. Celestial Mechanics 1, 12–30.
  * Ferraz-Mello [2007] Ferraz-Mello, S. 2007. Canonical Perturbation Theories - Degenerate Systems and Resonance. Series: Astrophysics and Space Science Library, Vol. 345. Springer US..
  * Henrard and Murigande [1987] Henrard, J., Murigande, C. 1987. Colombo’s top. Celestial Mechanics 40, 345–366.
  * Laskar and Robutel [1993] Laskar, J., Robutel, P. 1993. The chaotic obliquity of the planets. Nature 361, 608–612.
  * Millholland and Batygin [2019] Millholland, S., Batygin, K. 2019. Excitation of Planetary Obliquities through Planet-Disk Interactions. The Astrophysical Journal 876.
  * Batygin and Adams [2013] Batygin, K., Adams, F. C. 2013. Magnetic and Gravitational Disk-Star Interactions: An Interdependence of PMS Stellar Rotation Rates and Spin-Orbit Misalignments. The Astrophysical Journal 778.
  * Lithwick and Wu [2011] Lithwick, Y., Wu, Y. 2011. Theory of Secular Chaos and Mercury’s Orbit. The Astrophysical Journal 739.
  * Batygin et al. [2015] Batygin, K., Morbidelli, A., Holman, M. J. 2015. Chaotic Disintegration of the Inner Solar System. The Astrophysical Journal 799.
  * Peale [1976] Peale, S. J. 1976. Orbital resonance in the solar system.. Annual Review of Astronomy and Astrophysics 14, 215–246.
  * Borderies and Goldreich [1984] Borderies, N., Goldreich, P. 1984. A Simple Derivation of Capture Probabilities for the J+1:J and J+2:J Orbit-Orbit Resonance Problems. Celestial Mechanics 32, 127–136.
  * Plum [2016] Plum, M. A. 2016. Beam Loss in Linacs. arXiv e-prints, arXiv:1608.02456
  * Plum [2012] Plum, M. A. 2012. Beam Loss Mechanisms in High-Intensity Linacs. Proceedings of ICFA Advanced Beam Dynamics Workshop on High-Intensity and High-Brightness Hadron Beams, 36
  * Boine-Frankenheim et al. [2001] Boine-Frankenheim, O. and 6 colleagues 2001. Emittance Coupling in High Intensity Beams Applied to the SNS Linac. 19th Particle Accelerator Conference (PAC 2001), 2902–2904.
  * Miura et al. [2016] Miura, A., Hasegawa, K., Miyao, T., Maruta, T., Liu, Y., Horioka, K. 2016. Beam profile monitor for intense, negative, hydrogen-ion beams in the J-PARC linac. Journal of Korean Physical Society 69, 1005–1013. doi:10.3938/jkps.69.1005
  * Miyao [2019] Miyao, T. 2019. Beam profile measurement with carbon nanotube wire at J-PARC Linac. Proceedings of J-PARC 10th anniversary commemorative symposium.
  * Liu et al. [2019] Liu Y., Fang Z., Futatsukawa K. and 9 colleagues 2019. Progress of J-PARC LINAC Commissioning. 10th International Particle Accelerator Conference
  * Reiser and Brown [1993] Reiser, M., Brown, N. 1993. Thermal distribution of relativistic particle beams with space charge. Physical Review Letters 71, 2911–2914.
  * Henrard [1982] Henrard, J. 1982. Capture Into Resonance - an Extension of the Use of Adiabatic Invariants. Celestial Mechanics 27, 3–22.
  * Neishtadt [1984] Neishtadt, A. I. 1984. Separation of motions in systems with a rapidly rotating phase. Prikladnaia Matematika i Mekhanika 48, 197–204.
  * Chirikov [1960] Chirikov, B. V. 1960. Resonance processes in magnetic traps. Journal of Plasma Physics 1, 253–260.
  * Chirikov [1979] Chirikov, B. V. 1979. A universal instability of many-dimensional oscillator systems. Physics Reports 52, 263–379.
  * Bruhwiler and Cary [1989] Bruhwiler, D. L., Cary, J. R. 1989. Diffusion of particles in a slowly modulated wave. Physica D Nonlinear Phenomena 40, 265–282.
  * Henrard and Henrard [1991] Henrard, J., Henrard, M. 1991. Slow crossing of a stochastic layer. Physica D Nonlinear Phenomena 54, 135–146.
  * Elskens and Escande [1991] Elskens, Y., Escande, D. F. 1991. Slowly pulsating separatrices sweep homoclinic tangles where islands must be small: an extension of classical adiabatic theory. Nonlinearity 4, 615–667.
  * Wisdom [1983] Wisdom, J. 1983. Chaotic behavior and the origin of the 3/1 Kirkwood gap. Icarus 56, 51–74.
  * Reiser [1994] Reiser, M. 1994. Theory and Design of Charged Particle Beams. Theory and Design of Charged Particle Beams, Wiley-VCH
  * Cincotta and Simó [2000] Cincotta, P. M., Simó, C. 2000. Simple tools to study global dynamics in non-axisymmetric galactic potentials - I. Astronomy and Astrophysics Supplement Series 147, 205–228.
  * Benettin et al. [1980] Benettin, G., Galgani, L., Giorgilli, A., Strelcyn, J.-M. 1980. Lyapunov characteristic exponents for smooth dynamical systems and for Hamiltonian systems - A method for computing all of them. I - Theory. II - Numerical application. Meccanica 15, 9–30.