Probing quasi-integrability of the Gross-Pitaevskii equation in a harmonic-oscillator potential

Previous simulations of the one-dimensional Gross-Pitaevskii equation (GPE) with repulsive nonlinearity and a harmonic-oscillator trapping potential hint towards the emergence of quasi-integrable dynamics -- in the sense of quasi-periodic evolution of a moving dark soliton without any signs of ergodicity -- although this model does not belong to the list of integrable equations. To investigate this problem, we replace the full GPE by a suitably truncated expansion over harmonic-oscillator eigenmodes (the Galerkin approximation), which accurately reproduces the full dynamics, and then analyze the system's dynamical spectrum. The analysis enables us to interpret the observed quasi-integrability as the fact that the finite-mode dynamics always produces a quasi-discrete power spectrum, with no visible continuous component, the presence of the latter being a necessary manifestation of ergodicity. This conclusion remains true when a strong random-field component is added to the initial conditions. On the other hand, the same analysis for the GPE in an infinitely deep potential box leads to a clearly continuous power spectrum, typical for ergodic dynamics.


Introduction
Integrability, relaxation, and thermalization of many-body systems are intricately-linked key topics of the modern theory of non-equilibrium dynamical systems. Although, strictly speaking, a closed quantum system should exhibit no thermalization in the usual sense, non-integrable closed systems can nonetheless mimic relaxation to thermal equilibrium through dephasing occurring within the eigenstate thermalization hypothesis [1,2].
The investigation of these issues has recently become a core activity in studies of dynamics of ultracold gases [3], due to the uniquely precise experimental control achieved in this field. Such settings can be engineered in both weak-and strong-interaction regimes, the effectively one-dimensional (1D) realizations being of particular relevance, as the respective model equations may be able to support integrable dynamics. In this context, pioneering experiments with ultracold atoms in the effectively 1D regime have revealed evidence for a long-term absence of thermalization [4], attributed to the expected integrability of the underlying (Lieb-Liniger) model of strong interactions. Subsequent works, however, have predicted time scales for the breakdown of integrability in experimentally relevant geometries, with thermalization possible through virtual excitation of higher radial modes [5,6]. These findings were reported to be consistent with both the previous experiments [4] and other relevant observations [7], subsequent work also addressing the emergence of pre-thermalization [8], in which a closed system loses part of its initial information.
In the idealized setting described by an integrable equation, which possesses an infinite number of conserved quantities, the trajectories are weakly sensitive to initial conditions, lying on invariant tori in the phase space, realistic systems often exhibit 'weak integrability breaking', in the sense that one can construct and probe 'quasi-conserved' quantities. One should here distinguish between two different issues: the perceived presence of (quasi-)integrability of a given physical system as probed in experiments, and the emergence of integrability in the equations believed to accurately describe the physical system, which is usually probed through numerical simulations. The fundamental equation describing ultracold atoms in the weakly-interacting regime is the nonlinear Schrödinger equation (NLSE), with cubic nonlinearity arising from inter-atomic collisions, alias the Gross-Pitaevskii equation (GPE). This equation is the workhorse of the theoretical studies of ultracold atoms, with an impressive portfolio of successes in predicting experimental phenomena to high accuracy, including the static characteristics of the ultracold gases, their modes, nonlinear waves, dynamical instabilities, etc [9][10][11]. The most common case to which we limit our study here is when the effective interactions are repulsive (i.e., the respective nonlinearity is defocusing). Such an equation is known to be integrable in the 1D free space (including the case of periodic boundary conditions) [12][13][14][15], but not in the presence of the harmonicoscillator confining potential, which is relevant for modeling actual experiments. Even in this case, however, long time simulations of the 1D GPE have revealed no conclusive evidence of chaotization [16,17], which is believed to originate in the experiment from the coupling to transverse degrees of freedom, beyond the limits of the 1D approximation [6]. On the other hand, a single particle in the harmonic-oscillator trap is commonly known to be integrable. The question then arises under what conditions, and to what extent, features of the integrability may be approximately preserved in many-body systems trapped by this potential.
The closest many-body state which exhibits some particle-like properties is a solitonic excitation-specifically, a dark soliton in the case of repulsive interactions, which is thus a natural candidate to use as a probe of the integrability. Importantly, previous studies of the motion of dark solitons in the harmonic-oscillator-trapped 1D GPE lead to a quasi-periodic evolution, revealing no evidence of chaotization (ergodicity) in the evolution of the mean-field wave function, unlike certainly non-integrable settings, corresponding to other (anharmonic) probed trapping potentials [18,19]. This observation suggests an 'apparent quasi-integrability' of the 1D GPE in the harmonic-oscillator trap, with regard to the motion of a dark soliton, which was predicted to perform shuttle motion, as a classical particle, with a well-defined oscillation amplitude and frequency [20][21][22][23][24][25][26] . This behavior is consistent with experiments which have generated dark solitons and demonstrated their motion in elongated quasi-1D BECs [27,28]. However, the presence of any potential, including the harmonic-oscillator trap, is known to break the integrability of the underlying GPE, and, in particular, to trigger the emission of small-amplitude excitations ('sound waves') from dark solitons moving with acceleration [18,22,23,26,[29][30][31]. This mechanism of the decay of dark solitons into radiation is similar to that known for optical dark solitons governed by the NLSE [32][33][34].
To reconcile these apparently contrasting predictions, one implying the presence of the effective quasi-integrability, and the other referring to the non-integrability of the GPE with the harmonic-oscillator potential, it was proposed that the emission of sound waves might be reversible, i.e., that the dark soliton may reabsorb the emitted waves, thus stabilizing itself against the systematic decay [19,26,30] . This effect may even be employed to preferentially stabilize dark solitons in states with selected energies [35]. The reversibility effect has been shown to be crucial over time scales shorter than those imposed by other non-integrability factors (for example, those related to thermal dissipation and coupling to the transverse dimensions) in harmonic-oscillator-trapped BECs [20,21,35,36]. In turn, the sound emission reversibility suggests that the harmonicoscillator potential may maintain quasi-integrability of the system. Further evidence to support this conjecture comes from simulations which reveal a systematic decay when the harmonic-oscillator potential is altered, and the quasi-integrability is clearly broken, e.g., by the addition of dimple traps [22,26], an optical lattice [18,37], or a localized obstacle [38,39]. Another sign of the quasi-integrability in the presence of the harmonicoscillator confinement is an essentially elastic character of collisions between two trapped dark solitons, observed in direct simulations and verified experimentally [28,40].
To gain insight into the presumably quasi-integrable dynamics, we here develop a finite-mode approximation for the 1D GPE with the harmonic-oscillator potential, known as the Galerkin approximation [41]: the wave field is expanded over the full set of eigenmodes of the linear Schrödinger equation with the harmonic-oscillator potential, thus replacing the underlying cubic GPE by a chain of nonlinearly coupled ordinary differential equations for the evolution of amplitudes of the eigenmode expansion. The chain is truncated for a finite set of M modes, sufficient to provide an accurate approximation for the global evolution of the mean-field wave function governed by the GPE, including relevant features such as the above-mentioned sound emission and absorption by the dark soliton. A similar expansion approach was developed for various nonlinear models [42], including multicomponent and multi-dimensional GPE systems [43,44].
The finite-mode Galerkin expansion is also at the heart of the projected Gross-Pitaevskii equation (PGPE) [45], which has been extensively applied to model Bose gases at finite temperatures. However, there are several contextual differences between our study and those using the PGPE.
Specifically, we seek to approximate the zero-temperature GPE wave field, not a thermal field, with the key point being that we can very accurately capture the soliton dynamics and aspects of quasi-integrability by employing M=16 modes.
The aim of our analysis, performed in the framework of the suitably truncated finite-mode dynamical system, is to highlight the degree of the quasi-integrability of the underlying cubic GPE including the harmonic-oscillator potential. Specifically, we find, with high numerical accuracy, that the power spectrum of all dynamical trajectories remains quasidiscrete in the course of the indefinitely long evolution, corresponding to a quasi-periodic motion, rather than to chaotic dynamics. This observation strongly suggests that the Galerkin approximation system with a finite number of the degrees of freedom has almost all its trajectories spanning invariant tori, in accordance with the Kolmogorov-Arnold-Moser theorem [46,47]. Such a finding provides an adequate explanation of the effective quasi-integrability featured by the underlying harmonically-trapped GPE in the previously reported direct simulations [26].
Although the analysis reveals strong evidence of the repeated reversible cycles of the emission/absorption of radiation from/by the dark soliton, there is no straightforward way to isolate the soliton and sound modes through the Galerkin approximation in the condensate trapped in the harmonic-oscillator potential. To demonstrate the role of this process in a more explicit form, we also develop a similar analysis for the GPE in a potential box with zero boundary conditions (i.e., an infinitely deep rectangular potential, which can be experimentally realized using electromagnetic fields [48], although the box walls in the experiment are softer than the ideal impenetrable ones). The Galerkin approximation for the potential box can be naturally built on the basis of the underlying sine and cosine eigenfunctions [49]. Systematic simulations of the GPE in the infinitely deep box show that a moving dark soliton shuttles back and forth in a stable manner (although the soft walls may cause an instability and sound emission [50]). Interestingly, and perhaps somewhat unexpectedly, we find in this case that the power spectrum of generic trajectories is continuous, in direct contrast to the quasi-discrete spectrum found in the harmonic-oscillator potential, which clearly suggests chaotization (ergodicity) of the dynamics in the box, rather than evolution guided by invariant tori. Such behavior is also wholly captured by our finite-mode expansion (without the need for using the known exact box eigenstates of the nonlinear equation [51]).
It is relevant to mention that, strictly speaking, the NLSE in a finite interval with zero boundary conditions belongs to the class of integrable equations [52]. This fact seems to be in contradiction with the above-mentioned ergodicity revealed by the Galerkin approximation for the potential box. However, it is known that there are two types of integrability, strong and weak [53]. The contradistinction between them is based on the relation between the numbers of degrees of freedom and dynamical invariants. Indeed, a Liouville-integrable dynamical system with a finite number of degrees of freedom must have it equal to the number of dynamical invariants [54]. In the limit of the infinite number of degrees of freedom (integrable PDEs), the set of dynamical invariants is also infinite, but in the case of weak integrability the set is incomplete (not 'sufficiently infinite'), which allows the system to feature nonintegrable dynamics, such as fission and merger of solitons in the weakly integrable three-wave system [55], another known example of weak integrability being provided by the Kadomtsev-Petviashvili-I equation. Models of this type, in spite of their formal integrability, readily admit chaotic dynamics-in particular, in the form of wave turbulence in the free space [53]. Of course, the finite-mode truncation most plausibly breaks the strong and weak integrability of the underlying partial differential equation alike, but the concept of the weak integrability suggests a possible explanation to the fact that the truncation, derived for a weakly integrable model, may feature ergodicity: if the underlying model admits chaotic dynamics, the truncated version may feature it too. Concerning the GPE in the finite-size box, the issue of its strong/weak integrability is not explored yet, to the best of our knowledge. This issue may be a subject for a separate study, which is definitely beyond the scope of the present work.
The rest of the paper is organized as follows. In section 2, we summarize the Galerkin approximation for both the harmonic-oscillator and box traps and demonstrate its success in capturing both the ground state solutions and dark soliton motion, with only a small number of modes in the truncation (technical details are given in appendices A and B). This finding enables us to use the motion of the dark soliton as a probe for the quasi-integrability of the 1D harmonicallyconfined GPE, focusing in section 3 on the distinction between the quasi-discrete and continuous spectra of the evolution of complex amplitudes of the Galerkin truncation, which are found, respectively, in the harmonic-oscillator and box traps. Our findings are summarized in section 4.

The Galerkin approximation and its validity
Our analysis starts from the well-known 1D GPE, written in the presence of an arbitrary time-independent potential, V(x) Here Ψ(x, t) is the mean-field wave function of the BEC, normalized to the number of particles dx , and g is the coefficient of the cubic nonlinearity, induced by the van der Waals interactions between atoms which make up the BEC. The characteristic energy, length and time scales are the chemical potential μ of the BEC, the healing length is the speed of sound, and m is the atomic mass. Using these scales to define dimensionless energy, position and time variables, equation (1) can be rewritten as and σ=+1 and −1 corresponds to the repulsive and attractive nonlinearities, respectively. The dimensionless wave function Thus, for typical experimental parameters the atom number will correspond to N 10 4 ~. From this point on, we drop the tilde notation for dimensionless variables; the exception is in figures, where x and t are presented in the dimensional form.
In this work, we are interested in repulsive interactions which admit dark solitons trapped in the external potential [22], therefore we fix σ=1. We consider harmonic-oscillator and box potentials, which are defined, respectively, as In the former case, ω x ≡1 is fixed by rescaling. It is relevant to note that the infinite-box potential, which gives rise to zero boundary conditions, ψ (x=0)=ψ (x=L)=0, directly applies, in addition to BEC, as the model of a metallic conduit for microwaves [56]. The Galerkin approximation takes two different forms, depending on the potential considered. In each case, the wave field is approximated by an M-mode linear combination of time-dependent eigenmodes of the corresponding linear Schrödinger equation, with each eigenmode subject to the unitary normalization. In the harmonic-oscillator case, the corresponding ansatz is where a n (t) are complex amplitudes, which are slowly varying functions of time, in comparison with i n t exp , and amplitudes a n (t) being slowly varying functions in comparison with iE t exp n -( ). Note that it would also be possible to conduct the expansion using the nonlinear modes of the system. However, while the nonlinear solutions in the 1D box potential are known [51], they are not known in analytical form for the harmonic potential. Moreover, the nonlinear box solutions would give rise to extremely cumbersome differential equations. Hence for the purpose of this work we focus on expansion into linear modes, which are clearly demonstrated to capture the key features of the underlying GPE they are approximating.
Evolution equations for amplitudes a n can be readily derived by means of the variational principle [57,58]. To this end, we use the Lagrangian of equation (2) The substitution of ansätze defined by equations (4) and (5) into the Lagrangian leads to the following result and function f is a combination of quartic terms, depending on the number of modes kept in the Galerkin approximation. Accordingly, the dynamics are governed by the Euler-Lagrange equations derived from the Lagrangian This is a mechanical system with M degrees of freedom and two dynamical invariants, H and the total norm We employed a Crank-Nicolson method to find stationary and dynamic solutions to the GPE in MATLAB. Typical simulation parameters are (in the scaled units): spatial discretization Δ x=0.05, simulation-box length L=20, and time step Δt =0.001. The finite-mode dynamical system based on equation (9) was also solved with the help of MATLAB. The largest mode number considered in this work is M=16, the consideration of still larger M being technically possible, but not really necessary, as shown by the results presented below. The MATLAB codes used to solve equation (9) for M=4, 8 and 16, in the presence of the harmonic-oscillator and box potentials, are available online [59].
In appendix B we evaluate accuracy of the Galerkin approximation in capturing the ground state of the system. In the case of the harmonic-oscillator trap and M=16, the agreement is almost perfect for all norms considered. The agreement is almost as good for the box trap: for M=16, nearly perfect agreement is obtained for the case of low norm N10. Similarly, with M=16 modes, the truncated Galerkin approximation is sufficient to reproduce in a virtually exact form the evolution of the mean-field wave function in the presence of harmonic confinement over indefinitely long times, thus accurately reproducing the oscillations of a dark soliton in the trap which fully incorporates the previously-characterized emission and re-absorption of sound waves by the dark soliton. Note that, contrary to the harmonic-oscillator model, the agreement for the infinitely deep box does eventually deteriorate, but only in the course of very long time evolution.

Dark soliton dynamics as the testbed for the validity of the Galerkin approximation
Next we test the validity of the Galerkin approximation for the dark soliton through its comparison to the numerically exact GPE solution. In the case of the repulsive nonlinearity (σ>0), which we deal with in this work, the free space GPE (no trap) has a commonly known family of dark soliton solutions, written here in the unscaled units [12] x t n x x vt i v c e , t a n h .
x 0 is the initial position and v the soliton's velocity. Stationary (alias black) solitons, with v=0, have a zero-density notch with a phase slip of π across it. The soliton's energy decreases with increasing speed [33], emulating a particle with a negative effective mass [22]. Figure 1 shows oscillations of a dark soliton in the harmonic-oscillator and box trapping potentials (left and right columns, respectively). The first row displays the evolution predicted by the Galerkin approximation, as produced by the solution of equation (9), with M=16 modes. In the box trapping potential, the soliton trajectory closely approximates a triangular wave, as would be expected. In the harmonicoscillator potential, one might expect that the soliton would closely follows a sinusoidal trajectory; however, the interaction of the soliton with the dynamical background condensate significantly distorts the trajectory. The deformation of the background field reveals evidence for the interaction of the soliton with the sound (propagating excitations), whereas, in the potential box, we observe chaotization of the soliton dynamics at very long evolution times.
It is relevant to mention that the above-mentioned oscillations of the background condensate trapped in the harmonic-oscillator potential (panel (a)) are excited by the sound emission from the accelerating soliton. The wavelength of the sound is comparable to the system size, hence it becomes visible as a dipole oscillation of the cloud (note that the direct visualization of the sound pulse is challenging due to the immediate reflection of the sound from boundaries and re-interaction with the soliton, motivating the use of dimple traps elsewhere to overcome this issue [19,26]). A quasisteady state is thus established, wherein the background modes and soliton maintain constant average amplitude. This equilibrium is attributed to the balance between emission and re-absorption of sound by the soliton [19,26]. Fluctuations in this energy balance are evident in the quasi-periodic acceleration and deceleration of the soliton, visible in the panel (a). Figures 1(b) and (e) display the soliton's center-of-mass motion, with overlaid results produced by the predictions of the Galerkin approximation with M=16 and GPEsimulations. The Fourier transform of these center-of-mass oscillations is displayed in figures 1(c) and (f), where the oscillation frequency of the soliton, ω s , is highlighted, along with frequency ω d (in panel (b)) corresponding to the dipole mode of  (11)) for the harmonic-oscillator trap; N=5, L=20, v=0.6 and x 0 =10 for the box.
small excitations of the condensate as a whole. The latter mode, with (3) for the definition of ω x ), is excited by the motion and sound emission of the dark soliton traversing the condensate [26]. The dark soliton in the harmonic-oscillator potential is known to oscillate at frequency 2 , 13 as shown theoretically [22] and experimentally [28], deep in the Thomas-Fermi limit, corresponding to large N in our notation. The role of the total number of modes, M, of the Galerkin approximation is addressed in appendix B, demonstrating perfect dynamical accuracy of the approximation for M=16 modes. Moreover, the dark soliton's oscillation frequency in the harmonic-oscillator potential indeed approaches, as expected, the value 2 x w as N increases.

Probing quasi-integrability in the 1D harmonicoscillator potential
In this section, we address the challenging issue of detecting quasi-integrability of the GPE with the harmonic-oscillator potential, which is the main reason why the above analysis was undertaken. As is well known, in strictly integrable dynamical systems the power spectrum of the time dependence of dynamical variables (the complex amplitudes, in the present case), a a t where F stands for the Fourier transform, is truly discrete, corresponding to the generic quasi-periodic motion on a surface of an invariant torus, while non-integrable systems feature a conspicuous continuous component in the spectrum, as a result of destruction of the tori [46,60]. We have applied this criterion of the integrability, by analyzing the motion displayed in figure 1, extending the computation of the spectra to a hundred oscillations of the dark soliton. Figure 2 depicts our main findings, both for the evidently discrete spectrum in the harmonic-oscillator trap, and the case of the box potential (left and right images, respectively). The results are represented by power spectra a 0 w ( ) (top panels) and a 1 w ( ) (bottom panels) of the first and second amplitudes of the Galerkin expansion, which are overlaid on the corresponding results of the GPEsimulations. The GPE spectra were produced by computing the corresponding amplitudes as per equation (B.1), and then calculating their power spectra. In both cases, the agreement between the Galerkin approximation and full GPE simulations is impressive, a feature which is also true for higher-order a j w ( ) coefficients (not shown here). A crucial finding is the stark difference in the results for the dark soliton's motion in the two potentials: in the harmonic-oscillator case (left) we obtain a spectrum consisting of extremely sharp peaks, which may be definitely categorized as a practically discrete spectrum, thus representing quasiintegrable dynamics. The tallest peaks in the spectrum can be immediately identified as located at the above-mentioned frequencies ω d and ω s of the dipole mode of the excitations of the condensate as a whole (see equation (13)), and shuttle oscillations of the dark soliton (see equation (12)). The surrounding peaks can be readily identified as combinational frequencies produced by mixing of these two modes. In stark contrast to this, the spectrum in the box trap exhibits a broad peak, which clearly represents a continuous spectrum, typical to non-integrable systems, that give rise to dynamical chaos and ergodicity [60]. Simulations of the dynamical Galerkin system for the potential box model with M<16 produce a similar behavior, but with a growing noise component.
The finite width of the peaks representing the harmonicoscillator potential is attributed to numerical accuracy, the inherent frequency resolution of the discrete Fourier transform for total simulation time T being Δω=2π/T. In the present case, Δω≈0.0005, and, indeed, the width of the peaks is equal to 2Δω .
A physically relevant situation, when a random initial condition is generated for each coefficient j=1, K, M such that a x iy 0 exp where x j is an observation from the random variable X∼U(−1, 1) and y j from Y∼U(0, 2π) (taking care to renormalise according to equation (10)), has been tested too. Typical examples of power spectra found in this case are displayed, for both the harmonic-oscillator and box traps, in figure 3. It is seen that both models keep the character of their dynamics, corresponding to the quasi-discrete and continuous spectra, respectively, in the presence of the strong random component in the input. Thus, the quasiintegrability of the GPE with the harmonic-oscillator potential is a robust property. Figure 3 also shows the density and phase profiles of the initial condition for both wave-fields. This initial condition is akin to a soliton-gas configuration, generated by summation of several dark soliton solutions with random position, phase and velocity. Using this comparison we can describe the nature of the arising peaks. Considering the dynamics in the harmonic-oscillator case, shallow (fast) solitons have an oscillation frequency close to ω s ≈ω x , whereas deep (slow) solitons near the condensate center have an oscillation frequency near to 2 s x w w » , as shown previously. In the dynamics ensuing from this highly nonequilibrium initial state the power spectra displays a complicated mixing of these modes, with an envelope of spectra centered around their average, ω≈0.85 ω x . Similar to figure 2 the envelopes at larger ω are due to mixing of these frequencies.

Discussion
Previous numerical simulations based on the 1D GPE have revealed shuttle oscillations of dark solitons in the harmonicoscillator potential. In the course of the periodic motion, the dark soliton reversibly emits small-amplitude waves ('sound'), being able to fully reabsorb them. No chaotization was observed in the course of indefinitely long simulations of this model. On the contrary to that, GPE simulations with other types of trapping potentials exhibit irreversible evolution and an eventual trend to the onset of dynamical chaos (wave-function 'turbulence') [19,50]. To explain this phenomenology, we have first derived a finite-mode dynamical system, in the form of the Galerkin approximation, based on the truncated expansion of the wave function, governed by the GPE with the repulsive cubic nonlinearity, over the set of eigenmodes of the corresponding linear Schrödinger equation. The comparison of results produced by the Galerkin approximation to those of full GPE simulations shows that the Galerkin approximation for the model with the harmonic-oscillator potential, with M=16 modes, reproduces the full solutions virtually exactly (with fidelity indistinguishable from 1, see appendix B) for indefinitely long evolution times. In the case of the box potential, the Galerkin approximation with M=16 also provides a high accuracy, although, eventually, there emerges a deviation from the GPE solutions at large evolution times. The main finding is that generic trajectories of the Galerkin approximation derived for the model with the harmonic-oscillator potential produce a discrete power spectrum (up to the accuracy of the numerically implemented Fourier transform), which is a remarkable manifestation of the conjectured quasi-integrability. This finding (which remains true in the presence of a strong random-noise component in the input) strongly suggests that, in the underlying dynamical system, virtually all trajectories wind upon invariant tori, only an extremely small share of the tori (if any) being destroyed. It remains a challenge to understand the quasi-integrability of the GPE with the harmonic-oscillator potential at a deeper mathematical level than the explanation offered by the present analysis.
On the other hand, both direct simulations of the GPE and the Galerkin approximation for the model with the box potential produce a continuous power spectrum, in the form of a very broad peak, which clearly implies that the latter system is subject to the (rather slow) onset of chaotization. Thus, this work puts forward an open question concerning the nature of the non-integrable dynamics in the truncated version of the formally integrable system [52].
It may also be interesting to perform a similar analysis to the one performed here for the case of attractive nonlinearity (σ=−1 in equation (6)) focussing on shuttle oscillations of a trapped bright soliton (see also the related work in [61]), as well as for recurrent collisions between two (or several) solitons (the latter setting was experimentally realized in the self-attracting BEC [62] ). Furthermore, for both cases of the self-repulsion and attraction, the analysis may be extended to a two-component GPE with equal strengths of the self-and cross-interactions, which, in the free space, corresponds to the integrable Manakov's system [63]. This system remains integrable too if it includes the Rabi coupling, i.e., linear interconversion between the components [64], which is thus also an appropriate subject for the consideration. The Manakov's system finds the well-known realization in terms of the two-component BEC mixtures [10]. Moreover, this methodology may enable insight into the stability and dynamics of dark solitons within the nonlocal dipolar GPE, an equation which can be realized experimentally through BECs of atoms which possess strong magnetic dipoles; while the nonlocality breaks the integrability of this governing mean-field equation, dark solitons were also found to show quasi-integrable dynamics, both in homogeneous [65][66][67] and trapped systems [68].
Interesting questions are also expected to arise in the development of the Galerkin approximation for the twodimensional (2D) GPE with an isotropic harmonic-oscillator potential (see also [44] for multi-dimensional Schrödinger equations with generalized nonlinearities and damping). In particular, it is known that the 2D model with the attractive nonlinearity and harmonic-oscillator trapping potential makes the trapped fundamental solitons completely stable (against the critical collapse in the 2D space [69]), and provides for partial stabilization of vortex solitons with topological charge 1 against the collapse and splitting [70]. The investigation of the 2D model may be interesting also for the reason that the 2D GPE in the free space is not integrable, the question being if the harmonic-oscillator confinement may induce a quasiintegrability in this case.
Data supporting this work is openly available under an 'Open Data Commons Open Database License' [59].

Acknowledgments
We acknowledge valuable discussions with Tom Billam, Maxim Olshanii, and Alexander Its. NGP acknowledges funding from the Engineering and Physical Sciences Research Council (Grant No. EP/M005127/1). TB acknowledges support from Engineering and Physical Sciences Research Council. The work of BAM on this project was carried out in the framework of the visiting professorship provided by the Newcastle University. This author also acknowledges support provided by grant No.2015616 from the joint program in physics between the NSF and Binational (US-Israel) Science Foundation, and by grant No.1286/17 from the Israel Science Foundation.

Appendix A. The four-mode truncated Hamiltonian and dynamical equations
In this appendix we provide an explicit example of the dynamical system produced by the Galerkin approximation with M=4 modes in the model with the harmonic-oscillator potential. The Hermite polynomials required to construct the Calculation of the quartic term in the corresponding Lagrangian (6) leads to Hamiltonian (8) in the following form: a  a a   a a  a a  a a  a a   a a  a a a a  a a a a   a a a  a a a   a Finally, substituting this Hamiltonian in Euler-Lagrange equation (9), we arrive at the following dynamical system with four degrees of freedom: a e a  a a e a   a a e a  a a e a  a a a   a a a  a a a   a a e  a a e   a a a e  a a a   where  is the actual range of x for the harmonic-oscillator potential, and interval 0<x<L for the box. Using the amplitudes given by equation (B.1) to construct the Galerkin approximation wave function ψ GA (x, 0), as per equations (4) and (5), we define the fidelity F of the approximation as where F=1 (F=0) corresponds to two identical (mutually orthogonal) wave functions. Figure B1 (top row) shows how the number of modes affects the fidelity of the initial Galerkin approximation wave function, for different norms N. In the case of the harmonicoscillator trap and M=16, the fidelity is virtually exactly F=1 , implying an almost perfect GPE-Galerkin approximation overlap for all norms considered. For the box trap, the fidelity is still good but poorer than for the harmonic-oscillator; even in the case of N=1 (weak nonlinearity), ψ GPE for the box is not perfectly approximated by the truncation with M<10.
Finally we explore the validity of the Galerkin model for recreating dynamical simulations of the GPE. Figure B1 (bottom row) explores the role of the mode truncation in the Galerkin approximation, by direct comparison of the oscillation frequencies, extracted from the GPE simulations, and their counterparts, predicted by the Galerkin approximation, as a function of the total norm. The results clearly show that the increasing number of modes, M, improves the Galerkin approximation accuracy for all norms considered.