Quasiperiodic oscillations and homoclinic orbits in the nonlinear nonlocal Schr\"odinger equation

Quasiperiodic oscillations and shape-transformations of higher-order bright solitons in nonlinear nonlocal media have been frequently observed in recent years, however, the origin of these phenomena was never completely elucidated. In this paper, we perform a linear stability analysis of these higher-order solitons by solving the Bogoliubov-de Gennes equations. This enables us to understand the emergence of a new oscillatory state as a growing unstable mode of a higher-order soliton. Using dynamically important states as a basis, we provide low-dimensional visualizations of the dynamics and identify quasiperiodic and homoclinic orbits, linking the latter to shape-transformations.


Introduction
Bright solitons are particle-like nonlinear localized waves , that keep their form while evolving due to a compensation of diffraction or dispersion of the medium by the nonlinear self-induced modification of the medium [1]. Usually, solitons are studied in systems exhibiting local nonlinearities, where the guiding properties of the medium at a particular point in space depend solely on the wave intensity at that particular point [2]. Here, we consider nonlocal nonlinearities, i.e. situations in which the nonlinear response of the medium at a point depends on the wave intensity in a certain neighborhood of that point, where the extent of this neighborhood is referred to as degree of nonlocality. Nonlocal nonlinearities are ubiquitous in nature, for example, when the nonlinearity is associated with some sort of transport process, such as heat conduction in media with thermal response [3][4][5], diffusion of charge carriers [6,7] or atoms/molecules in atomic vapors [8,9]. Nonlinearities are also nonlocal in case of long-range interaction of atoms in Bose-Einstein condensates (BEC), such as in case of dipolar BEC [10][11][12][13] or BEC with Rydberg-mediated [14,15] interactions. In addition, long-range interactions of molecules in nematic liquid crystals also result in nonlocal nonlinearities [16][17][18][19].
The balance between diffraction and nonlinearity may lead to stable solitons withstanding even strong perturbations. In particular, it has been shown, that nonlocal nonlinearities crucially modify stability properties of localized waves. With respect to bright solitons, they lead to a much more robust evolution as compared to its local counterpart [20,21]. This is due to the fact, that nonlocality acts like a filter by averaging or smoothing-out effect on perturbations which would otherwise grow in case of local response of the medium [22]. For example, higher-dimensional solitons would collapse for systems exhibiting local nonlinearities, whereas they can be stabilized by nonlocality [23][24][25].
In this work, we investigate the linear stability and nonlinear dynamics of higherorder solitons. In particular, we study the quadrupole soliton Q and the second-order radial soliton R 2 (a hump with a ring), as sketched in Fig. 1. For those solutions, a quasiperiodic shape transformation between states of different symmetries has been observed recently in [26,27]. However, a complete understanding of this spectacular phenomenon is still missing. One difficulty in the analysis of the shape transformations is that they cannot be described solely in terms of linear perturbation because they are not small [27]. Nevertheless, here we show that in spite of the fact that we are dealing with a highly nonlinear phenomenon, deeper understanding can be gained from the linear stability analysis of the corresponding Bogoliubov-de Gennes (BdG) equations. In other words, solutions of the linear stability analysis of the solitons are used to describe wave dynamics in the neighborhood of a soliton solution. Moreover, in order to fully understand nonlinear dynamics, we employ and further develop techniques recently introduced in dynamical systems studies of dissipative partial differential equations (PDE) [28,29]. These methods employ projection of PDE solutions from a functional infinite space onto a finite number of important physical states or dynamically relevant directions. Here, the relevant directions are mainly the unstable and stable internal modes of the solitons. The introduction of these low-dimensional projections will allow us to interpret the non-periodic soliton oscillations as indication of homoclinic connections. Moreover, we are able to understand how different solutions, including quasiperiodic oscillations, are organized by this homoclinic connection. The same analysis should also work for a larger variety of higher-order solitons of this nonlocal system, such as those presented e.g. in [26,27].
The paper is organized as follows: In Sec. 2, we introduce the governing equations of motion. In Sec. 3, we solve the BdG equation to find the internal modes of the quadrupole soliton Q as well as the second-order radial soliton R 2 . In Sec. 4, we discuss nonlinear soliton propagation, introduce low-dimensional projections and study homoclinic and quasiperiodic trajectories in this representation. Finally, we will conclude in Sec. 5.

Model equations
The underlying model equation for our subsequent considerations is the nonlocal nonlinear Schödinger equation (NLS) where ∆ = ∂ xx + ∂ yy denotes the transverse Laplacian. Depending on the actual context, |ψ(r, t)| 2 can be identified with either the intensity of an optical beam in scalar, paraxial approximation, or the density of a two-dimensional BEC within mean field approximation. The nonlinearity θ is given by the convolution integral where the kernel K is determined by the physical system under investigation, and r = (x, y). If K(r) = K(|r|), then Eq. (1) is invariant under rotation and the angular momentum is conserved. This is the case here, as we consider the Gaussian nonlocal model, for which quasiperiodic oscillations have been originally observed [26,27]: Even though there is no actual physical system associated with the Gaussian model, it is commonly used in the literature as a toy model for nonlocal nonlinearities. Note that without loss of generality the width of the kernel K has been set to unity, in order to have the same scaling as used in [26,27].

Linear stability analysis of higher-order solitons
Let Φ be a bright solitonic solution to our governing equation (1) where λ is the propagation constant or chemical potential for the case of optical beam or BEC, respectively, and ψ S denotes the stationary profile of the soliton. Because we will not consider solitons carrying angular momenta (e.g., vortices), we can choose ψ S (r) to be real. In order to find numerically exact stationary profiles ψ S (r), we use variational solutions as input to an iterative solver [30]. Typically, we use a grid of 400 × 400 points to determine ψ S (r). This transverse resolution is also employed for numerical integration of Eq. (1), i.e., for beam propagation or time evolution of the two-dimensional BEC. Figure 2 shows solitonic family curves or the two higher order solitons we choose to study here, the second-order radial state R 2 and the quadrupole Q. Apart from the total angular momentum, there are two conserved functionals, i.e. the Hamiltonian H[ψ] associated with invariance with respect to time-translations and the mass M[ψ] due to a global U(1) phase-invariance: Obviously, the family curves for the R 2 and Q solitons are quite close to each other, which was used in [26] to explain the observed quasiperiodic shape transformations (energy crossing). However, we will see in the following analysis of projected propagation dynamics in Sec. 4 that this very intuitive picture does not hold. Let us first recall that the linear stability of solitonic solutions can be studied as an eigenvalue problem as follows. We introduce a small perturbation δψ(r, t) to our solitonic solution ψ S (r) via ψ(r, t) = [ψ S (r) + δψ(r, t)] e iλt Plugging Eq. (7) into the governing equation Eq. (1) and retaining only first order terms in δψ, yields the following (linear) evolution equation for δψ: With the ansatz for the perturbation we can derive the eigenvalue problem (BdG equation) Real-valued eigenvalues κ of Eq. (10) are termed orbitally stable and the corresponding eigenvector (δu, δv) can be chosen real-valued. On the other hand, complex eigenvalues with negative imaginary part indicate exponentially growing instabilities. We note that due to the special structure of Eq. (10) [which has its origins in the Hamiltonian structure of Eq. (1)], if κ is an eigenvalue, then −κ as well as ±κ * are also eigenvalues. Next, we solve Eq. (10) in order to obtain the internal modes of the second-order radial soliton R 2 and the quadrupole Q, respectively. A trivial solution to this problem is always given by (δu, δv) = ±(ψ S , −ψ S ) with eigenvalue κ = 0. This so-called trivial phase mode is linked to the phase invariance of solitons. Derivatives of this trivial phase mode with respect to x or y are also trivial eigenvectors ‡ with eigenvalue κ = 0, and thus the eigenvalue κ = 0 is degenerate. Moreover, due to symmetry properties of the system trivial modes appear twice in the spectrum, i.e., we expect sixfold degeneracy of the eigenvalue κ = 0. However, when solving the discretized version of Eq. (10) numerically, this degeneracy may be lifted. Thus, degenerate eigenvectors with zero eigenvalues may in fact become slightly complex without actually indicating an instability. In other words, their nonzero imaginary part is a numerical artefact of the discretization and occurs because the full eigenspace has to be spanned by the eigenvectors. The actual computation of the linear eigenvalue problem Eq. (10) is numerically expensive, since the matrix we have to diagonalize is full, i.e. all entries are nonzero. In order to achieve reasonable computation times, we usually reduce the grid-size to 100 × 100 points only. Then, the matrix we have to diagonalize has 4 × 10 8 non-zero elements.
In Fig. 3, we show the spectrum of the linear stability analysis (BdG equation) for the second-order radial soliton R 2 and the quadrupole soliton Q [a) resp. b)] for mass M = 200. Note that for modes with purely imaginary eigenvalue κ = iIm κ, Eq. (9) reads δψ(r, t) = [δu(r) + δv * (r)] e −Im (κ)t , and it makes sense to definê Only the second order radial state R 2 is unstable, and we name the two unstable internal modesê 1 ,ê 2 . The unstable modesê 1 ,ê 2 ought to be degenerate for symmetry reasons, the small splitting of the eigenvalues (κ 1 ≈ −2.7i, κ 2 ≈ −2.5i) is again a numerical artefact due to the discretization of the eigenvalue problem Eq. (10). Interestingly, the shape of the unstable eigenmodesê 1 (r),ê 2 (r) resembles quadrupoles. In fact, for practical purposes (see next section) as well as to verify these findings we furthermore solved the eigenvalue problem Eq. (10) for R 2 on a radial grid [31] with eightfold resolution. Then, instead of two stable and unstable quadrupoles, one finds one stable and unstable vortex with topological charge m = ±2 and |κ| ≈ 2.74. The vortices corresponding to m = 2 and m = −2 can be superposed to again yield the quadrupolesê 1 ,ê 2 found already with the full 2D solver, but with much higher precision. Because Eq. (10) is linear, the amplitudes of theê j are not fixed, and we normalize the latter according to The quadrupole soliton Q in Fig. 3b) is stable, because all complex eigenvalues correspond to trivial modes and hence the complex form of these eigenvalues is a numerical artefact as discussed above. However, the quadrupole becomes linearly unstable for M 90, as indicated in Fig. 2 by dashed lines. In Fig. 4, we show the results of our numerical stablility analysis for the quadrupole soliton Q with mass M = 85. Interestingly, the unstable modeê 1 with κ 1 ≈ −1.2i resembles the secondorder radial soliton R 2 , i.e., a hump with a (modulated, i.e. not rotationally symmetric) ring.

Projected nonlinear dynamics
The typical dynamics for R 2 (here for M = 200) as an initial condition is shown in Fig. 5 a). To determine the shape of R 2 , we use the iterative solver mentioned above on a grid containing 400 × 400 points, and we use the same grid for the actual propagation.  The unstable eigenmodeê 1 (r) resembles the shape of R 2 (see inset), but is of course not rotationally symmetric. Again, the degeneracy of the trivial modes is lifted, a numerical artefact due to the discretization of the eigenvalue problem Eq. (10). For sake of readability, the inset shows the absolute square |ê(r)| 2 = |δu(r) + δv * (r)| 2 only.
As we have seen in Sec. 3, the second-order radial soliton R 2 is unstable over the whole range of mass M and therefore any perturbation, that has a non-zero overlap with the unstable internal modesê 1 ,ê 2 will lead to an exponential growth of the latter. Practically, the residual in numerical determination of R 2 as well as the propagation algorithm based on the Fourier split-step method [1] lead to inevitable numerical noise when propagating and therefore trigger the instability without adding any additional perturbation. In our case, however, we added the eigenmodeê 1 as initial perturbation with tiny amplitude ∼ 10 −4 to the soliton R 2 to control the breakup in a preferred direction. For small times the dynamics is governed by the exponential growth of the unstable internal modeê 1 , while for later times the evolution becomes highly non-linear, exhibiting oscillations between R 2 [see inset (α) in Fig. 5 a)] and a state that resembles the quadrupole soliton Q [see inset (β) in Fig. 5 a)] [26]. This state (β) we will call the Figure 5. a) Evolution of the peak-intensity of the second-order radial soliton R 2 with mass M = 200 (upper blue curve). As expected from the stability analysis Fig. 3 b), the peak-intensity of the quadrupole soliton Q with same mass (lower red curve) is constant during propagation.  (19)]. In this three-dimensional projection, the distance between the quadrupole Q (γ) and the "turning point" (β) becomes apparent. For reasons of clarity, the 3Ddynamics (blue) is additionally projected into (S, w)-plane (black), and the orbit of the the quadrupole is again shown in red. The three insets show snapshots of the nonlinear dynamics.
"turning point". In the following we will examine in detail the origin and properties of these oscillations.

Projection methods
Let us now introduce the projection method mentioned in the introduction [28,29] and adopt it to our problem. To this end, we recall the scalar product of two complex functions f and g, defined as Obviously, the (unstable) internal modesê j of R 2 introduced before [see Fig. 3] are not orthogonal to their complex conjugate (stable)ê * j (j = 1, 2) counterparts with respect to this inner product. In other words, stable and unstable eigenspaces E s and E u spanned by eigenfunctions {ê * 1 ,ê * 2 } resp. {ê 1 ,ê 2 } are not mutually orthogonal. Thus, straightforward projections ontoê j andê * j do not help to elucidate the propagation dynamics. To overcome this difficulty we introduce a set of functions which is biorthogonal toê j ,ê * j using a Gram-Schmidt-like technique as follows. First, we define e j⊥ =ê j − ê * j ,ê j ê * j , which is simply the projection of the unstable eigenmodeê j onto the orthogonal complement of the stable eigenmodeê * j . Second, we note that (e j⊥ ) * =ê * j − ê j ,ê * j ê j corresponds to projection of the stable eigenmodeê * j onto the orthogonal complement Figure 6. Schematic sketch of the relation betweenê j ,ê * j , e j⊥ , and (e j⊥ ) * . By construction, e j⊥ is orthogonal to the stable eigenvectorê * j , and (e j⊥ ) * is orthogonal to the unstable eigenvectorê j . It is worth to notice that e j⊥ and (e j⊥ ) * are not orthogonal to each other. b) and c) show the modulus squared of the internal modê e 1 and e 1⊥ , respectively. of the unstable eigenmodeê j . Then, it is easy to verify biorthogonality of e j⊥ , (e j⊥ ) * with respect toê j ,ê * j : ê j , (e j⊥ ) * = 0 ê j , e j⊥ = 0 (16) In Fig. 6a) a schematic sketch of the relation betweenê j ,ê * j , e j⊥ , and (e j⊥ ) * is depicted, and Fig. 6b-c) showê 1 and e q⊥ explicitly. It is worth to notice that e j⊥ , (e j⊥ ) * = 0, i.e., e j⊥ and (e j⊥ ) * are not orthogonal to each other.
In order to analyze the propagation dynamics of a solution ψ(x, t) of Eq. (1), we introduce the quantities U j = e j⊥ , ψ , S j = (e j⊥ ) * , ψ .
By construction, U j is associated with the unstable eigenmode only (e j⊥ is orthogonal to the stable one), while S j is associated with the stable eigenmode only. Finally, for R 2 , the two unstable eigenvectorsê 1 ,ê 2 are degenerate (due to rotational symmetry about the origin), therefore we introduce the rotationally invariant projected variables Then, any pair of wavefunctions ψ 1 (x, t) and ψ 2 (x, t) related through a rotation amounts to the same value of U(t) and S(t). For a rigorous proof, see Appendix Appendix A.  Fig. 5 a) in the variables S(t), U(t) introduced in Eq. (18). We clearly see the second-order radial soliton R 2 (α) decaying into a quadrupole-like state (β), the "turning point", and then coming back to R 2 . In the vicinity of R 2 , the decay starts via the local unstable eigenspace E u (i.e., U(t) > 0, S(t) ≈ 0), and the revival of R 2 happens via the local stable eigenspace E s (i.e., S(t) > 0, U(t) ≈ 0). The fact that the system repeatedly returns (close) to its initial state R 2 and remains at this point some finite, non-constant time with (nearly) zero velocity, hints at the existence of a homoclinic connection. A homoclinic connection is a solution which is asymptotic to R 2 both in the t → ∞ and t → −∞ limit. The time-span, in which the solution remains close to its initial state R 2 , i.e. the homoclinic point (α) in Fig. 5 b), with practically zero velocity, corresponds to intervals with maximum (nearly) constant peak-intensities in Fig. 5 a). Because we added a small perturbation in the direction of the eigenmodeê 1 to the initial condition R 2 , and the presence of numerical noise in general, we do not see the exact homoclinic connection in our numerical simulations; as the trajectory comes back towards R 2 along E s , there is always a small perturbation along the unstable eigenspace E u and the trajectory leaves the neighborhood of R 2 to return to it later on. We want to stress here that the existence of homoclinic connections is by no means anticipated in general; our numerical results however indicate the existence of such homoclinic connections and their persistence along a large range of the mass M.

Indication of homoclinic connections
To further illustrate that the "turning point" (β) is indeed well-separated from the quadrupole soliton Q (γ), we introduce a third variable w by projecting the solitonic wave function ψ onto the radial soliton R 2 , Obviously, for ψ = R 2 we find w = 1, while for ψ = Q for symmetry reasons we have w = 0. Figure 5 c) shows the resulting projected dynamics on the variables U, S, w. We clearly recognize similarities with Fig. 5 b), however, it becomes much more clear how the solution evolves from its origin (α) and becomes much more "quadrupole-like" in (β). In particular, the important separation between the quadrupole-like "turningpoint" (β), which still maintains a nonzero projection on R 2 and the quadrupole soliton Q (γ) becomes evident.

Quasiperiodic motion
In the previous section, we argued that due to numerical limitations, we cannot actually track the homoclinic orbit precisely, but what we find are trajectories that are very close to the homoclinic connection. In the present section we will further probe the dynamical importance of the homoclinic orbit by studying trajectories adjacent to it. In a sense, the "turning point" (β) of the homoclinic orbit is a state "in between" R 2 and Q. Here we will investigate the dynamics of such "in between" states obtained by perturbing the homoclinic orbit at the "turning point" (β). The perturbations we will consider are not necessarily small and, as we will see, they typically lead to quasiperiodic oscillations. A homoclinic orbit is obtained by (slightly) perturbing the initial wavefunction of R 2 in the direction of one of the unstable modes (e.g. ofê 1 ) and integrating Eq. (1) forward in time. Choosing the direction of the initial perturbation fixes the "orientation" of the subsequent dynamics, and we can thus decompose the wavefunction at the turning point [point (β) in Fig. 5] t t into a part parallel to the quadrupole soliton Q and a remainder L where c Q = Q, ψ(r, t t ) / Q, Q was introduced §. Perturbed wavefunctions ψ Γ (r) are then constructed through where Γ parametrizes mixed states between R 2 and Q, and, in Eq. (22), the wavefunction was normalized. Clearly, for Γ = 1, the homoclinic trajectory of R 2 can be recovered, whereas of Γ = 0, the quadrupole soliton is recovered. In the following the time evolution of the function ψ Γ will be studied. Let us first consider the dynamics for Γ = 1.01 as shown in Fig. 7, which indicates quasiperiodic behavior for small times (up to t ≃ 25). The time spend by this orbit close to R 2 is much smaller than for the homoclinic connection, of the previous section. This becomes apparent when comparing the peak-intensity evolution in Fig. 7a) with the one in Fig. 5a). In the (U, S, w) projection, this fact results in a smoother curve close to the origin (whereas for a homoclinic connection a kink appears as R 2 is approached, while the "velocity" approaches zero). On the other hand, in the intensity representation, Fig. 7c-h), the difference between homoclinic and quasiperiodic behavior is much harder to discern. Propagation in Fig. 7a) and b) is shown until t = 35, when the dynamics already deviates from the quasiperiodic orbit, indicating that the latter is unstable. This behavior hints to the existence of some chaotic region in state-space, an issue that will be studied elsewhere.
On the other hand, the dynamics for Γ = 0.99, shown in Fig. 8, appear again quasiperiodic (see also the discussion of the Fourier spectra in Sec. 4.4), but in this case the orbit appears stable, as it persists at least up to t = 1500. The qualitatively different behaviour for Γ = 1.01 and Γ = 0.99 with respect to stability further corroborates the importance of the homoclinic solution Γ = 1.00 (R 2 ). In a certain sense, the homoclinic orbit "organizes" regions of stability in parameter space. However, the homoclinic orbit should not be seen as a kind of "boundary" between regions of different stability behaviours, because it is just a one-dimensional line in the highly-dimensional parameter space.
Let us finally consider the trajectory in Fig. 9 which is far away from both the quadrupole soliton as well as from the "turning point" (β) by letting Γ = 0.5. The dynamics is still quasiperiodic and stable (at least up to t = 1500), but involves multiple frequencies. Interestingly, the dominant frequency of oscillation with period T ≈ 2.6 can be related to a stable eigenvalue of the quadrupole soliton Q for M = 200. In the (stable) eigenvalue spectrum of Q shown in Fig. 3b), the internal mode with κ ≈ 2.6 resembles a (modulated) ring with a hump (not shown). The duration of one period T would then be given by T = 2π/κ ≈ 2.4, which is what we find when we slightly § More generally, if the direction of the breakup is arbitrary, one may generalize Eq. (20) by decomposing ψ(t t ) into two quadrupoles Q 1 , Q 2 , where Q 1 is rotated by π/4 with respect to Q 2 , via ψ(r, t t ) = c Q1 Q 1 (r) + c Q2 Q 2 (r) + L(r).   Fig. 5, where the blue curve again represents the actual 3D dynamics and the black curve its projection on the (s, w)plane, and the red curve is the orbit of the quadrupole.
perturb the quadrupole soliton Q by this mode. Moreover, for Γ = 0.1 (not shown) we also find an oscillation with period T ≈ 2.4. In both case, the propagation dynamics resemble the one shown in Fig. 9 for Γ = 0.5. Thus, even though for Γ = 0.5 we are no longer in the region where perturbation analysis of the quadrupole soliton Q holds, we still find qualitatively similar dynamics. We note that in the same system Eq. (1), quasiperiodic nonlinear solutions (so-called azimuthons) linked to stable internal modes of solitons were reported earlier [31,32].
To sum up, we have identified a family of stable quasiperiodic solutions to Eq. (1), starting from ψ Γ given in Eq. (21) and 0 < Γ < 1. The two limiting solutions are the stable quadrupole solitons Q (Γ = 0) and the homoclinic orbit linked to the unstable radial solitons R 2 (Γ = 1). We want to emphasize here that for lower masses, where the quadrupole soliton Q becomes unstable (e.g., M = 85), we were not able to find stable quasiperiodic solutions by the same construction.

Fourier spectrum
Further insight can be gained by considering the Fourier spectrum of the above trajectories. Given a trajectory ψ(r, t) we compute the modulus of the Fourier transform F of the wavefunction at a fixed point in space (in our case the origin r = 0): For a bright soliton solution of the form Eq. (4), one would expect f (ω) to comprise of a single sharp peak at ω = λ. On the other hand, in the case of quasiperiodic dynamics with vibration frequency Ω and propagation constant λ, one would expect peaks at λ + m Ω, where m is integer. This is readily verified for the orbits with a = 0.99 and a = 0.5, as can be seen in Fig. 10, where we see sharp peaks associated with these orbits. On the other hand, there is no well defined periodicity associated with the homoclinic orbit, since the time spent in the vicinity of R 2 is in principle infinite. In practice, this time is greatly affected by numerical noise and the spectrum appears continuous [see Fig. 10a)]. Even if it is possible to associate a dominant frequency Ω with the homoclinic orbit, f (ω) around Ω is much broader than in the case of quasiperiodic orbits for Γ = 0.99 and Γ = 0.5 [see Fig. 10b) and c)] .
A limitation on the spectral resolution for f (ω) for the homoclinic orbit appears due to the fact that dynamics become unstable around t = 520. Here, we used the interval t = [0 : 500] to compute Thus, the Fourier spectra yield an additional indication of the qualitatively different nature of the dynamics of Sec. 4.2 from the quasiperiodic motion of Sec. 4.3, providing further support for the conjectured existence of an underlying homoclinic connection in the former case.

Conclusions
In previous works, an oscillatory shape-transformation of modes in nonlocal media has been observed [26,27]. In this paper, we approached this phenomenon by means of linear stability analysis and projection techniques borrowed from dynamical systems studies of dissipative PDEs. By studying the linear stability of the quadrupole soliton Q and the second-order radial soliton R 2 , we found that the former becomes linearly stable for mass M 90, whereas the latter remains linearly unstable for all masses. The initial stage of the shape-transformations under consideration, i.e. the emergence of a new state on top of R 2 , can be understood in terms of this linear instability, which is triggered by the unavoidable numerical noise. However, the most striking feature of the dynamics, i.e. the return to the initial state, is inherently nonlinear, as it occurs only after the linear instability saturates. To study this phenomenon, we introduced a low-dimensional representation of the dynamics, through a projection to dynamically important states, which were constructed from the radial soliton R 2 itself and its unstable/stable eigenmodes. Projecting the time evolution of the wavefunction ψ(r, t) (obtained by integrating the NLS) onto these states allows a visualization of oscillatory shape-transformations in terms of trajectories, revealing that shape-transformations can be interpreted as a homoclinic orbit leaving and re-approaching R 2 . Moreover, in the neighborhood of this homoclinic orbit we found quasiperiodic solutions, which for small enough perturbations resemble the homoclinic connection. This indicates that the homoclinic connection provides a basic recurrence mechanism around which quasiperiodic dynamics is organized, as is common in lower-dimensional dynamical systems [33]. We were also able to construct and identify a whole family of stable quasiperiodic orbits when the quadrupole soliton Q is stable.
The projection method introduced here allows a compact representation of the dynamics, dual to the commonly used intensity plots. Moreover, in certain cases it helps to uncover features of the dynamics that are not apparent in snapshots of the intensity evolution. We expect that similar studies can be carried out for other states exhibiting similar dynamics [26] and that our projection method (or similar extensions of the methods of Refs. [28,29]) could be applied to a variety of high-and infinitedimensional conservative systems. the spectrum. Thus, compared to the other two spectra shown in Fig. 10, where the propagation was performed until t = 1500, the spectral resolution is coarser by a factor of three. Comparing with Eq. (A.2) we conclude that α(θ) = α * (θ) and thus our representation is real, and that α(θ − θ 0 ) = −α(θ + θ 0 ). ¶ Using Eqs. = | e 1⊥ , ψ | 2 + | e 2⊥ , ψ | 2 = U 2 (t) .