# Hamiltonian chaos for one particle with two waves: Self-consistent dynamics

Matheus J. Lazarotto*

Aix-Marseille Université, CNRS, UMR 7345 PIIM, F-13397, Marseille cedex 13, France and
Instituto de Física, Universidade de São Paulo, Rua do Matão 1371, São Paulo 05508-090, Brazil

Iberê L. Caldas†

Instituto de Física, Universidade de São Paulo, Rua do Matão 1371, São Paulo 05508-090, Brazil

Yves Elskens‡

Aix-Marseille Université, CNRS, UMR 7345 PIIM, F-13397, Marseille cedex 13, France

(Dated: January 28, 2025)

A simple model of wave-particle interaction is studied in its self-consistent form, that is, where the particles are allowed to feedback on the waves dynamics. We focus on the configurations of locked solutions (equilibria) and how the energy-momentum exchange mechanism induces chaos in the model. As we explore the system, we analyse the mathematical structure that gives rise to locked states and how the model’s non-linearity enables multiple equilibrium wave amplitudes for equilibrium. We also explain the predominance of regularity as we vary the control parameters and the mechanism behind the emergence of chaos under limited parameter choices.

Keywords: Hamiltonian System; Chaos; Plasma Physics; Wave-particle interaction

## I. INTRODUCTION

When describing a plasma on a microscopic scale, a fundamental process is the interaction between its particles (electrons and ions) and waves propagating in the plasma [1, 2]. In the particular limit of collisionless plasmas, that is when the time scale of particle ‘collisions’ is much longer than the time scale of the interaction with waves, the energy and momentum exchange between particles and waves becomes relevant to the point of generating instabilities, favoring wave amplification or particle acceleration [3].

For instance, this scenario can be experimentally reproduced in traveling wave tubes (TWT), where a beam of electrons interacts with electrostatic waves in a controllable environment [4, 5]. This setup enables the analysis of chaotic interactions due to the energy-momentum exchange that may lead to instabilities and further turbulence in plasmas, which in turn correlates to the emergence of chaotic regions in phase space. This chaotic transition can be caused by resonance overlap of waves, nonlinear synchronization by a single nonresonant wave or the ‘devil’s staircase’ mechanism, causing a spread of the velocity distribution of the beam [6, 7]. Stability islands on the other hand can correlate with coherent particle acceleration (beam trapping), and the breaking of barriers (invariant lines in the one-particle phase space) allowing for new ways of chaos control.

When modeling such scenarios, the Langmuir waves, that is the collective vibration of electrons in presence of much heavier ions that neutralize the system, are found

to be coupled with quasi-resonant particles (i.e. electron’s velocity being close to waves phase velocity) [8]. Thus, for simplicity, the usual Hamiltonian models initially consider only the effect of waves on the particles. To give a step ahead, the re-interaction of the particles with the waves can be taken into account, therewith making the dynamics self-consistent [1, 2]. Consequently, electrons (and ions) do not play just the role of test particles but instead promote energy-momentum exchange with the waves, allowing for the aforementioned nonlinear phenomena. To model it, one can consider *M* independent Langmuir waves behaving as harmonic oscillators and *N* quasi-resonant free charged particles, with the addition of a coupling term for each wave-particle pair [9–11].

Gomes *et al.* [12] considered two of the simplest cases, namely that of one wave and one particle (*M* = 1, *N* = 1), which is analytically treatable and shows the electron either being trapped (strong resonance) or passing through wave potential wells [13, 14], and the case of one wave interacting with two particles (*M* = 1, *N* = 2), that can present chaotic behavior. We direct the reader to the work of Gomes *et al.* [12] for a simple yet broad introduction on the experimental context, modeling and applications of wave-particle dynamics.

In this sense, our work is complementary to Gomes’ work, for we consider the remaining case of two waves interacting with one particle (*M* = 2, *N* = 1), aiming to a simple description of the interaction between waves and particles in contexts such as TWTs. By studying these simplest cases, one may shed light on the microscopic mechanisms that induce chaos in the macroscopic system and how to prevent (or allow) it. In the current work, for instance, it is shown that equilibrium solutions only come in locked configurations allowing for multiple unstable wave amplitudes combinations. The emergence of chaos in the system is predominantly due to separatrix

<footnote>* matheus.jean.l@hotmail.com</footnote>

<footnote>† ibere@if.usp.br</footnote>

<footnote>‡ yves.elskens@univ-amu.fr</footnote>
chaos related to the parameter limit where the non-linear coupling produces energy-momentum exchange.

In what follows, Sec. II A starts by presenting the generic wave-particle self-consistent Hamiltonian model and Sec. II B the particular case studied here (for $M = 2, N = 1$), along with its mathematical simplification to reduce the number of degrees of freedom and number of parameters. Sec. III presents a global analysis of the equilibrium solutions in parameter space. Sec. IV discusses the emergence of chaos in phase space and the integrable limits of the system. Appendix sections provide more details on the model parameter simplifications, the stability of equilibrium solutions, and the decoupling limit of the system equations of motion.

## II. THE WAVE-PARTICLE HAMILTONIAN

### A. The generic model

Within the context of Langmuir waves interacting with ions in plasmas, one can model the self-consistent dynamics of $N$ identical particles traveling in a periodically bounded interval of length $L$ via a coupling with $M$ longitudinal waves. Each particle is described by generalized coordinate-momentum pairs $(x_i, p_i)$ and each wave written in phasor formulation as $Z_j = X_j + iY_j = \sqrt{2I_j} e^{-i\theta_j}$, with either cartesian $(X_j, Y_j)$ or polar $(\theta_j, I_j)$ canonical coordinates. Within this interval, wave numbers are given by $k_j = \frac{2\pi j}{L}$, for $j \in \mathbb{Z}$, and natural frequencies by $\omega_{0j}$. The system's Hamiltonian is then [1, 2]

$$H_{\text{sc}}^{N,M} = \sum_{i=1}^{N} \frac{p_i^2}{2m} + \sum_{j=1}^{M} \omega_{0j} \frac{X_j^2 + Y_j^2}{2} + \epsilon \sum_{i=1}^{N} \sum_{j=1}^{M} \frac{\beta_j}{k_j} (Y_j \sin(k_j x_i) - X_j \cos(k_j x_i)), \quad (1)$$

or equivalently in polar-phasor coordinates

$$H_{\text{sc}}^{N,M} = \sum_{i=1}^{N} \frac{p_i^2}{2m} + \sum_{j=1}^{M} \omega_{0j} I_j - \epsilon \sum_{i=1}^{N} \sum_{j=1}^{M} \frac{\beta_j}{k_j} \sqrt{2I_j} \cos(k_j x_i - \theta_j), \quad (2)$$

with $\beta_j$ as a coupling constant with the jth-wave and $\epsilon$ an overall coupling scale. For the wave degrees of freedom, $X_j$ and $\theta_j$ are treated as coordinates with $Y_j$ and $I_j$ as their respective conjugate momenta.

The Hamiltonian $H_{\text{sc}}^{N,M}$ comprises three contributions: the free motion (kinetic energy) of particles ($m > 0$); the harmonic oscillation of waves ($\omega_{0j} > 0$); and the coupling between particles and waves. Furthermore, the Hamiltonian $H_{\text{sc}}^{N,M}$ is invariant under translations in time and in space, so that the total energy $E = H_{\text{sc}}^{N,M}$ and the

total momentum $P = \sum_{i=1}^{N} p_i + \sum_{j=1}^{M} k_j I_j$ are conserved. The latter constant reveals that the growth or decay of a wave with positive phase velocity ($\frac{\omega_{0j}}{k_j} > 0$) is directly balanced with the slowing down or acceleration of particles. On interchanging the time and space variables, this Hamiltonian also captures the key physics in the dynamics of traveling wave tubes [2, 4]. For instance, in the paradigmatic two-wave scenario with a slaved particle, given parameters $\epsilon\beta_1 = \epsilon\beta_2 = \beta, k_1 = k_2 = k, m$, the Chirikov overlap parameter [15, 16] is $s = 2\sqrt{\beta k/m}((2I_1)^{1/4} + (2I_2)^{1/4})/|\omega_1 - \omega_2|$.

### B. The single particle with two waves

To address the key physics of the system, one can start looking at the simpler scenarios, namely those with few waves and particles interacting with each other. The simplest case, that of $M = N = 1$, is integrable and was treated before in [12–14]; thanks to its integrability, it provides reference information for the more complex cases where chaos emerges. With the addition of a particle ($N = 2, M = 1$), as studied by Gomes et al. [12], the non-linear term in Hamiltonian (1) implies the emergence of chaos in phase space not only around separatrices but also near elliptic points.

Here, we consider the ‘mirror’ case of one particle coupled to two waves, that is the case with $N = 1, M = 2$ in Hamiltonian (1) (2). The system thus has 7 free parameters: $k_1, k_2, \epsilon\beta_1, \epsilon\beta_2, \omega_{01}, \omega_{02}$ and $m$ (along with 3 free scales, leaving 4 dimensionless free parameters). For simplicity, the $\epsilon\beta_j$ are selected to be equal and unity: $\epsilon\beta_1 = \epsilon\beta_2 = 1$, as well as the particle mass: $m = 1$, which can be set without loss of generality by rescaling variables (see Appendix A for more details). Given the translational invariance, only the relative positions $\theta_j - k_j x$ are dynamically relevant; thus, for commensurate $k_1$ and $k_2$, the particle position $x$ can meet periodic boundary conditions, $x \in [0, L]$, compatible with both waves. For incommensurate wave numbers, the position space is necessarily the full real line. By setting equal wave numbers and given a proper rescaling of space, we select $k_1 = k_2 = 1$ (in order to keep free phase-velocities, the frequencies $\omega_{0j}$ are kept different).

The resulting 3-degrees-of-freedom Hamiltonian is autonomous and translation-invariant; in cartesian-phasor coordinates it reads

$$H_{\text{sc}}^{1,2} = \frac{p^2}{2} + \omega_{01} \frac{X_1^2 + Y_1^2}{2} + \omega_{02} \frac{X_2^2 + Y_2^2}{2} + (Y_1 + Y_2) \sin(x) - (X_1 + X_2) \cos(x), \quad (3)$$

and likewise in polar-phasor form

$$H_{\text{sc}}^{1,2} = \frac{p^2}{2} + \omega_{01} I_1 + \omega_{02} I_2 - \sqrt{2I_1} \cos(x - \theta_1) - \sqrt{2I_2} \cos(x - \theta_2). \quad (4)$$
To further simplify the model and better visualize solutions in phase space, it is convenient to reduce the number of degrees of freedom, which is obtained with the following canonical trasformations. The particle's degree of freedom can be suppressed with the aid of the total momentum constant P using $p = P - I_1 - I_2$. Besides, a Galileo transform to a reference frame traveling with the first wave, provided by the generating function

$$G(x, \theta_1, \theta_2, \bar{p}, \bar{I}_1, \bar{I}_2) = (x - \omega_{01}t)(\bar{p} + \omega_{01}) + (\theta_1 - \omega_{01}t)\bar{I}_1 + (\theta_2 - \omega_{01}t)\bar{I}_2 + \frac{\omega_{01}^2 t}{2},$$

and the definition of the relative angle $\phi_j := \bar{\theta}_j - \bar{x} = \theta_j - x$, provided by the generating function

$$F(\bar{x}, \bar{\theta}_1, \bar{\theta}_2, p', I'_1, I'_2) = I'_1(\bar{\theta}_1 - \bar{x}) + I'_2(\bar{\theta}_2 - \bar{x}) + p'\bar{x},$$

simplify Hamiltonian (4) to

$$H' = \frac{(I'_1 + I'_2)^2}{2} - \bar{P} (I'_1 + I'_2) + \Delta_{\omega} I'_2 - \sqrt{2I'_1} \cos(\phi_1) - \sqrt{2I'_2} \cos(\phi_2), \quad (5)$$

where prime symbols denote the final transformed variables $(x, \theta_j, p, I_j) \stackrel{G}{\mapsto} (\bar{x}, \bar{\theta}_j, \bar{p}, \bar{I}_j) \stackrel{F}{\mapsto} (x', \phi_j, p', I'_j)$ and will be dropped from now on.

From old to new variables, the amplitudes are unmodified $(I'_j = I_j)$ and the new angles are the relative position between particle and wave $\phi_j = \theta_j - x$. The new total momentum is now $\bar{P} = P - \omega_{01}$ and we define the detuning parameter between the waves: $\Delta_{\omega} := \omega_2 - \omega_1$, eliminating one of the frequencies. This final form has 2 degrees of freedom $(\phi_1, \phi_2)$ and 3 control parameters $(H, \bar{P}, \Delta_{\omega})$; the Hamiltonian $H$ itself being a constant of the motion. The frequency $\omega_{01}$ can be set to zero without loss of generality.

Complementary, in cartesian form the Hamiltonian (3) in simplified form reads

$$H' = \frac{1}{8} (u_1^2 + v_1^2 + u_2^2 + v_2^2)^2 - \frac{\bar{P}}{2} (u_1^2 + v_1^2 + u_2^2 + v_2^2) + \frac{\Delta_{\omega}}{2} (u_2^2 + v_2^2) - u_1 - u_2. \quad (6)$$

where

$$u_i = \sqrt{2I_i} \cos(\phi_i), \quad v_i = \sqrt{2I_i} \sin(\phi_i), \quad (7)$$

for $i = 1, 2$.

It may be of interest to keep both forms since cartesian coordinates provide a smooth transition for the limit $I_j \to 0$, though the polar ones provide a more intuitive picture of the dynamics.

It is worth mentioning that Hamiltonian (6) is fourth degree in both $u_i$ and $v_i$, making the visualization of phase space via Poincaré sections less straightforward. Since the intersection of trajectories with the section may

yield two positive (or negative) roots, the usual method of fixing a position (say $u_2 = 0$) and a sign condition on the conjugate momentum (say $v_2 > 0$), is not enough to produce a uniquely defined map. Since extra roots may exist and appear over the points of map made by the first set of roots, invariant curves appear overlapping. A solution can be achieved by filtering these two sets of roots with an extra condition as $\dot{u}_2 > 0$.

### III. LOCKED SOLUTIONS

As a reference for the system's global dynamics, one can look at the equilibrium solutions (i.e. $(\dot{I}_1, \dot{I}_2, \dot{\phi}_1, \dot{\phi}_2) = \mathbf{0}$) of the equations of motion provided by Hamiltonian (5)

$$\begin{align}
\dot{I}_1 &= -\partial_{\phi_1} H' = -\sqrt{2I_1} \sin(\phi_1), \nonumber \\
\dot{I}_2 &= -\partial_{\phi_2} H' = -\sqrt{2I_2} \sin(\phi_2), \nonumber \\
\dot{\phi}_1 &= \partial_{I_1} H' = (I_1 + I_2) - \bar{P} - \frac{\cos(\phi_1)}{\sqrt{2I_1}}, \tag{8} \\
\dot{\phi}_2 &= \partial_{I_2} H' = (I_1 + I_2) - \bar{P} + \Delta_{\omega} - \frac{\cos(\phi_2)}{\sqrt{2I_2}}. \nonumber
\end{align}$$

In the present context, equilibrium solutions are also referred to as locked states, given that, when in equilibrium, the amplitude equations $(\dot{I}_1, \dot{I}_2)$ yield the condition that the relative phases are constant $(\phi_i^* \in \{0, \pi\}$, for $i = 1, 2$) meaning that both waves and the particle travel at the same velocity with the latter constrained to be at either the minimum or maximum of one of the waves while the other is in phase or anti-phase. It is worth mentioning that the last two equations in (8) prevent the existence of zero-wave solutions $(I_1(t) = I_2(t) = 0)$ for both waves simultaneously for generic $(\bar{P}, \Delta_{\omega})$. For single null waves, i.e. $I_1^* = 0 \neq I_2^*$ (or vice-versa), it is possible to obtain such solutions for $\Delta_{\omega} = 0$. In such case, the null wave has indefinite phase $\phi_i$ and the remaining one must have fixed amplitude given by the total momentum $\bar{P}$. Other considerations on the $\Delta_{\omega} = 0$ case are discussed in Appendix B.

From the locked relative phases $\phi_1^* = n_1\pi$ and $\phi_2^* = n_2\pi$, for $n_i = 0, 1$ and $i = 1, 2$, the equilibrium amplitudes $I_1^*$ and $I_2^*$ thus must satisfy

$$\frac{\cos(n_2\pi)}{\sqrt{2I_2^*}} - \frac{\cos(n_1\pi)}{\sqrt{2I_1^*}} = \Delta_{\omega} \quad (9)$$

Therewith, the relation between amplitudes enables one to find $(I_1^*, I_2^*)$ by substituting eq. (9) into $\dot{I}_1 = 0$, yielding

$$\begin{align}
R(I_2, n_2) = &I_2 \sqrt{2I_2} \left(1 + \left((-1)^{n_2} - \Delta_{\omega} \sqrt{2I_2}\right)^2\right) + \nonumber \\
&\left(\left(\Delta_{\omega} - \bar{P}\right) \sqrt{2I_2} - (-1)^{n_2}\right) \left((-1)^{n_2} - \Delta_{\omega} \sqrt{2I_2}\right)^2 \tag{10}
\end{align}$$
Solutions are then given by the condition $R(I_2^*, n_2) = 0$, with each one amounting for two locked states $(I_1^*, I_2^*, n_1\pi, n_2\pi)$, namely one with $n_1 = 0$ and the other with $n_1 = 1$. The existence of roots to condition (10) was evaluated numerically by scanning over the parameter space $(\bar{P}, \Delta_\omega)$, with the result shown in figure 1 for both cases $n_2 = 0$ and $n_2 = 1$.

<figure><img src="image_1.png" /><figcaption>FIG. 1. Color map with the number of equilibrium solutions for different $n_2$. (Left) case $\phi_2^* = 0$; (Right) case $\phi_2^* = \pi$.</figcaption></figure>

<figure><img src="image_2.png" /><figcaption>FIG. 2. Surface plot of $R(I_2^*, \bar{P}, \Delta_\omega) = 0$ with its contour lines at different constant values of $\bar{P}$. Dashed (resp. bold) lines correspond to $\bar{P} < 0$ (resp. $\bar{P} > 0$). (Left) $\phi_2^* = 0$; (Right) $\phi_2^* = \pi$.</figcaption></figure>

In order to clarify the dependence of the number of roots as the control parameters change, figure 2 shows the surface $R(I_2^*) = 0$ along with the $I_2^*$ axis. In this plot, roots will exist whenever $I_2^*$ crosses the surface for fixed control parameters.

For either case of $\phi_2^*$, the set of solutions is composed of two non-intersecting surfaces. One of them has a cup-like form, which asymptotically converges to $I_2^* \to 0$ as $\Delta_\omega \to \infty(-\infty)$ for $\phi_2^* = 0$ ($\phi_2^* = \pi$). The second surface branch, completely detached from the cup, presents different concavity in each scenario of $\phi_2^*$. When $\phi_2^* = 0$, the branch is convex and therefore allows for one crossing throughout all parameter space, whereas for $\phi_2^* = \pi$ the concave surface allows only for crossings up to the limit $\Delta_\omega \lesssim 0$. The appearance of 2 extra roots, in both cases, is due to the cup-like branch, since any intersection with

it implies two extra solutions, although these are only possible for $\bar{P} \gtrsim 1$. It is worth pointing out that the reasoning above does not consider any constraint on the Hamiltonian $H$ value.

From both figures 1 and 2, the case $\phi_2^* = 0$ guarantees that it is always possible that the particle is placed at the minimum of one of the waves and that it is locked to travel at the same speed, regardless of $\bar{P}, \Delta_\omega$. Moreover, if the total momentum of the system is large enough ($\bar{P} \gtrsim 1$), new locked states become possible in which the waves may have enhanced or suppressed amplitude to provide balance of forces acting on the particle.

In addition to the solutions themselves, the stability is promptly obtained from the system's jacobian, which yields a simple biquadratic form (see Appendix C). Figures 3 to 5 show bifurcation diagrams over the contour lines shown in figure 2; in them, stable points (eigenvalues $\lambda = ib$, for $b \in \mathbb{R}_{\neq 0}$) are colored blue, purely unstable points ($\lambda \in \mathbb{R}_{\neq 0}$) are colored red, and spirally unstable ($\lambda = a + ib$ for $a, b \in \mathbb{R}_{\neq 0}$) are colored green.

When analyzing stability, we now consider the combinations of locked phases including $n_1$, thus resulting in four possible scenarios of $(n_1, n_2)$. For $(n_1, n_2) = (0, 0)$, all solutions are stable.

For $(n_1, n_2) = (0, 1)$, figure 3 shows that the concave surface has its lower (upper) branch always stable (unstable). This indicates that at the expense of suppressing the second wave momentum, it is possible for wave 1 to carry the particle at its minimum. Similarly, as the total momentum increases and the cup bifurcation takes place, the upper (lower) half of the branch is unstable (stable). For $\bar{P} > 0$, the presence of spirally unstable points amidst stable ones, in both branches, prevents the second wave amplitude growth while the particle is locked to wave 1.

For $(n_1, n_2) = (1, 0)$, figure 4 shows that higher values of $I_2^*$ are now stable, indicating the opposite behavior than the one found previously for $(n_1, n_2) = (0, 1)$. Now, with the increase of the second wave amplitude, the stability of its minimum allows for the particle to be placed at the maximum of the first wave at the same time. Again, as the total momentum increases, stable branches become spirally unstable, but now preventing the locking for small $I_2^*$ amplitudes, with small stability windows, particularly for $\bar{P} = 0$ and $\bar{P} = 8$. Also, the cup-like surface provides a stable branch for detuning $\Delta_\omega$ higher enough.

For $(n_1, n_2) = (1, 1)$, with the particle placed at the maxima of both waves simultaneously, the expected instability is found for most of the parameter space. However, figure 5 shows that for total momentum negative enough, e.g. $\bar{P} = -4$, or high enough, e.g. $\bar{P} = 4$, branches with low $I_2^*$ are stable, with no spiral instability found. This indicates that at the expense of the second wave attenuation, the particle can be kept at an unexpected unstable position.
<figure><img src="image_3.png" /><figcaption>FIG. 3. Bifurcation diagram for $n_1 = 0, n_2 = 1$ for different values of fixed $\bar{P}$, corresponding to different horizontal lines from the right frame in figure 2. Colours indicate stability, with purely unstable points in red, stable points in blue, and spirally unstable in green.</figcaption></figure>

<figure><img src="image_4.png" /><figcaption>FIG. 4. Bifurcation diagram for $n_1 = 1, n_2 = 0$ for different values of fixed $\bar{P}$, corresponding to different horizontal lines from the left frame in figure 2. Colours indicate stability, with purely unstable points in red, stable points in blue, and spirally unstable in green.</figcaption></figure>

## IV. POINCARÉ SECTIONS AND CHAOS

Once the 3-degrees-of-freedom Hamiltonian (3) (or in polar form (4)) is reduced to the equivalent 2-degree-of-freedom form (6) (or in polar form (5)), trajectories can be visualized in phase space. To this end, we made use of the Poincaré section defined as

$$ \Sigma = \{ (\phi_1, \phi_2, I_1, I_2) \in \mathbb{R}^4; \phi_1 = \phi_2; \dot{\phi}_2 > 0 \}. \quad (11) $$

From the selected Poincaré section, we aimed to measure the area of phase space occupied either by stability islands or chaotic regions as a function of the control parameters $(\bar{P}, \Delta_\omega, H)$, once fluctuations are expected given the nonlinear nature of the model. For this purpose, the chaotic/regular areas were measured over the section $\Sigma$ (eq. (11)) via a smaller alignment index (SALI) method. We refrain from detailing the method here and direct the reader to the original reference as presented by Ch. Skokos [17, 18]. Briefly, the algorithm integrates a single trajectory to numerically determine whether it is stable or chaotic, along with two deviation vectors evolved in tangent space according to the linearized equations of motion. Based on the local linearized geometry of phase space around the trajectory, the two deviation vectors align with each other along the unstable manifold direction in case the trajectory is chaotic. On the other hand, in case it is stable, the local geometry is a plane tangent to a torus to which the vectors become parallel while keeping their relative angle non null. Therefore, the orbit is deemed chaotic if the deviation vectors align or stable in case they do not.

Running the SALI algorithm for each trajectory on a grid of initial conditions over the Poincaré section $\Sigma$ for a given set $(\bar{P}, \Delta_\omega, H)$ enables one to evaluate the percentage of stable area in phase space. By repeating this while scanning the parameters $(\bar{P}, \Delta_\omega)$, we obtain the profiles seen in figure 6, where each panel corresponds to a different energy $H$ value. It is worth noting that Hamiltonian (6) defines compact energy surfaces in the 4-dimensional phase space, thus ensuring that the fractional area $A$ is well defined.

At first, the predominance of stability in the system
<figure><img src="image_5.png" /><figcaption>FIG. 5. Bifurcation diagram for $n_1 = 1, n_2 = 1$ for different values of fixed $\bar{P}$, corresponding to different horizontal lines from the right frame in figure 2. Colours indicate stability, with purely unstable points in red, stable points in blue, and spirally unstable in green.</figcaption></figure>

becomes clear as, for most of the parameter space, chaotic area is nearly zero (dark brown regions – $A \approx 1$). This should hold beyond the parameter space limits chosen here ($\bar{P}, \Delta_\omega \in [-12, 12] \times [-25, 25]$) as for either large

<figure><img src="image_6.png" /><figcaption>FIG. 6. Color map of the stable area portion in parameter space for multiple energies $H$. In the color bar, total chaos is given by $A=0$ and total regularity by $A=1$. Grid size is 75 $\times$ 75.</figcaption></figure>
$|\bar{P}|$ or $|\Delta_{\omega}|$, the equations of motion become

$$ 
\begin{aligned}
\dot{I}_1 &\approx -\sqrt{2I_1} \sin(\phi_1), \\
\dot{I}_2 &\approx -\sqrt{2I_2} \sin(\phi_2), \\
\dot{\phi}_1 &\approx \pm |\bar{P}| - \frac{\cos(\phi_1)}{\sqrt{2I_1}}, \\
\dot{\phi}_2 &\approx \pm |\bar{P}| \pm |\Delta_{\omega}| - \frac{\cos(\phi_2)}{\sqrt{2I_2}};
\end{aligned} 
\quad (12) $$

noting that $I_1$ and $I_2$ are limited by the total energy $E = H$ in eq. (5). If either $I_1$ or $I_2$ is small (in view of the denominator for $\dot{\phi}_i$), their amplitude have limited growth rate once $\phi_i$ will evolve rapidly and the average increment in $\dot{I}_1$ is zero or small, therefore generating a stable scenario. In case only $|\Delta_{\omega}|$ is large, $\phi_2$ and $I_2$ become regular as well, resulting in a regular evolution for $\phi_1$. Therefore, for high enough values of $|\bar{P}|$ and $|\Delta_{\omega}|$ stability should still dominate phase space, as the dynamics of waves uncouples and the system reduces to the paradigmatic $M = N = 1$ case which is analytically treatable [12–14]. More details on the dynamics for weak coupling are given in Appendix D.

Despite this predominance of stability, a prominent chaotic region appears at lower values of total momentum and frequency detuning, which decreases for large energy and eventually vanishes for $H \approx 10$ (figure 6). Similarly to what was previously discussed, the suppression of chaos for increasing energy is due to the dominance of the waves amplitude energy term in Hamiltonian (5). For high $H$, the spatial coupling terms become negligible, as their amplitude scales linearly with the waves amplitude whereas the ‘kinetic’ term scales quadratically. As a result, the systems approaches a free wave limit.

The emergence of these chaotic regions in phase space was seen to be of a common type, namely where chaos appears over and around separatrices, as illustrated by the series of portraits in figures 7 and 8, for $H = 0.5$. Figure 7 shows portraits for varying $\bar{P}$, in a parameter range over the first quadrant ($\bar{P}, \Delta_{\omega} > 0$) in figure 6, while figure 8 shows portraits for varying $\Delta_{\omega}$ for values in the third quadrant ($\bar{P}, \Delta_{\omega} < 0$), as the system goes from regions with full regularity ($A = 1.0$) through regions of co-existence with chaos ($A < 1.0$). Normally, chaos occurs within limited regions of phase space and it coexists with invariant lines that are often seen for higher values of amplitude $I_1$.

In more detail, the chaotic layer as shown in figure 7 emerges around the main island centered at $\phi_1 = 0$ (panel A), and increases as secondary bifurcations disrupt the island’s invariant lines. Simultaneously, a pitchfork bifurcation further promotes chaos within the main island but from the inside out (panels B, C and D). In panel C, it is seen that the chaotic domain, although fully connected, has a slower chaotic regime around the bifurcated pair of islands, with stickiness concentrated around this area. Finally, panels E and F show that the secondary resonant islands disappear, leaving a rather uniform chaotic

sea (panel E) that further disappears, returning to full regularity (panel F). In a slightly different manner, the transition shown in figure 8 corresponds to a chaotic layer emerging in between invariant regular sets (panels B and C), although with the same process of secondary bifurcations taking place.

<figure><img src="image_7.png" /><figcaption>FIG. 7. Phase space portraits over the Poincaré section $\Sigma$ (eq. (11)). A) $\bar{P} = 0.5$, B) $\bar{P} = 1.0$, C) $\bar{P} = 1.5$, D) $\bar{P} = 1.75$, E) $\bar{P} = 3.0$, F) $\bar{P} = 5.0$. In all cases, $H = 0.5$ and $\Delta_{\omega} = 2.0$.</figcaption></figure>

V. CONCLUSIONS

In the simpler scenario of one particle interacting with two waves, although initially a 3 degrees-of-freedom system, we show that simple canonical transformations yield a 2 degrees-of-freedom Hamiltonian with clear connection to the original one. On top of the transformations, the only restrictive hypotheses made were the equal coupling parameter between each wave and the particle, and the equal wavenumber for both waves. The other assumptions made, i.e. the unit mass and unit coupling coefficients, can be made by rescaling variables without loss of generality.

For this simplified Hamiltonian, it was shown that the number of equilibrium solutions has a non-trivial dependence on the control parameters, namely the total momentum $\bar{P}$ and frequency detuning $\Delta_{\omega} = \omega_2 - \omega_1$. These equilibria correspond to any locked configuration, where
<figure><img src="image_8.png" /><figcaption>FIG. 8. Phase space portraits over the Poincaré section Σ (eq. (11)). A) Δ<sub>ω</sub> = −0.5, B) Δ<sub>ω</sub> = −4.0, C) Δ<sub>ω</sub> = −6.0, D) Δ<sub>ω</sub> = −15.0. In all cases, H = 0.5 and P̄ = −1.0.</figcaption></figure>

the particle travels along with the waves while aligned with their extremum points. For any of these spatially locked states, the nonlinear nature of the wave-particle coupling can yield up to 4 combinations of wave amplitudes allowing for equilibrium. Moreover, the highest number of equilibria combinations was found for P̄ > 0, indicating that the waves require a minimum amount of momentum in order to carry the particle along with them. This comes as a consequence of the topology of the waves amplitude equilibrium function as two disjoint surfaces with concavity based on the type of locked position selected.

Despite the interaction between one particle and two waves being nearly integrable throughout most of the domain of its control parameters, the auto-consistent coupling term added in our model was able to significantly induce chaos. In the limit where both the total momentum P̄ and the frequency detuning are small (up to order of magnitude P̄ ≈ 10), the energy exchange between the waves, as mediated by the particle, is enough to destabilize trajectories. This possibility of energy-momentum exchange promoted by the non-linear coupling prevents a dynamics where the particle is either free or enslaved to one of the waves, as in the cases for large |P̄| or |Δ<sub>ω</sub>|. The emergence of chaos was seen around separatrices and higher order bifurcations, as usually observed in a KAM description. However, an ergodic limit, that could be related to a destabilization of the particle motion, was not found in the equal wave limit tested here, which could naturally be investigated as an extension of the current work.

Beyond the first exploratory results shown here, the extent to which the coupling strength parameters β<sub>i</sub> can differ in order to promote enough energy-momentum ex-

change to produce chaos is still unclear. In the same way, a difference in wavenumbers could induce chaotic behavior, although the analysis in this case becomes cumbersome as the equations lose spatial periodicity and the number of free parameters increases. More generally, one could evaluate how to extend the energy-momentum exchange mechanism to the limit of many particles and many waves, thus approaching the complete model.

ACKNOWLEDGMENTS

We acknowledge the financial support from the scientific agencies: São Paulo Research Foundation (FAPESP) under Grant No. 2018/03211-6; Conselho Nacional de Desenvolvimento Científico e Tecnológico (CNPq) under Grants No. 200898/2022-1 and 304616/2021-4. Coordenação de Aperfeiçoamento de Pessoal de Nível Superior (CAPES) and Comité Français d'Évaluation de la Coopération Universitaire et Scientifique avec le Brésil (COFECUB) under Grants CAPES/COFECUB 8881.143103/2017-1 and COFECUB 40273QA-Ph908/18. Centre de Calcul Intensif d'Aix-Marseille is also acknowledged for granting access to its high-performance computing resources.

Y. Elskens is grateful to Serge Aubry for his fascinating lectures in 1985 at Beg-Rohu, which revealed more facets of Hamiltonian dynamics and its broader range of applications, and for the enthusiasm Serge always kindly shared with his audience.

Appendix A: Parameter reduction via scalings

In the simplification step of Hamiltonian (2) to (4), we impose equal wave number for both waves (k<sub>1</sub> = k<sub>2</sub> = k) as well as the coupling parameters with the particle (β<sub>1</sub> = β<sub>2</sub>), allowing one to set ϵβ<sub>1</sub> = ϵβ<sub>2</sub> = β. To eliminate k from within the coupling term, the position can be normalized to x' = kx. For the remaining parameters, we set the global scaling factor

$$
\alpha = (\beta^2/m)^{1/3} \qquad (A1)
$$

and rescale variables according to

$$
\begin{align}
t' &= \alpha t, \nonumber \\
\omega_i' &= \alpha^{-1} \omega_i, \tag{A2}
\end{align}
$$

$$
p'/p = I_1'/I_1 = I_2'/I_2 = P'/P = (m\beta)^{-2/3},
$$

and

$$
H' = (m\beta^4)^{-1/3} H. \qquad (A3)
$$

The system in primed variables obtained from (2) is the reduced model (4) (for N = 1 and M = 2). This scaling also acts on the action so that S' = (mβ)<sup>-2/3</sup>S.
## Appendix B: Degenerate limit of null detuning

In Hamiltonian (3) or (4), the case $\omega_1 = \omega_2$ is degenerate (or equivalently for $\Delta_\omega = 0$ in Hamiltonian (6) or (5)), as both waves have the same wavenumber and the same phase velocity. This is better investigated with the change of variables

$$ \begin{aligned} \xi_1 &= \frac{X_1 + X_2}{\sqrt{2}}, & \xi_2 &= \frac{X_2 - X_1}{\sqrt{2}}, \\ \eta_1 &= \frac{Y_1 + Y_2}{\sqrt{2}}, & \eta_2 &= \frac{Y_2 - Y_1}{\sqrt{2}}, \end{aligned} \quad \text{(B1)} $$

which is canonical. Then, the Hamiltonian can be expressed as follows

$$ H_{\text{sc}}^{1,2} = H_0(p, \eta_1, x, \xi_1) + H_2(\eta_2, \xi_2) + \frac{\Delta_\omega}{2} (\xi_1\xi_2 + \eta_1\eta_2) \quad \text{(B2)} $$

with

$$ H_0(p, \eta_1, x, \xi_1) = \frac{p^2}{2} + \bar{\omega} \left( \frac{\xi_1^2 + \eta_1^2}{2} \right) + \sqrt{2} (\eta_1 \sin(x) - \xi_1 \cos(x)), \quad \text{(B3)} $$

and

$$ H_2(\eta_2, \xi_2) = \bar{\omega} \left( \frac{\xi_2^2 + \eta_2^2}{2} \right), \quad \text{(B4)} $$

for $\bar{\omega} := (\omega_1 + \omega_2)/2$.

One recognizes in $H_0$ the integrable model $H_{\text{sc}}^{1,1}$ (one wave, one particle) with frequency $\bar{\omega}$ and coupling $\sqrt{2}$, and in $H_2$ a harmonic oscillator (with an eigenfrequency matching that of the wave in $H_0$). When $\Delta_\omega = 0$, these systems are uncoupled, and the limit $\Delta_\omega \to 0$ can be investigated using KAM-type theory (noting that the Hamiltonian $H_2$ is degenerate as it has no shear).

Additionally, a Galileo transformation can be used to set $\bar{\omega} = 0$. Then, the limit $\Delta_\omega \to 0$ can also lead to slow chaos [19–21].

## Appendix C: Stability of locked solutions

By assuming the locked phases solution $\phi_i^* = n_i\pi$ for $n_i = 0, 1$ and $i = 1, 2$ and corresponding amplitudes $I_1^*, I_2^*$, the jacobian for equilibrium solutions reads

$$ \mathcal{J} = \begin{pmatrix} 0 & 0 & f_1 & 1 \\ 0 & 0 & 1 & f_2 \\ g_1 & 0 & 0 & 0 \\ 0 & g_2 & 0 & 0 \end{pmatrix} $$

where $f_i$ and $g_i$ are defined as

$$ f_i = f(I_i^*, n_i) := 1 + \frac{(-1)^{n_i}}{2\sqrt{2(I_i^*)^3}} \quad \text{and} $$

$$ g_i = g(I_i^*, n_i) := -(-1)^{n_i} \sqrt{2I_i}, \quad \text{for } i = 1, 2. $$

Given the biquadratic form of the jacobian's characteristic polynomial

$$ \lambda^4 - (g_1 f_1 + g_2 f_2) \lambda^2 + g_1 g_2 (f_1 f_2 - 1) = 0, $$

one may solve it analytically for its eigenvalues $\lambda_i$ as

$$ \lambda_{\pm, \pm} = \pm \sqrt{\Lambda_{\pm}}, $$

with

$$ \Lambda_{\pm} = \frac{1}{2} \left( g_1 f_1 + g_2 f_2 \pm \sqrt{(g_1 f_1 - g_2 f_2)^2 + 4g_1 g_2} \right), $$

and where all the four possible sign combinations are considered in the notation $\pm$.

## Appendix D: Decoupling limit

When analysing the system with weak coupling, that is in the limit $\beta = m \ll 1$, one recovers the paradigmatic 1.5-degree-of-freedom model of a particle slaved to two free waves. As done in appendix A, we set $k_1 = k_2 = k$ and $\epsilon\beta_1 = \epsilon\beta_2 = \beta$ in Hamiltonian (2), for the case $N = 1, M = 2$. Setting $p' = p/m$ in the equations of motion then yields

$$ \ddot{x} = -\sqrt{2I_1} \sin(x - \theta_1) - \sqrt{2I_2} \sin(x - \theta_2) $$

$$ \dot{I}_i = \beta \sqrt{2I_i} \sin(x - \theta_i) $$

$$ \dot{\theta}_i - \omega_i = \beta(2I_i)^{-1/2} \cos(x - \theta_i), $$

so that the particle evolves with the waves whereas the waves themselves evolve under an $O(\beta)$ feedback, becoming free waves in the limit $\beta \to 0$. This scale separation does not satisfy the standard KAM hypotheses, as the waves dynamics are degenerate (harmonic oscillators).

An alternative way to apply KAM theory to the current model would be to consider Hamiltonian (5) in the form

$$ H' = H_1(I'_1, \phi_1) + H_2(I'_2, \phi_2) + I'_1 I'_2, $$

where

$$ H_1 = I_1'^2/2 - \bar{P}I_1' - \sqrt{2I_1'} \cos \phi_1 \quad \text{(D1)} $$

and

$$ H_2 = I_2'^2/2 - (\bar{P} - \Delta_\omega)I_2' - \sqrt{2I_2'} \cos \phi_2, \quad \text{(D2)} $$

where both $H_1$ and $H_2$ reduce by a mere Galileo transformation to the $N = M = 1$ integrable model, (eq. (9) in Gomes et al [12]). We leave these analyses for future work.
[1] Y. Elskens and D. Escande. *Microscopic Dynamics of Plasmas and Chaos*. IOP Publishing, Bristol, 2003.

[2] D. F. Escande and Y. Elskens. Microscopic dynamics of plasmas and chaos: the wave-particle interaction paradigm. *Plasma Phys. Control Fusion*, 45:A115–A124, 2003.

[3] R. Balescu. *Statistical Mechanics of Charged Particles*. Wiley-Interscience, London, 1963.

[4] F. Doveil; A. Macor and Kh. Auhmani. Wave-particle interaction investigated in a travelling wave tube. *Plasma Phys. Control Fusion*, 47:A261–A271, 2005.

[5] M. C. de Sousa; F. Doveil; Y. Elskens and I. L. Caldas. Wave-particle interactions in a long traveling wave tube with upgraded helix. *Physics of Plasmas*, 27(9):093108, 2020.

[6] A. Macor; F. Doveil and Y. Elskens. Electron climbing a “devil’s staircase” in wave-particle interaction. *Physical Review Letters*, 301(P-2.056), 2006.

[7] E. Ott. *Chaos in Dynamical Systems*. Cambridge University Press, New York, 1993.

[8] G. R. Smith and N. R. Pereira. Phase-locked particle motion in a large-amplitude plasma wave. *Phys. Fluids*, 21:2253–2262, 1978.

[9] T. M. O'Neil; J. H. Winfrey and J. H. Malmberg. Non-linear interaction of a small cold beam and a plasma. *The Physics of Fluids*, 14(6):1204–1212, 1971.

[10] I. N. Onishchenko; A. R. Linetskii; N. G. Matsiborko; V. D. Shapiro and V. I. Shevchenko. Contribution to the nonlinear theory of excitation of a monochromatic plasma wave by an electron beam. *JETP Letters*, 12(8):281–285, 1970.

[11] H. E. Mynick and A. K. Kaufman. Soluble theory of non-linear beam-plasma interaction. *Phys. Fluids*, 21:652–663, 1978.

[12] J. V. Gomes; M. C. de Sousa; R. L. Viana; I. L. Caldas and Y. Elskens. Low-dimensional chaos in the single

wave model for self-consistent wave-particle Hamiltonian. *Chaos*, 31, 2021.

[13] J. C. Adam; G. Laval and I. Mendonça. Time-dependent nonlinear Langmuir waves. *Phys. Fluids*, 24:260–267, 1981.

[14] D. del Castillo-Negrete. Dynamics and self-consistent chaos in a mean field hamiltonian model. *Dynamics and Thermodynamics of Systems with Long-Range Interactions*, Les Houches, 18–22 February 2002, edited by T. Dauxois, S. Ruffo, E. Arimondo, and M. Wilkens:407–436.

[15] B. V. Chirikov. A universal instability of many-dimensional oscillator systems. *Physics Reports*, 52(5):263–379, 1979.

[16] D. F. Escande. Stochasticity in classical hamiltonian systems: Universal aspects. *Physics Reports*, 121(3-4):165–261, 1985.

[17] C. Skokos; T. Bountis; C. G. Antonopoulos and M. N. Vrahatis. Detecting order and chaos in hamiltonian systems by the SALI method. *Journal of Physics A Mathematical and General*, 37:6269–6284, 2004.

[18] G. A. Gottwald C. H. Skokos and J. Laskar. *Chaos Detection and Predictability*. Springer-Verlag, Berlin Heidelberg, 2015.

[19] Y. Elskens and D. F. Escande. Slowly pulsating separatrices sweep homoclinic tangles where islands must be small – an extension of classical adiabatic theory. *Nonlinearity*, 4:615–667, 1991.

[20] V. V. Sidorenko A. I. Neishtadt and D. V. Treschev. Stable periodic motions in the problem on passage through a separatrix. *Chaos*, 7(2):1–11, 1997.

[21] A. I. Neishtadt; C. Simó; D. Treschev and A. Vasiliev. Periodic orbits and stability islands in chaotic seas created by separatrix crossings in slow-fast systems. *Discrete and continuous dynamical systems series B*, 10(2&3):621–650, 2008.