Chaotic dynamics driven by particle-core interactions

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-14 . a) Electronic mail: kbatygin@caltech.edu; author to whom correspondence should be addressed. b) Electronic mail: batygin@lanl.gov One illustrative approach towards characterizing phase-space growth of space-charge-dominated beams rests on the particle-core interaction model of halo formation 7,9 . 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 phasespace have been routinely observed in numerical simulations of the particle-core interaction model, the phasespace 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.

A. Pulsation of the Core
The dynamics of the beam core's radius is governed by the Hamiltonian 9 : where r = x c /R is a dimensionless core radius and p 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 bothx-andŷ-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 where m is the particle mass, c is the speed of light, p is the particle momentum, I is the beam current, I c = 4π 0 m c 3 /q is the characteristic beam current, 0 is the permittivity of free space, and is the normalized beam emittance. Note that instead of b, Wangler et al. 7 parameterize H c in terms of the in-core space-charge tune-depression ratio where σ and σ 0 are the phase advances of transverse particle oscillations per focusing period, with and without space-charge forces, respectively. Moreover, notice that because H c is expressed in terms of dimensionless variables, the unit of time must also be appropriately nondimensional, and here, t = τ Ω, where Ω is the frequency of particle oscillation in a uniform focusing channel in absence of space-charge, and τ is the dimensional time. The Hamiltonian (1) possesses a global minimumwhich 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 ρ = r − 1, p ρ = 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 ρ, we have: where the characteristic frequency of core oscillation is ω = 2 (2 + b)/(1 + b), ζ = (6 + b)/(1 + b), and we have dropped a dynamically irrelevant constant term from the expression. 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 ρ 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): where ∆ is the amplitude of oscillation, and ϕ is an arbitrary phase. For consistency with previous studies 7,11,16 , we set ϕ = π such that r is minimal and p r is null at t = 0. The solution (5) is shown on a (ρ, p ρ ) plane in Figure  (1) for a series of values of ∆ ∈ (0, 0.4) and b ∈ (0.1, 10)a physically relevant parameter range that we will adopt for the remainder of the paper. We note that these bounds on ∆ and b are motivated by experiment: even the most high-intensity machines in operation are characterized by b 10, while beams with oscillation amplitude in excess of ∆ 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 H c is satisfactory for ∆ 0.4 and improves for higher values of b. In fact, this dependence can be expected, since (ζ/ω) 2 ranges from 3/2 for b → 0 to 1/2 for b 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.

B. 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.

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 7,9,16 : where u and p u refer to the particle's position and momentum, respectively. Because H 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,7,10,11,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: 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 ∆ = 0, and can even exceed unity for ∆ = 0.4 and |u| 1 + ∆. 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 ∼ 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 B 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+ρ and expand the resulting expression as a power-series in ∆. For the purposes of our analysis, it suffices to truncate this expansion at linear order in ∆ for two reasons: first, retention of higher-order terms does not generate any slowlyvarying 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 ∆, meaning that a more precise representation of (1/r 2 ) is unlikely to yield more accurate results.
With these approximations in place, the Hamiltonian takes the form: To proceed further, we introduce action-angle coordinates (ψ, Ψ), implicitly defined in terms of (u, p u ): 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 F 1 = −(1/2) u 2 tan −1 (ψ). After some rearrangement, equation (8) becomes: In order to identify the slowly-varying harmonic in equation (10), let us consider H p in the limit of vanishing b and ∆. In this regime, it is trivial to see thaṫ ψ = ∂H p /∂Ψ → 1. Simultaneously, equation (4) dictates that ω = 2 for b = 0. This means that the harmonic cos(2 ψ − ω 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 (φ, Φ) that reflect this property of the dynamics.
To accomplish this, we extend the phase-space of the system, by adding a dummy action, T , conjugate to time 24 , such that the new Hamiltonian becomes: Next, we employ a canonical contact transformation of variables, stemming from the type-2 generating function: The transformation equations then trivially yield: In preparation for canonical averaging, let us write the transformed HamiltonianH p asH p =Ȟ p +Ĥ p , wherě represents a trivially integrable kernel of the dynamics, and is a harmonic perturbation 23 .
The terms inĤ p that exhibit an explicit dependence on ξ 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, χ (in this case, the un-averaged and averaged variables are related to one another through Φ = S t χΦ , φ = S t χφ , where S t χ is the Lie series operator 24 ). Accordingly, carrying out the averaging procedure to first order, the homologic equation reads: where {, } denote the Poisson bracket. For the problem at hand, equation (16) is satisfied by: sin(4 φ + 2 ω ξ) After this reduction, particle dynamics are governed by the integrable Hamiltonian where we have dropped the now-constant action Ξ, and the averaged nature of the coordinates (Φ, φ) is implied. This expression agrees with equation (18) of ref. 16, although our derivation is different. The functional form of χ illuminates an important feature of the resonant structure of the dynamics encapsulated by Hamiltonian (7) -namely that conditions for χ 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. 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 H 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.

of Morbidelli
This interpretation is consistent with the numerical results in refs. 6 and 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 {{Ȟ p , χ}, χ}, {{{Ȟ p , χ}, χ}, χ}, 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.

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-31 , secular interactions 32,33 , as well as orbital resonances 18,34 .
In general, Andoyer Hamiltonians have the form 35 : where k is an integer greater than zero, (γ, Γ) are actionangle variables, and δ 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 phasespace, and for the problem at hand, k = 2. Moreover, the topological character of the phase-space portrait of K is regulated exclusively by δ, and it is straightforward to show that depending on the value of δ, the qualitative picture lies in one of three regimes: (i) for δ < −1, 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 δ < 0, the equilibrium point at Γ = 0 is rendered hyperbolic, and a figure-8 resonant separatrix emerges; (iii) for δ 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 Let us rewrite the particle-core interaction model in the form of equation (19), with an eye towards expressing δ in terms of b and ∆. 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 Γ = Φ/η and K =H p /η. Upon direct substitution, we require the coefficient upfront the Γ 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 η: which reduces to η = 32 ∆/3. Finally, we set K = K (2(1 + b))/(b ∆), which reduces the Hamiltonian to the form given by equation (19), with the resonance proximity parameter given by: It is illuminating to examine the limiting regimes of equation (21). In the space-charge dominated case of b 1, the resonance proximity parameter δ ≈ ( , meaning that the existence of a stable inner circulation region of phasespace occupied by the core -granted by the condition that δ > 0 -is also guaranteed as long as ∆ < 1/2. Although the possibility that δ can become negative for sufficiently large values of the beam pulsation amplitude -causing the Φ = 0 fixed to become hyperbolically Within the parameter range of interest, δ 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 ( unstable -is intriguing, we remind the reader that our entire analysis is predicated on the assumption that ∆ 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 ∆. This suggests that the potential for δ to become smaller than zero insinuated by equation (21) is an artifact of the ∆ 1 assumption. Figure (3) shows the resonance proximity parameter for ∆ = 0.1 − 0.4 as a function of b. Notably within our adopted range of beam parameters, δ is bracketed by 0 δ 2. Importantly, the fact that δ 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, ∆, 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 highintensity 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 ≈ 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 focus-ing 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 particlecore 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.

A. A Semi-Analytic Model
The electromagnetic field of a beam with a Gaussian charge distribution has the form: 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 Λ ∝ u for u r to Λ ∝ 1/u for u r. Nevertheless, a powerseries expansion Λ ≈ (b/(1 + b))(2 u/r 2 − 2 u 3 /r 4 + ...) does not yield an adequate approximation for equation (22) for u 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 Λ in closed form.
The Hamiltonian corresponding to particle motion within a Gaussian beam is: where T is a dummy action conjugate to time and Ei is the exponential integral. Changing variables in accord with canonical transformations outlined in equations (9) and (13), the above expression becomes: and we remind the reader that ξ = t.
Rather than attempting to reduce H p to a one-degreeof-freedom Hamiltonian through the Lie series approach (which entails solving the homologic equation 16 as in FIG. 5. The field of a Gaussian beam, Λ, evaluated setting the core radius to its equilibrium value of r = 1, is shown as a black curve. Polynomial approximations to Λ of various degrees are shown with multi-colored curves. Clearly, powerseries expansions fail to capture the form of Λ for u 1, even at a crude level. section II B 1), here we instead take advantage of separation of timescale inherent to the (Φ, φ) and (Ξ, ξ) degrees of freedom. In particular, becauseφ ∼ 1 − ω/2 ξ = 1, the evolution of (Φ, φ) 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 ξ: Because this procedure reducesH p to an integrable system semi-analytically, in order to illuminate the associated dynamics, we can simply project the level curves ofH p onto the (Φ, φ) plane. Figure (6) shows a sequence of phase-space portraits ofH 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 ∆ 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 H 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 ∆ = 0 (such that r = 1, which renders H p autonomous), and examine its equilibrium. In particular, it is instructive to consider the Hessian matrix of H p , evaluated at the (u, p u ) = (0, 0) fixed point. Since there is no u − p u coupling in the Hamiltonian, ∂ 2 H p /∂u ∂p u = 0. It is further trivial to see that ∂ 2 H p /∂p 2 u = 1. The determinant of the Hessian matrix is then simply 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 H p ) to a hyperbolic fixed point (corresponding to a saddle point of H 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 FIG. 6. Numerically averaged phase-space portraits ofHp, depicting the resonance in (Φ, φ), 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 δ > 0.
at |u| > 0, p u = 0 is shown in Figure (7), where we plot the phase-space portrait of H p for b = 1/3, 1, 3, 10, setting ∆ = 0. An important consequence of the appearance of this resonance is that when ∆ 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 ∆ = 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: H p = p 2 u /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 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 ∆, this layer con-nects 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.

B. 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 H 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 Λ at cubic order: 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 B 1. We  FIG. 7. Phase-space portraits of u = √ 2 Ψ sin(ψ), pu = √ 2 Ψ cos(ψ) dynamics within a Gaussian beam for ∆ = 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 ∆ > 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 < δ < 0.
substitute equation (5) for r, but to keep the model as compact as possible, we expand the resulting expression to leading order in ∆. Moreover, in spirit of a seriesexpansion approach, upon transforming to action-angle variables (Ψ, ψ) via equation (9), we neglect all harmonics proportional to Ψ 2 , retaining only those that are linear in Ψ. This reduces the Hamiltonian to the following simplified form: 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 ω. 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 fullfledged deterministic chaos. For this reason, the modulated pendulum is routinely adopted as a paradigm for stochastic dynamics 24,[47][48][49] .
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 H p upon Ψ, it is easy to see how sinusoidal variation its strength would give rise to a total of three resonances. In particular, because (1 + 2 ∆ cos(ω t)) cos(2 ψ) = cos(2 ψ) + ∆ cos(2(ψ − ω t/2)) + ∆ cos(2(ψ + ω t/2)), the fixed points of these resonances would be separated in action space (equivalently, frequency space) by ω/2. Ac-cordingly, global chaos would ensue when the half-widths of these resonant multiplets become comparable to ω/2 such that the separatrixes 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 Ψ 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 Ψ > 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 σ/σ 0 < 1/ √ 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 ψ (which only exists for b > 1 but is largely independent of ∆), the width of the resonance corresponding to φ = ψ−ω/2 depends on ∆. Consequently, a reasonable way to formulate the Chirikov criterion is to equate the distance between the stable equilibrium points of the ψ and φ resonances to the half-width of the resonance in φ. Let us now compute the relevant quantities. The equilibrium action of the unperturbed (Ψ, ψ) dynamics, [Ψ], follows FIG. 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 ψ) and the "breathing" resonance (in φ) 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.
from Hamilton's equations in the ∆ → 0 limit: From inspection, it is evident that the equilibrium value of the angle lies at [ψ] = π/2, 3π/2 such that cos(2[ψ]) = −1. The value of the unperturbed equilibrium action is then The functional form of [Ψ] is yet another reminder of the fact that chaotic dynamics can only appear for b 1: for values of b smaller than unity, resonance overlap is precluded by [Ψ] 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 [Φ] can be computed in a similar manner. Upon transforming variables according to equations (13) and dropping (i.e., "averaging away") ξ-dependent terms, Hamilton's equations give: where we have set cos(2 [φ]) = −1 by analogy with equation (29). After mild rearrangement, the distance between resonant equilibria in action-space reads: The remaining step is to compute the resonance halfwidth δΦ. In the pendulum approximation -which is justifiable for the φ-resonance as long as δΦ [Φ] 50 -we can substitute [Φ] for the value of the action upfront the cosine term in H p . This gives 24 : Recalling that the existence of the resonance in the (ψ, Ψ) degree of freedom requires that b > 1, the Chirikov criterion ([Φ] − [Ψ])/δΦ 1 can now be written as Although it is possible to express b as a function of ∆ (or alternatively ∆ 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 (∆, b) plane. Upon inspection, two physically important characteristics of this threshold are immediately evident. First, for diminishing values of ∆, progressively higher values of b are necessary to drive chaotic motion. This relationship stems from the ∼ 1/ √ ∆ dependence in equation (34), and simply insinuates that in the limit of vanishing ∆, H p reduces to an integrable form, implying perfectly regular dynamics.
Second, for ∆ > ( √ 3 − 1)/3 ≈ 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 ∼ 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.

C. 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.
To populate each surface of section, we chose 30 initial conditions lying on the u−axis, spanning the range u 0 ∈ (0, 1 − ∆). 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ṕ u = ,ú = u 0 + , where = 10 −6 ) through integration of linearized equations of motion (see appendix B). Keep-ing track of the distance between the primary and shadow trajectories in phase-space, ν = (p u −ṕ u ) 2 + (u −ú) 2 , we compute the Lyapunov coefficient, λ, of the trajectory. To be more specific, we computed the so-called mean exponential growth factor of nearby orbits 52 : which known to converge more rapidly than the standard renormalization method for computing λ 53 . To this end, we note that the running time-average of Y tends towards Y ∼ 2 if the trajectory is regular, and towardsȲ ∼ λ t/2 if the trajectory is chaotic.
If all initial conditions on the surface of section yield quasi-period dynamics, we record λ = 0 for the given combination of ∆ and b. If, on the other hand, some trajectories result in chaotic motion, we record the mean value of λ, 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 λ 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.
Beyond the binary question of whether or not a given combination of beam parameters yields chaotic motion or not, the value of λ itself tells an intriguing story. That is, the distribution of λ is clearly bimodal, with moderate values of ∆ and b corresponding to λ ∼ 0.1 − 0.3 and large values of ∆ and b corresponding to λ 0.5. Crucially, the boundary between these two domains appears rather sharp, and corresponds to the onset of chaotic linkage of the separatrixes of φ and ψ 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 φresonance (computed using semi-analytical averaging as in section III A) 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.
FIG. 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, λ, 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 ≈ 2.3 for the lowest numerically tested beam pulsation amplitude of ∆ = 0.01, this curve turns sharply towards b → ∞ as ∆ → 0, where the dynamics becomes integrable. Finally, the threshold for chaotic mixing between the core and the halo is outlined with a black curve.

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 phasespace 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 B and III A 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 spacecharge 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-∞ 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 ∆, 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.

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