Approximate nonlinear wave solutions of the coupled two-component Gross–Pitaevskii equations with spin–orbit interaction

Recent experimental observations of spin–orbit coupling (SOC) in Bose–Einstein condensates (BECs) open the way for investigating novel physics of nonlinear waves with promising applications in atomic physics and condensed matter physics. The interplay between atomic interactions and SOC are crucial for the understanding of the dynamics of nonlinear waves in BECs with SOC. Here, in the small linear coupling regime, an approach is presented which allows us to derive an infinite number of novel approximate solutions of the Gross–Pitaevskii equations (GPEs) in one and two dimensions including SOCs, time-dependent external potentials, and nonlinearities leading to breathers and periodic as well as quasiperiodic nonlinear waves. To verify the theoretical predictions we perform numerical simulations which show for several cases a very good agreement with the analytics. For the case of one spatial dimension, it is shown that functions describing the external potential and nonlinearities cannot be chosen independently. The management of the solutions is clarified along with some important physical properties such as Josephson oscillations and Rosen–Zener oscillations.

In recent years, important advancements in the creation and manipulation of laser fields culminated in the production of synthetic gauge fields allowing the emergence of very rich novel physics as explored, e.g., in ultracold atoms [1]. Important examples are the realization of synthetic electric [2] and magnetic [3] fields, the spin Hall effect [4], SOC [5], etc. The observation of an equal strength of Rashba and Dresselhaus SOC types in one-dimensional (1D) BECs has stimulated intensive studies on new effects involving SOCs. In particular, the interplay between nonlinearities and SOCs on the dynamics of nonlinear matter waves such as solitary waves (stripe phases, half vortex gap solitary waves) [6][7][8][9] and exotic complexes (periodic waves) [10] have been investigated. In addition, several other novel phenomena have been analyzed such as the Stueckelberg interferometry [11], Josephson oscillations [12,13], and momentum-space Josephson effects [14], to name just a few. The recent experimental demonstration of SOCs in two-dimensional (2D) BECs [15] promises to realize and explore numerous novel nonlinear matter waves supporting a huge variety of interesting quantum physics. For this goal, the management of nonlinear waves in 1D and 2D BECs with SOCs will be a crucial task and thus needs to be explored.
In this work, in the linear coupling regime where the Raman coupling is small, we predict several new types of approximate analytical nonlinear wave solutions of the spin-orbit coupled Gross-Pitaevskii equations (GPEs). We present a simple and effective mathematical hierarchy to unveil and study an infinite number of new nonlinear wave solutions for the GPEs with SOCs, time-dependent nonlinearities and external potentials in one and two spatial dimensions, respectively. Even though in our analytical 1. Results

Model
Nonlinear waves with time-dependent external potential, nonlinearity, and Raman coupling are described by the GPEs with SOCs with j = 1, 2. Some aspects of this important problem have been studied recently (see [11][12][13][14] where the external potential and nonlinearities are constant but the Raman coupling varies in time and [16] where the Raman coupling is set to zero but the nonlinearities and the external potential vary in time), however, so far no complete treatment is available. Equation (1) are written in dimensionless form [17] in D-dimensions (D = 1, 2, 3) in which Ψ j is a D-dimensional wavefunction depending on t and X ∈ R D . The coordinate X has the explicit form: X = (x, y, z) T for D = 3, X = (x, y) T for D = 2, and X = x for D = 1. As external potential we use the typical harmonic potential with a time dependence of the form V( being the aspect ratios of the trapping potential, ω i the trap frequencies, and ω ⊥ = max{ω x , ω y , ω z } the transversal trap frequency, k 0 and δ are the SOC strength and the detuning of the Raman transition, while the time-dependent functions g(t) and Ω(t) represent the strengths of the nonlinearities (two-body interactions) and the Rabi frequency of the Raman coupling, respectively. Time, space and energy are measured in units of 1 ω ⊥ , a ⊥ = 2mω ⊥ , and E = ω ⊥ , respectively. In BECs temporal variations of nonlinearities can be achieved by using the Feshbach resonance management technique [18] and time modulations of the Raman coupling are accessible in experiment with a high tunability [19]. Equation (1) with V(x, t) = 0 and constant parameters also appear in several other physical contexts such as nonlinear optics (fiber optics gratings [20,21], birefringent optics fibers [22], coupled optics waveguides [23]), field theory and the massive Thirring model [24].

Analytical results
It is well known that equation (1) are not integrable by the inverse scattering transform method mainly due to the presence of the Raman coupling terms and the external potential. Up to now only perturbative and variational methods have been used to study their solitary wave solutions [6][7][8][9][10]. The main idea proposed in this work is that equation (1) allow to construct 'a sea' of approximate solutions for the 1D case by allowing solutions of the integrable case that correspond to the absence of Raman coupling to rotate in time according to a special transformation. The integrability in the absence of Raman coupling for the 1D case is shown by mapping the system to an integrable Manakov system [25] and for the 2D case to an integrable higher-dimensional time-gate Manakov system [26]. In the small Raman coupling approximation, we therefore use transformations to treat the Raman coupling term and the external potential in equation (1). To achieve the abovementioned mapping, we follow a three-step hierarchy which is described in the following.
(i) Removal of Raman coupling terms: let us write equation (1) in matrix form where I is the identity matrix, σ 1 and σ 3 the Pauli matrices and Ψ = (Ψ 1 , Ψ 2 ) T . The right-hand side of equation (2) represents its inhomogeneous part and the left-hand side its homogeneous part [27,28] which is solved by where Ψ 0 = (Ψ 01 , Ψ 02 ) T is a constant, S(t) = 1 which is equation (1) but now without the Raman coupling terms. The transformation equation (3) leads to equation (4) provided that Ω(t) and s(t) both tend to zero (which also means that both quantities must be small) in addition to Ψ 01 , Ψ 02 being small. This becomes clear after plugging equations (3) into (2) and one realizes that cos[s(t)] → 1 and sin[s(t)] → 0 (a full derivation of equation (4) is given in the appendix below). The limit of such an assumption is that the smallness of Ω(t) and s(t) remains unclear. Nevertheless, approximate values of Ω(t) may be found during numerical simulations. One notices that the unitary transformation U(t) along with the conditions Ω(t), s(t) → 0 and small Ψ 01 , Ψ 02 only 'removes' the Raman coupling terms from equation (1) but does not affect the other terms. It is easy to show that equation (4) is invariant under unitary transformations, U(t) being similar to the unitary transformation of an electron spin. Solutions of equation (2) could be seen as rotating in 'spin space', the rotation being here U(t). Hence, one should expect time periodic variations of solutions of equation (1) induced by the presence of the Ω, even though the latter parameter remains small as in the current work. The transformation used is valid only for Manakov nonlinearity (any difference in self-and cross-modulation terms would lead to nonlinear x-dependent potentials [27,28] Ψ 0 = Ψ 01 ,Ψ 02 T and inserting equation (5) into equation (4) yields Equation (6) represents the coupled GPE model with neither SOC nor Raman detuning terms.
where ξ = x l(t) , dT dt = 1 l 2 with the constraints which converts equation (6) to the 1D Manakov system in the rescaled variables ξ, T with Φ = (Φ 1 , Φ 2 ) T and = ±1. The transformations (7) and (8) hold provided that the function g(t) is the solution of the ordinary differential equation Equation (10) serves as an integrability condition for the solution of equation (6). Hence, approximate solutions of equation (1) may be written as with Φ being a solution of the 1D Manakov system many of which are available in literature. In equation (11), one has dt and l(t) = 1 |g(t)| derived from equation (8). Furthermore, as the nonlinear differential equation equation (10) has depending on the geometry of the external potential many solutions this extends the number of solutions of equation (11). For example, choosing the expression of the nonlinearity management function g(t) according to experimental situations, one determines that of α x (t) from equation (10) and the explicit form of f(t). Another important aspect is amplitude time-modulations of the densities proportional to |g(t)| in addition to periodic time modulations induced by Raman coupling. It should be mentioned that it is possible to obtain a relationship between g(t) and α x (t) by assuming a particular form of α x (t) [29] or using appropriate similarity transformations [30,31].
(iii-(b)) 2D case: in higher spatial dimensions, D = 2, 3, equations (7) and (8) lead to V(X, t) = 0 and |g| = 1 and equation (6) reduces to a higher-dimensional Manakov system. For D = 2, the 2D Manakov system is completely integrable if |g(t)| = 1 as proven in [26]. Moreover, when V(X, t) = V(X) = 0, an infinite number of approximate solutions, i.e., excitations, of equation (6) were found in [26] including solitary waves and periodic and quasiperiodic waves. For D = 3, the 3D Manakov system is to the best of our knowledge not integrable, however, for this case one may use its approximate solutions that are available. The 2D (3D) solutions (approximate solutions) of equation (1) are written as whereΨ 0 is a 2D (3D) solution (approximate solution) of the Manakov system. The transformation equation (5) only removes the SOC and the Raman detuning terms and combined to the unitary transformation U(t) maps equation (1) to a Manakov system (D = 1, 2) when V(X, t) = 0.
Focusing on solutions, one may consider equation (6) and its large classes of solutions in D = 1, 2, 3 spatial dimensions ranging from solitary waves, vortices, to periodic and quasiperiodic waves in order to obtain solutions of equation (1) in the presence of external potentials with different geometries [32][33][34][35][36][37]. This implies that the hierarchy proposed here is valid for a wide range of physical situations some of which have been investigated in binary BECs [38].
The densities n j = |Ψ j | 2 are given by withΨ 0j = |Ψ 0j | exp(ıθ j ). Equations (13) and (14) show that the SOC contributes to spatial oscillations of the densities with frequencies (−1) j k 0 π (j = 1, 2) if θ 1 = θ 2 , which means that both states have opposite spatial oscillations. Such spatial oscillations are responsible for fringes observed in experiment [5]. On the other hand, solutions without SOC-induced spatial oscillations are obtained for θ 2 − θ 1 = 2k 0 x + mπ, m ∈ Z. This implies that solutions without SOC-induced fringes are possible if and only if the solutions of the related GPEs (equation (6)) are out of phase, θ 2 − θ 1 = 0. Stationary solutions can also be obtained ifΨ 01 =Ψ 02 are stationary solutions of equation (6) or ifΦ 01 =Φ 02 are stationary solutions of the 1D Manakov system. Many solutions of the Manakov system including solitary, periodic, and quasiperiodic waves are given in [39,40] and references therein.
Oscillations: equations (13) and (14) show that the presence of the Raman coupling term induces time-periodic oscillations of the densities of each wavefunction with an exchange of atoms between both pseudospin states as was also previously predicted [27,28]. With N j = n j dX (here, integration with respect to each coordinate is implied: dX = dxdydz for D = 3, dX = dxdy for D = 2, dX = dx for D = 1) the number of atoms in each pseudospin state and N = N 1 + N 2 the total number of atoms, the oscillations or transition dynamics from one state to another is determined by the transition rate function R(t) whereñ 0j = Ψ 0j dX. The novelty here is the presence of the SOC term in the expression of R(t) meaning the SOC term alters the exchange of atoms. Thus SOC may be used to precisely manage the transitions from one state to another. Depending on the form of the Raman coupling term, different types of oscillations can be observed.
The period of Raman-induced oscillations is 2π Ω for both the densities of each state and the transition rate R(t).
Case 2: Rosen-Zener oscillations. This situation corresponds to time modulations of Ω. One example of a time-dependent Raman coupling term has recently been realized in BECs with SOCs [11] where Ω(t) = Ω 0 + Ω M cos(2πf mod t), Ω 0 represents the unmodulated Raman coupling term, Ω M the modulation amplitude, and f mod is the modulation frequency. In this case, the oscillation rate, though rather complex, is also given by equation (16) for 2S(t) = [Ω 0 t + Ω M 2πf mod sin(2πf mod t)]. It is clear that our approach is able to provide and thus clarifies the interferometry pattern and rotations of the spin polarization of BECs during transport through the modulation-dressed band observed in [11].

Numerical results
Examples. We now present the results of numerical simulations to verify and analyze several of the analytical predictions described above. We restrict ourselves in this work to the one dimensional case where exact analytical solutions are available. Extension to higher dimensional solutions (2D and 3D) for which only approximate solution exist will be carried elsewhere. In the numerics we use the split-step Fourier method [41,42] with random amplitude perturbations of 0.01 strength of the inserted analytical solutions at initial time. Some of the solutions found above are taken as initial conditions of equation (1)    using for example phase-and density-engineering methods [44] (a wave function engineering technique [45,46]), two Raman lasers at 804.1 nm intersecting with an angle of 90 • are used to create a weak SOC as in [5].
Panels (a)-(e) of figure 1 (for |Ψ 1 | 2 ) and figure 2 (for |Ψ 2 | 2 ) display snapshots at t = 20, 40, 60, 80, 100 of comparisons between our analytical solutions and the numerical ones for the densities of both components, respectively. Panels (e) in figures 1 and 2 show that the very similar time evolution of the analytical and numerical results of the number of atoms of each component, respectively. One notices a very good agreement between our analytical predictions of densities and number of atoms versus results from direct numerical simulations. From figure 3 one observes that both analytical and numerical results fit well with a very high accuracy for the transition rate R(t) (a), while the agreement for the total number of atoms remains also quite close. Clearly, the numerical results on the spatiotemporal evolution of the densities shown in figures 3(c) and (d) demonstrate time periodic oscillations of period 20π (10 ms) induced by the Raman coupling in accordance with our analytical results. SOC-induced spatial oscillations amount to 10π (26.70 μm) and are too large to be observed as the width of the solutions is 0.125 (0.106 μm). The breathers remain stable up to a simulation time of 30 ms for attractive condensates of 7 Li atoms a value which corresponds to 189 in our dimensionless time units and thus may be observable using current experimental techniques. Additionally, figure 3(e) confirms that the total density is a conserved quantity within the transformation used in the above hierarchy.   (9) mentioned above. This is a solution that, to the best of our knowledge, has not been presented previously. However, it is a known exact solution of the Manakov system derived in [34]. Snapshots of densities at t = 27, 40, 60, 80, 100 shown in panels (a)-(e) of figure 4 (|Ψ 1 | 2 ) and figure 5 (|Ψ 2 | 2 ) show a very good agreement between analytical and numerical results; panels (f) in figures 4 and 5 confirms that our analytical predictions capture the numerical ones with a very high accuracy regarding the temporal variations of the number of atoms N 1 (t) and N 2 (t), respectively. Moreover, the transition rate R(t), the total number of atoms and the spatiotemporal evolution of the total density drawn in figures 6(a), (b) and (e) prove that our analytical predictions corroborate with high accuracy to the numerical results. Once again, the spatiotemporal evolutions of both densities exhibited in figures 6(c) and (b) show time periodic oscillations of period 20π (10 ms) induced by the Raman coupling as predicted by our analytical results. Once again, the numerical solutions remain stable up to 30 ms, a value large enough for our solutions to be observed in current experiments. It is obvious from comparisons of figures 1(a)-(e) and 3(c) and also figures 2(a)-(e) and 3(d), where the densities of the numerical solutions are smaller than those of the analytical ones, that the analytical bright-bright initial conditions have some excess energy which is released during the evolution. When such an excess energy called energy radiation is present in solitons of integrable systems it disperses as time evolves and only the solitons or solitary waves persist at longer times [47,48]. Similar energy radiations are also observable when comparing figures 4(a)-(e) and 6(c) and figures 5(a)-(e) and 6(d) for the double bright-bright initial condition.

Discussion
The reduction algorithm presented above enables, starting from a BEC with SOCs, to obtain a Manakov system. Such an idea was developed earlier in [49] for single nonlinear Schrödinger equations with time-dependent dispersion, nonlinearity and gain, then extended to single nonlinear Schrödinger equations with time-dependent dispersion, nonlinearity and external potentials in [50]. For the case of BEC, the reduction algorithm of [50] works with the same integrability condition as ours given by equation (10). The ideas developed in [49,50] were used in many works in the recent past years [29] (and references therein) and extended to single and binary BECs with time and spatial variations of nonlinearities and external potentials [51][52][53][54][55]. Though elaborating on the idea developed in [49][50][51][52][53][54][55], one notices some important differences to our current work. First of all, we used above a system of extended coupled nonlinear Schrödinger equations with the presence of both linear Rabi and SOC terms. The latter terms are not considered in [49][50][51][52][53][54][55]. Secondly, the solutions of each equations are given as superpositions of two solutions of the Manakov system in this work while other approaches were considered in previous works [49][50][51][52][53][54][55]. Thirdly, one must follow, in the order provided, the three-step algorithm presented above to obtain analytical solutions. Another striking difference is that in this work, the amplitudes of the solutions are modulated by g(t) which is related to the strength of the external potential while in [49,50] the dispersion enters the modulation of the amplitudes. Some examples of time modulations of nonlinearities and strengths of harmonic potentials are given by (i) [56].
In the simulations presented above, we focus on breathers (time-periodic solitons or solitary waves). Solitons of BECs with SOCs were also considered in other works. For example, solitons were investigated in BECs with SOC in the absence of external trapping in [57][58][59], in optical lattices [60][61][62][63] and in random trapping [64,65] just to name a few. The authors in [66] analyzed BECs with SOC in harmonic potential with repulsive interactions and focus on the effect of Zeeman splitting. They consider weakly interacting stationary soliton solutions within a bifurcation analysis starting from solutions of the noninteracting system while the solutions considered here are dynamical ones at least in one dimension and are available both for attractive and repulsive interactions. This stems from the fact that our model takes into account linear Raman couplings between the pseudospins which is absent in the work in [65]. Even in the case where one considers the reduction of our model to the one used in [66] by setting Ω = 0, our solutions remain dynamic ones and time-modulated by |g(t)|, a situation quite different from the spatial modulations obtained in [66].
It is well known that solitons or solitary waves including breathers are long lasting nonlinear waves. The simulations shown above present stable breathers after 30 ms that correspond at least for repulsive 87 Rb  BECs to the experimental time duration of collective oscillations observed in [67]. The unitary transformation allows to get time periodic oscillations as observed in experiments [67] that are in pretty good agreement with the numerics. Such agreements suggest that the approximate analytical solutions are close to the exact solutions, hence meaning that our choice of allowing Raman-induced time periodic oscillations works well. Such inclusion of time oscillations to solutions though not exact captures very well the behavior observed in experiment in addition to allowing to open the way to new behaviors and applications. The model has several parameters that may affect the stability of the solutions such as the nonlinearities, external potential, Raman coupling term, and spin-orbit coupling (SOC) strength. For the simulations presented in figures 1-6, we observed during intensive simulations carried that the agreement between analytical and numerical solutions increases for smaller trapping aspect rations ( ω x ω ⊥ ) in addition to Ω 0.1. In the course of the dynamics the solutions get rid of the radiation energy and remain stable at relatively longer times. Figures 7(a)-(c) display stable breather solutions for attractive 7 Li atoms after 100 ms with same parameters as in figure 1. Figures 7(d)-(f) portray a example of stable repulsive 87 Rb atoms after 214.86 ms using realistic experimental parameters (ω x = 1.5 Hz, ω ⊥ = 140 Hz) [67]. As the interplay among the system parameters may alter the stability of solutions, a thorough perturbation analysis of solutions found with the hierarchy proposed here will be carried elsewhere.
The solutions found and analyzed above are well suited to describe important physical situations. For example, we predict above that our hierarchy is able to describe slow Josephson oscillations between atomic states that are controlled by the strength of the Raman coupling Ω. This is confirmed for both solutions with very high accuracy by the numerical simulations shown in figures 3(a) and 6(a), respectively. Josephson oscillations of BECs with SOC may be used to explore and test relevant physical systems, such as topological insulators, topological superconductors, and superfluids [68][69][70], and may be relevant for applications especially the creation of atomic semiconductors with SOC and the fabrication of spin-based atomtronic devices [71]. An atomic state interchange controlled by the Raman coupling strength Ω takes place. A complete atomic state interchange is realized after a duration that is equal to an oscillation period. A striking example of this feature is observed by comparing for example the time evolution of the densities of the double bright-bright solitary waves plotted in figure 8(a) for |Ψ 1 (x, 0)| 2 ( figure 8(b) for |Ψ 2 (x, 0)| 2 )  and figures 4(a)-(e) (figures 5(a)-(e)). The state interchange predicted in this work may be used to store and transfer spin-orbit information from one state to another one. Let us assume one starts from, e.g., a bright solitary wave in the pseudospin up state while the pseudospin down state is empty (absence of atoms). After one period, the pseudospin down state will be filled with atoms, and, as both states satisfy equation (11), the pseudospin state down contains the spin-orbit information. Such a transfer mechanism of information should be useful for applications in quantum information processing and spin-based quantum computation with ultracold atoms [71]. The density oscillations predicted above and presented in figures 3(c) and (d), 4(c) and (d), and 7(a), (b), (d) and (e) may be used to explain spin current generation and collective oscillations observed in [67]. Moreover, periodic solutions of the Manakov system [39,40] not considered here may also be used to explore new resonances that may occur in different cases, for example when Rosen-Zener oscillations are considered. Extensions to 2D and 3D solutions shall allow to address many aspects like new localized solitary waves including multipoles and half-vortex solitary waves, new time-periodic structures on the hand in addition to analyzing new phase transitions and phase diagrams.
As stated above, approximate analytical solutions found using the hierarchy presented in this work are valid provided that Ω(t) and the exact solutions of the integrable case be small. While the condition on solutions may be easily satisfied by taking for example small amplitudes of solutions of the Manakov system, it is useful to relate the condition on Ω to a physical parameter of the system. As the SOC strength k 0 may be chosen proportional to Ω [2], both Ω and k 0 should be smaller than 1. To test this idea, we ran numerical simulations with same parameters as in figure 4 except Ω = k 0 = 1 and compared the results with analytical solutions. Figure 9 displays snapshots of comparisons between analytical and numerical solutions for densities at t = 27, 40, 80. Panels (a)-(c) show the comparisons for |Ψ 1 | 2 while panels (d)-(f) correspond to |Ψ 2 | 2 . In this case, there indeed is a large discrepancy between numerical and our analytical solutions. Conversely, when Ω, and k 0 are much smaller than 1, our analytical solutions corroborate the numerical ones with a higher accuracy. This feature is exemplified in figure 10 where panels (a)-(c) and panels (d)-(f) show comparisons of analytical and numerical solutions for |Ψ 1 | 2 and |Ψ 2 | 2 , respectively.

Conclusions
In summary, the proposed analytic three-step procedure provides a systematic way to unveil and analyze a large number of families of approximate analytical solutions of the GPEs with SOC valid only at small Ψ 01 , Ψ 02 and Raman coupling regime. The nonlinearities, external potential, and the Raman coupling are allowed to vary in time. One important finding is that the GPEs with SOC under the previous assumptions, admit a large number of approximate analytical solutions that may be constructed from known (exact 1D or approximate 2D, 3D) solutions of the Manakov, by following a simple and straightforward hierarchy. Furthermore, Focusing on the 1D case, the stability of a few analytical solutions has been confirmed by numerical simulations. The used mathematical procedure allows us to find new families of solutions including breathers and periodic and quasi periodic waves. Our findings not only apply to BECs but are also relevant for several other fields including nonlinear optics and condensed matter physics. Especially, the hierarchy proposed in this work may be applied to investigate and extend the stability properties of one and higher dimensional nonlinear modes with general gauge fields as suggested by the works in [72,73], in addition to the fact that it also applies to higher dimensional SOCs [74] for which new modulated nonlinear modes may be unveiled.