Quasiperiodic oscillations and homoclinic orbits in the nonlinear nonlocal Schrödinger equation

Bright solitons are particle-like nonlinear localized waves, that keep their form while evolving due to a compensation of diffraction or dispersion 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–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–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–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 their local counterparts [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–25].

In this work, we investigate the linear stability and nonlinear dynamics of higher-order solitons. In particular, we study the quadrupole soliton Q and the second-order radial soliton R₂ (a hump with a ring), as sketched in figure 1. For those solutions, a quasiperiodic shape transformation between states of different symmetries has been observed numerically in [26, 27]. However, a complete understanding of this spectacular phenomenon, which is beyond conventional instability of higher-order solitons in local nonlinear media due to symmetry-breaking and modulational instability [28], is still missing. One difficulty in the analysis of the shape transformations is that they cannot be described solely in terms of linear perturbation analysis since perturbations 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 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) [29, 30]. These methods project 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 nonperiodic 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].

Zoom In Zoom Out Reset image size

The paper is organized as follows. In section 2, we introduce the governing equations of motion. In section 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₂. In section 4, we discuss nonlinear soliton propagation, introduce low-dimensional projections and in section 5 we study homoclinic and quasiperiodic trajectories in this representation. Finally, we will conclude in section 6.

The underlying model equation for our subsequent considerations is the nonlocal nonlinear Schrödinger equation

$$\begin{equation} \mathrm{i}\partial_t\psi + \Delta \psi + \theta \psi = 0, \end{equation} \tag{ 1 }$$

where Δ = ∂_xx + ∂_yy denotes the transverse Laplacian. Depending on the actual context, |ψ(r,t)|² can be identified with either the intensity of an optical beam in scalar, paraxial approximation or the density of a two-dimensional (2D) BEC within mean field approximation. The nonlinearity θ is given by the convolution integral

$$\begin{equation} \theta = \int K(\mathbf{r}-\mathbf{r}^\prime) |\psi(\mathbf{r}^\prime,t)|^2 \,\mathrm{d}^2 \mathbf{r}^\prime, \end{equation} \tag{ 2 }$$

where the kernel K is determined by the physical system under investigation, and r = (x,y). If K(r) = K(|r|), then equation (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]:

$$\begin{equation} K(\mathbf{r})=\mathrm{e}^{-\mathbf{r}^2}. \end{equation} \tag{ 3 }$$

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

Let Φ be a bright solitonic solution to our governing equation (1)

$$\begin{equation} \Phi(\mathbf{r},t)=\psi_{\mathrm{S}}(\mathbf{r}) \,\mathrm{e}^{\mathrm{i}\lambda t}, \end{equation} \tag{ 4 }$$

where λ is the propagation constant or chemical potential for the case of optical beams 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 [31]. Typically, we use a grid of 400 × 400 points to determine ψ_S(r). This transverse resolution is also employed for numerical integration of equation (1), i.e. for beam propagation or time evolution of the 2D BEC. Figure 2 shows solitonic family curves of the two higher order solitons we choose to study here, the second-order radial state R₂ and the quadrupole Q. Apart from the total angular momentum, there are other conserved functionals, namely the Hamiltonian $\mathcal {H}[\psi ]$ associated with invariance with respect to time-translations and the mass M[ψ] due to a global U(1) phase-invariance:

$$\begin{equation} \mathcal{H}[\psi_{\mathrm{S}}] = \int |\nabla \psi_{\mathrm{S}}|^2\,\mathrm{d}^2\mathbf{r} - \frac{1}{2}\int |\psi_{\mathrm{S}}(\mathbf{r})|^2 K(\mathbf{r}-\mathbf{r}^\prime) |\psi_{\mathrm{S}}(\mathbf{r}^\prime)|^2 \,\mathrm{d}^2\mathbf{r}^\prime \,\mathrm{d}^2\mathbf{r}, \end{equation} \tag{ 5 }$$

$$\begin{equation} M[\psi_{\mathrm{S}}] = \int |\psi_{\mathrm{S}}|^2\, \mathrm{d}^2\mathbf{r}. \end{equation} \tag{ 6 }$$

Obviously, the family curves for the R₂ and Q solitons are quite close to each other. It was already suggested in [26] that such 'energy crossing' is a necessary condition for shape transformations to occur. However, energy crossing is not sufficient for the latter. Moreover, we will see in the following analysis of projected propagation dynamics in section 4 that the solitons do not fully convert into each other.

Zoom In Zoom Out Reset image size

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

$$\begin{equation} \psi(\mathbf{r},t) = \left[\psi_{\mathrm{S}}(\mathbf{r})+\delta\psi(\mathbf{r},t)\right]\,\mathrm{e}^{\mathrm{i}\lambda t}. \end{equation} \tag{ 7 }$$

Plugging equation (7) into the governing equation (1) and retaining only first order terms in δψ, yields the following (linear) evolution equation for δψ:

$$\begin{eqnarray}&&\fl \left[ \mathrm{i }\partial_t -\lambda + \Delta+ \int K(|\mathbf{r}-\mathbf{r}^{\prime}|) \psi_{\mathrm{S}}^2(\mathbf{r}^{\prime})\mathrm{d}\mathbf{r}^{\prime} \right]\delta \psi(\mathbf{r},t) \nonumber\\ + \psi_{\mathrm{S}}(\mathbf{r}) \int K(|\mathbf{r}-\mathbf{r}^{\prime}|) \psi_{\mathrm{S}}(\mathbf{r}^{\prime}) \left[ \delta \psi(\mathbf{r}^{\prime},t) + \delta \psi^*(\mathbf{r}^{\prime},t) \right] \mathrm{d}^2\mathbf{r}^{\prime} &=0. \end{eqnarray} \tag{ 8 }$$

With the ansatz

$$\begin{equation} \delta \psi(\mathbf{r},t) = \delta u(\mathbf{r})\mathrm{e}^{\mathrm{i}\kappa t} + \delta v^*(\mathbf{r})\,\mathrm{e}^{-\mathrm{i}\kappa^* t} \end{equation} \tag{ 9 }$$

for the perturbation we can derive the eigenvalue problem (BdG equation)

$$\begin{eqnarray}&&\fl \left[\Delta -\lambda + \int K(|\mathbf{r}-\mathbf{r}^{\prime}|) \psi_{\mathrm{S}}^2(\mathbf{r}^{\prime})\,\mathrm{d}^2\mathbf{r}^{\prime} \right] \delta u (\mathbf{r}) \nonumber\\ &&+\psi_{\mathrm{S}}(\mathbf{r})\int K(|\mathbf{r}-\mathbf{r}^{\prime}|)\psi_{\mathrm{S}}(\mathbf{r}^{\prime}) \left[\delta u(\mathbf{r}^{\prime}) + \delta v(\mathbf{r}^{\prime}) \right] \,\mathrm{d}^2\mathbf{r}^{\prime} = \kappa \delta u(\mathbf{r}), \end{eqnarray} \tag{ 10a }$$

$$\begin{eqnarray}&&\fl -\left[\Delta -\lambda + \int K(|\mathbf{r}-\mathbf{r}^{\prime}|) \psi_{\mathrm{S}}^2(\mathbf{r}^{\prime})\,\mathrm{d}^2\mathbf{r}^{\prime} \right] \delta v (\mathbf{r}) \nonumber\\ -\psi_{\mathrm{S}}(\mathbf{r})\int K(|\mathbf{r}-\mathbf{r}^{\prime}|)\psi_{\mathrm{S}}(\mathbf{r}^{\prime}) \left[\delta v(\mathbf{r}^{\prime}) +\delta u(\mathbf{r}^{\prime}) \right] \mathrm{d}^2\mathbf{r}^{\prime} & = \kappa \delta v(\mathbf{r}). \end{eqnarray} \tag{ 10b }$$

Real-valued eigenvalues κ of equation (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 equation (10) (which has its origins in the Hamiltonian structure of equation (1)), if κ is an eigenvalue, then −κ as well as ±κ* are also eigenvalues.

Next, we solve equation (10) in order to obtain the internal modes of the second-order radial soliton R₂ 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 equation (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 artifact of the discretization and occurs because the full eigenspace has to be spanned by the eigenvectors. The actual computation of the linear eigenvalue problem equation (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⁸ nonzero elements.

In figure 3, we show the spectrum of the linear stability analysis (BdG equation) for the second-order radial soliton R₂ (a) and the quadrupole soliton Q (b) for mass M = 200. Note that for modes with purely imaginary eigenvalue κ = i Im κ, equation (9) reads $\delta \psi (\mathbf {r},t) = \left [\delta u(\mathbf {r}) + \delta v^*(\mathbf {r})\right ]\mathrm {e}^{-\mathrm {Im}\,(\kappa ) t}$, and it makes sense to write these modes as

$$\begin{equation} \hat{e}(\mathbf{r})=\delta u(\mathbf{r}) + \delta v^*(\mathbf{r}). \end{equation} \tag{ 11 }$$

Only the second-order radial state R₂ is unstable, and we name the two unstable internal modes $\hat {e}_1$, $\hat {e}_2$. The unstable modes $\hat {e}_1$, $\hat {e}_2$ ought to be degenerate for symmetry reasons, the small splitting of the eigenvalues (κ₁ ≈ −2.7i, κ₂ ≈ −2.5i) is again a numerical artifact due to the discretization of the eigenvalue problem equation (10). Interestingly, the shapes of the unstable eigenmodes $\hat {e}_1(\mathbf {r})$, $\hat {e}_2(\mathbf {r})$ resemble quadrupoles. In fact, for practical purposes (see next section) as well as to verify these findings we furthermore solved the eigenvalue problem equation (10) for R₂ on a radial grid [32] with eightfold resolution. Then, instead of stable and unstable quadrupoles, one finds stable and unstable vortices with topological charge m = ± 2 and |κ| ≈ 2.74. The vortices corresponding to m = 2 and −2 can be superposed to again yield the quadrupoles $\hat {e}_1$, $\hat {e}_2$ found already with the full 2D solver, but with much higher precision. Because equation (10) is linear, the amplitudes of the $\hat {e}_j$ are not fixed, and we normalize the latter according to

$$\begin{equation} \int \hat{e}_j^*(\mathbf{r}) \hat{e}_j(\mathbf{r}) \mathrm{d}^2\mathbf{r} = 1,\quad j=1,2. \end{equation} \tag{ 12 }$$

Zoom In Zoom Out Reset image size

The quadrupole soliton Q in figure 3(b) is stable, because all complex eigenvalues correspond to trivial modes and hence the complex form of these eigenvalues is a numerical artifact as discussed above. However, the quadrupole becomes linearly unstable for M ≲ 90, as indicated in figure 2 by dashed lines. In figure 4, we show the results of our numerical stability analysis for the quadrupole soliton Q with mass M = 85. The unstable mode $\hat {e}_1$ with κ₁ ≈ −1.2i resembles the second-order radial soliton R₂, i.e. a hump with a (modulated, i.e. not rotationally symmetric) ring. It is quite interesting to note that both solitons, R₂ and Q, feature internal modes resembling each other. Moreover, in unstable configuration, exactly those modes show the largest growth rate. In other words, R₂ tends to decay in the direction of Q, and Q, if unstable, tends to decay in the direction of R₂. This strong link between R₂ and Q is also reflected by successful variational description in terms of generalized Hermite–Laguerre–Gaussian (HLG) beams [33, 34]. In order to distinguish these variational HLG solutions from our numerically exact solutions, we introduced the notation R₂ and Q.

Zoom In Zoom Out Reset image size

The typical dynamics for R₂ (here for M = 200) as an initial condition is shown in figure 5(a). To determine the shape of R₂, we use the iterative solver mentioned above on a grid containing 400 × 400 points, and we use the same grid for the actual propagation. As we have seen in section 3, the second-order radial soliton R₂ is unstable over the whole range of mass M and therefore any perturbation, that has a nonzero overlap with the unstable internal modes $\hat {e}_1$, $\hat {e}_2$ will lead to an exponential growth of the latter. Practically, the residual in numerical determination of R₂ 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 $\hat {e}_1$ as initial perturbation with tiny amplitude  ∼ 10⁻⁴ to the soliton R₂ to control the breakup in a preferred direction. For small times the dynamics is governed by the exponential growth of the unstable internal mode $\hat {e}_1$, while for later times the evolution becomes highly nonlinear, exhibiting oscillations between R₂ (see inset (α) in figure 5(a)) and a state that resembles the quadrupole soliton Q (see inset (β) in figure 5(a)) [26]. This state (β) we will call the 'turning point'. In the following we will examine in detail the origin and properties of these oscillations.

Zoom In Zoom Out Reset image size

4.1. Projection methods

Let us now introduce the projection method mentioned in the introduction [29, 30] and adopt it to our problem. To this end, we recall the scalar product of two complex functions f and g, defined as

$$\begin{equation} \langle f,g\rangle = \int f^*(\mathbf{r}) g(\mathbf{r})\,\mathrm{d}^2 \mathbf{r}. \end{equation} \tag{ 13 }$$

Obviously, the (unstable) internal modes $\hat {e}_{j}$ of R₂ introduced before (see figure 3) are not orthogonal to their complex conjugate (stable) $\hat {e}_{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 $\{\hat {e}_{1}^*,\,\hat {e}_{2}^*\}$ respectively $\{\hat {e}_{1},\,\hat {e}_{2}\}$ are not mutually orthogonal. Thus, straightforward projections onto $\hat {e}_{j}$ and $\hat {e}_{j}^*$ do not help to elucidate the propagation dynamics. To overcome this difficulty we introduce a set of functions which is biorthogonal to ${\hat {e}_{j},\hat {e}_{j}^*}$ using a Gram–Schmidt-like technique as follows. Firstly, we define

$$\begin{equation} e_{j\perp}=\hat{e}_{j}-\langle \hat{e}_{j}^*,\hat{e}_{j}\rangle\hat{e}_{j}^*, \end{equation} \tag{ 14 }$$

which is simply the projection of the unstable eigenmode $\hat {e}_{j}$ onto the orthogonal complement of the stable eigenmode $\hat {e}_{j}^*$. Secondly, we note that $(e_{j\perp })^*=\hat {e}_{j}^*-\langle \hat {e}_{j},\hat {e}_{j}^*\rangle \hat {e}_{j}$ corresponds to projection of the stable eigenmode $\hat {e}_{j}^*$ onto the orthogonal complement of the unstable eigenmode $\hat {e}_{j}$. Then, it is easy to verify biorthogonality of e_j⊥,(e_j⊥)* with respect to ${\hat {e}_{j},\hat {e}_{j}^*}$:

$$\begin{equation} \langle \hat{e}_{j},(e_{j\perp})^*\rangle = 0, \end{equation} \tag{ 15a }$$

$$\begin{equation} \langle \hat{e}_{j}^*,(e_{j\perp})^*\rangle \neq 0, \end{equation} \tag{ 15b }$$

$$\begin{equation} \langle \hat{e}_{j},e_{j\perp}\rangle \neq 0, \end{equation} \tag{ 15c }$$

$$\begin{equation} \langle \hat{e}_{j}^*,e_{j\perp}\rangle = 0. \end{equation} \tag{ 15d }$$

In figure 6(a) a schematic sketch of the relation between ${\hat {e}_{j},\hat {e}_{j}^*},e_{j\perp }$ and (e_j⊥)* is depicted, and figures 6(b)–(c) show $\hat {e}_1$ and e_1⊥ 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.

Zoom In Zoom Out Reset image size

In order to analyze the propagation dynamics of a solution ψ(x,t) of equation (1), we introduce the quantities

$$\begin{equation} U_j = \langle e_{j\perp},\psi\rangle, \quad S_j = \langle (e_{j\perp})^*,\psi\rangle. \end{equation} \tag{ 16 }$$

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₂, the two unstable eigenvectors $\hat {e}_{1}$, $\hat {e}_{2}$ are degenerate (due to rotational symmetry about the origin), therefore we introduce the rotationally invariant projected variables

$$\begin{equation} U(t) = \sqrt{\sum_{j=1}^2 |U_j|^2}\,, \quad S(t) = \sqrt{\sum_{j=1}^2 |S_j|^2}. \end{equation} \tag{ 17 }$$

Then, any pair of wavefunctions ψ₁(x,t) and ψ₂(x,t) related through a rotation amounts to the same value of U(t) and S(t). For a rigorous proof of this statement, see the appendix.

4.2. Indication of homoclinic connections

Figure 5(b) illustrates the dynamics shown in figure 5(a) in the variables S(t),U(t) introduced in equation (17). We clearly see the second-order radial soliton R₂ (α) decaying into a quadrupole-like state (β), the 'turning point', and then coming back to R₂. In the vicinity of R₂, the decay starts via the local unstable eigenspace E^u (i.e. U(t) > 0, S(t) ≈ 0), and the revival of R₂ 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₂ 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₂ both in the t → ∞ and t → − ∞ limit. The time-span, in which the solution remains close to its initial state R₂, i.e. the homoclinic point (α) in figure 5(b), with practically zero velocity, corresponds to intervals with maximum (nearly) constant peak-intensities in figure 5(a). Because we added a small perturbation in the direction of the eigenmode $\hat {e}_1$ to the initial condition R₂, 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 toward R₂ along E^s, there is always a small perturbation along the unstable eigenspace E^u and the trajectory leaves the neighborhood of R₂ 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.

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₂,

$$\begin{equation} w(t)=\frac{\left|\langle R_2,\psi\rangle\right|}{\langle R_2,R_2\rangle}. \end{equation} \tag{ 18 }$$

Obviously, for ψ = R₂ 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 figure 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 'turning-point' (β), which still maintains a nonzero projection on R₂ and the quadrupole soliton Q(γ) becomes evident. Thus, the homoclinic orbit does not connect the two stationary nonlinear solutions R₂ and Q, as one would expect from previous variational calculations involving approximate HLG beams [26]. Instead, trajectories starting from (α), respectively, (γ) stay well separated (see red, respectively, blue curve in figure 5(a)).

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

5.1. Quasiperiodic motion

The homoclinic orbit is obtained by (slightly) perturbing the initial wavefunction of R₂ in the direction of one of the unstable modes (e.g. of $\hat {e}_{1}$) and integrating equation (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 figure 5) at time t_t into a part parallel to the quadrupole soliton Q and a remainder L

$$\begin{equation} \psi(\mathbf{r},t_{\mathrm{t}}) = c_Q Q(\mathbf{r}) + L(\mathbf{r}), \end{equation} \tag{ 19 }$$

where c_Q = 〈Q,ψ(r,t_t)〉/〈Q,Q〉 was introduced⁶. Perturbed wavefunctions ψ_Γ(r) are then constructed through

$$\begin{equation} \psi_\Gamma^\prime(\mathbf{r}) = c_Q Q(\mathbf{r}) + \Gamma L(\mathbf{r}), \end{equation} \tag{ 20a }$$

$$\begin{equation} \psi_\Gamma = \sqrt{\frac{\langle \psi,\psi\rangle}{\langle \psi_\Gamma^\prime,\psi_\Gamma^\prime\rangle}} \psi_\Gamma^\prime, \end{equation} \tag{ 20b }$$

where Γ parameterizes mixed states between R₂ and Q, and, in equation (20b), the wavefunction was normalized. For Γ = 1, we obtain the 'turning point' (β), 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 figure 7, which indicates quasiperiodic behavior for small times (up to t ≃ 25). The time spend by this orbit close to R₂ is much smaller than for the homoclinic connection of the previous section. This becomes apparent when comparing the peak-intensity evolution in figure 7(a) with the one in figure 5(a). 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 when R₂ is approached, while the 'velocity' approaches zero). On the other hand, in the intensity representation, figures 7(c)–(h), the difference between homoclinic and quasiperiodic behavior is much harder to discern. Propagation in figures 7(a) 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.

Zoom In Zoom Out Reset image size

On the other hand, the dynamics for Γ = 0.99, shown in figure 8, appear again quasiperiodic (see also the discussion of the Fourier spectra in section 5.2), but in this case the orbit appears to be stable, as it persists at least up to t = 1500. The qualitatively different behavior for Γ = 1.01 and 0.99 with respect to stability further corroborates the importance of the homoclinic solution Γ = 1.00 (β). 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 behaviors, because it is just a one-dimensional line in the highly dimensional parameter space.

Zoom In Zoom Out Reset image size

Let us finally consider the trajectory in figure 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 figure 3(b), 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 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 figure 9 for Γ = 0.5. Thus, even though for Γ = 0.5 we are no longer in the region where linear perturbation analysis of the quadrupole soliton Q holds, we still find qualitatively similar dynamics. We note that in the same system equation (1), quasiperiodic nonlinear solutions (so-called azimuthons) linked to stable internal modes of solitons were reported earlier [32, 35].

Zoom In Zoom Out Reset image size

To sum up, we have identified a family of stable quasiperiodic solutions to equation (1), starting from ψ_Γ given in equation (20) 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₂ (Γ = 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.

5.2. 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 $\mathcal {F}$ of the wavefunction at a fixed point in space (in our case the origin r = 0):

$$\begin{equation} f(\omega) = \left|\mathcal{F}\left[\psi(\mathbf{r}=0,t)\right]\right|^2. \end{equation} \tag{ 21 }$$

For a bright soliton solution of the form equation (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 0.5, as can be seen in figure 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₂ is in principle infinite. In practice, this time is greatly affected by numerical noise and the spectrum appears continuous (see figure 10(a)). 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 figures 10(b) and (c))⁷.

Zoom In Zoom Out Reset image size

Thus, the Fourier spectra yield an additional indication of the qualitatively different nature of the dynamics of section 4.2 from the quasiperiodic motion of section 5.1, providing further support for the conjectured existence of an underlying homoclinic connection in the former case.

In previous works, an oscillatory shape-transformation of modes in nonlocal media has been observed numerically [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₂, 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₂, 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₂ itself and its unstable/stable eigenmodes. Projecting the time evolution of the wavefunction ψ(r,t) (obtained by integrating the nonlocal nonlinear Schrödinger equation) 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₂. 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 [36]. 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 [29, 30]) could be applied to a variety of high- and infinite-dimensional conservative systems.

We would like to mention that during the review process, a very recent experimental observation of self-induced nonlinear optical mode-conversion in a nematic liquid crystal cell was brought to our attention [37].

Here we prove that the quantities U(t), S(t) are rotationally invariant, i.e. they have the same value if we substitute ψ(x,y,t) with $\mathcal {R}(\theta )\psi (x,y,t)=\psi (x\,\cos \theta -y\sin \theta ,x\,\sin \theta +y\cos \theta ,t)$, where $\mathcal {R}(\theta )$ is an SO(2) rotation.

The eigenproblem equation (10) for the ring soliton R₂ is rotationally symmetric and, as a result, its internal modes $\hat {e}_{1},\,\hat {e}_{2}$ transform according to

$$\begin{equation} \mathcal{R}(\theta) \hat{e}_{i} = \sum\limits_{j=1}^2D_{ji}(\theta)\,\hat{e}_{j}\,, \end{equation} \tag{ A.1 }$$

where D(θ) is a 2D matrix-representation of SO(2). The explicit representation D(θ) depends on the basis $\hat {e}_{j}$, but for our purposes it is sufficient to show that we have a real representation. We begin by noting that the constraints of orthogonality, D^TD = 1, and unit determinant, $\det (D)=1$, lead to the following general form

$$\begin{equation} D(\theta)= \left( \begin{array}{@{}ll@{}} \alpha(\theta) & \beta(\theta) \\ -\beta^*(\theta) & \alpha^*(\theta) \end{array} \right), \end{equation} \tag{ A.2 }$$

where the functions a(θ), β(θ) are related through

$$\begin{equation} \det\left(D(\theta)\right)= |\alpha(\theta)|^2+|\beta(\theta)|^2=1. \end{equation} \tag{ A.3 }$$

On the other hand, using $\hat {e}_{2}=\mathcal {R}(\theta _0)\hat {e}_{1}$, where θ₀ is the angle that rotates $\hat {e}_{1}$ onto $\hat {e}_{2}$, we can express all matrix elements $D_{ji}=\langle \hat {e}_{j},\mathcal {R}(\theta )\hat {e}_{i}\rangle$ in terms of D₁₁,

$$\begin{equation} D(\theta)= \left( \begin{array}{@{}cc@{}} \alpha(\theta) & \alpha(\theta+\theta_0) \\ \alpha(\theta-\theta_0) & \alpha(\theta) \end{array} \right). \end{equation} \tag{ A.4 }$$

Comparing with equation (A.2) we conclude that α(θ) = α*(θ) and thus our representation is real, and that α(θ − θ₀) = −α(θ + θ₀).⁸

Using equations (A.1)–(A.3) in definition equation (14), along with the relation $\langle \hat {e}_{1}^*,\hat {e}_{1}\rangle =\langle \hat {e}_{2}^*,\hat {e}_{2}\rangle$, one can show that

$$\begin{equation} \mathcal{R}(\theta) e_{i\perp} = \sum\limits_{j=1}^2D_{ji}(\theta)\,e_{j\perp}. \end{equation} \tag{ A.5 }$$

Then, using equations (A.2)–(A.5), it is easy to show that

$$\begin{eqnarray*}&&\fl \bar{U}^2(t) \equiv \left|\langle e_{1\perp},\mathcal{R}(\theta)\psi\rangle\right|^2+\left|\langle e_{2\perp},\mathcal{R}(\theta)\psi\rangle\right|^2 = \left|\langle \mathcal{R}(-\theta)e_{1\perp},\psi\rangle\right|^2+\left|\langle \mathcal{R}(-\theta)e_{2\perp},\psi\rangle\right|^2\\ &&= \left|\langle e_{1\perp},\psi\rangle\right|^2+\left|\langle e_{2\perp},\psi\rangle\right|^2 = U^2(t). \end{eqnarray*}$$

A similar proof holds for S(t).