The phase diagram and stability of trapped D-dimensional spin-orbit coupled Bose-Einstein condensate

By variational analysis and direct numerical simulation, we study the phase transition and stability of a trapped D-dimensional Bose-Einstein condensate with spin-orbit coupling. The complete phase and stability diagrams of the system are presented in full parameter space, while the collapse dynamics induced by the mean-filed attraction and the mechanism for stabilizing the collapse by spin-orbit coupling are illustrated explicitly. Particularly, a full and deep understanding of the dependence of phase transition and stability mechanism on geometric dimensionality and external trap potential is revealed. It is shown that the spin-orbit coupling can modify the dispersion relations, which can balance the mean-filed attractive interaction and result in a spin polarized or overlapped state to stabilize the collapse, then changes the collapsing threshold dependent on the geometric dimensionality and external trap potential. Moreover, from 2D to 3D system, the mean-field attraction for inducing the collapse is reduced and the collapse speed is enhanced, namely, the collapse can be more easily stabilized in 2D system. That is, the collapse can be manipulated by adjusting the spin-orbit coupling, Raman coupling, geometric dimensionality and the external trap potential, which can provide a possible way for elaborating the collapse dynamics experimentally.

is still not clear. Particularly, a full and deep understanding of the dependence of phase transition and stability mechanism on geometric dimensionality and external trap potential is still missing.
The purpose of this work is to study the distinct phases and their stability, and the stability mechanism of a trapped D-dimensional binary BEC with SOC. The collapsing threshold, the phase and stability diagrams and the collapse dynamics are presented by variational methods and confirmed by the direct numerical simulation of Gross-Pitaevskii equation in full parameter space. The main results are summarized in Fig. 1, where the complete stability and phase diagrams are depicted in intra-and inter-species interaction (g-g 12 ) plane for different SOC, RC and geometric dimensionality. For a conventional BEC (i.e., without SOC and RC), the collapse will occurs for the attractive interactions beyond the threshold. However, the strong RC Ω and weak SOC k 0 (i.e., k / 2 0 2 Ω > ) can stabilize the collapsed conventional BEC induced by strong intra-species attractive interaction g, then result in a stable phase I (i.e., the zero-momentum phase). As the increase of Ω k / 0 2 (see Fig. 1(a1,b1,a3,b3)), the critical intra-species attractive interaction value for collapse increases, thus the BEC with a stronger intra-species attraction is also stable, while the region of phase I becomes larger and the region of phase II (the plane wave phase) becomes smaller. On the other hand, the collapsed BEC with strong inter-species attractive interaction g 12 can be stabilized by weak RC and strong SOC (i.e., Ω < k / 2 0 2 ), which generates a stable phase II. As the decrease of k / 0 2 Ω (see Fig. 1(a2,b2,a3,b3)), the threshold of collapse shifts down to stronger inter-species attraction, thus the BEC with a stronger inter-species attraction is also stable, while the regions of phase I (phase II) decreases (increases). Furthermore, a full and deep understanding of the dependence of phase transition and stability mechanism on geometric dimensionality and external trap potential is presented explicitly. Compared with the 3D case ( Fig. 1(a1,a2,a3)), the regions of stable BEC are expanded, the collapse threshold is enhanced, the collapse speed is reduced, and the collapse can be more easily prevented in 2D external trap potential ( Fig. 1(b1,b2,b3)), which can provide a possible way for elaborating the collapse dynamics.

Results
Model. Based on the recent experiment related to the realization of BEC with SOC 5-8 , we consider a binary SOC-BEC loaded into a D-dimensional harmonic trap potential. According to the mean-field approach, the BEC is governed by the following dimensionless Gross-Pitaevskii (GP) equation [5][6][7][12][13][14][15] : 2 characterizes the interatomic two-body interaction with the dimensionless interaction constants g N a m l represents the atomic mass, the number of atom, and characteristic length along z direction with trapped frequency ω z , respectively. The positive (negative) s-wave scattering lengths a ij refers to the interatomic repulsion (attraction). The physical variables are rescaled as l z The single-particle Hamiltonian Figure 1. The stability and phase diagram in g-g 12 plane for D = 3 (a1-a3) and D = 2 (b1-b3). The solid and short dashed line, respectively, refers to stability and phase boundary by means of the variational methods. The dots represent the results of direct numerical simulation of Eq. (1). The solid line divides the interaction plane into two, one is collapse region the other is stable region, while the stable region is divided into phase I (zero momentum phase) and phase II (plane wave phase) by the short dashed line.
Scientific REPORTs | 7: 15635 | DOI:10.1038/s41598-017-15900-w x ∼  is the dimensionless RC constant accounting for the transition between the two spin states, 1 is the dimensionless SOC strength fixed by the momentum transfer of the two Raman lasers, the operator = − ∇ p i is the canonical momentum, σ i are the usual 2 × 2 Pauli matrices, while ω = V r ext 1 2 2 2 is a harmonic trap with 1 ω ∼ for D = 3 due to the spherical symmetry and ω ω ω = / 1 z for D = 2 due to the pancake-shaped cylindrical symmetry where ω  is the transverse trapped frequency. Note that, there is a peculiar property of violating both parity and time-reversal symmetry in Hamiltonian Eq. (2). Here, we have g 11 = g 22 = g in the absence of Zeeman splitting. Based on experiment 5,39,40 , for 87 Rb atoms, corresponding parameters can be widely adjusted, which has already been estimated in the previous work 37 .
From the scale analysis, the amplitude of normalized spinor wave function can be estimated, i.e., Ψ ∼ − L D/2 , where L is the dimensionless characteristic size of BEC. Therefore the dependence of the system's energy on the characteristic size is obtained as , where coefficient c kin , c trap and c soc are all positive, while positive (negative) c int intra and c int inter refers to intra-and inter-species repulsive (attractive) interaction, respectively. It is well known that the external trap potential changes the dispersion relations, which can balance the mean-field repulsive interaction and result in stable BEC instead of diffusion. Similarly, the SOC can also modify the dispersion relations, which suggests a way to neutralize the mean-field attractive interaction to stabilize the collapsed BEC dependent on the external trap potential and the geometric dimensionality. A full and deep understanding of the dependence of phase transition and stability mechanism on the coupling effects of SOC, RC, geometric dimensionality, external trap potential and interatomic interaction is presented explicitly in following by means of accurate variational analysis.
Variational analysis for phase transition and stability. In order to obtain the stability diagram and the collapse dynamics of trapped D-dimensional binary BEC with SOC, while centers around the zero momentum and plane wave states, the following normalized spinor order parameter is convenient and applicable 37 , , the variational rate of radius δ, the phase difference between two pseusospin states φ, and the average spin polarization s (−1 ≤ s ≤ 1), i.e., 〈σ z 〉 = s. Upon substitute the spinor wave function into the Lagrangian  Obviously, in the ground state, the BEC is located at a momentum state (k m , 0, 0), where k m = k 0 s. For the equilibrium state, one has δ = 0 and φ = 0, while s and R can be determined by the stationary equation of Eqs (8) and (9). The width of the condensate R is closely related to SOC and RC, thus the collapse can be manipulated by SOC and RC. The distinct phases of the ground state can be distinguished by the stationary solution of Eqs (7-9), while the breather dynamics of the BEC can be described by the evolution of Eqs (7-9).
On the other hand, in order to acquire a stable equilibrium state, there is not a negative eigenvalue in the corresponding Hessian matrix, i.e., and ∂ 2 E/∂R 2 > 0, which results in Hence, the stable equilibrium state can be obtained by the stationary equations of Eqs (8)(9)(10)(11). In an external trap potential, the instability only demonstrates as collapse induced by strong interatomic attractive interaction. Thus, the BEC is unstable for R → 0, i.e., collapse occurs. However, for a stable ground state, R should be finite, where some distinct phases are demonstrated, which depends on s.
Phase and stability diagrams. Here, we focus on distinct phases, the stability boundary and stabilizing mechanism of 3D system, while discuss the difference of phase transition and stability between 2D and 3D system. In 3D system, the ground state can not been obtained analytically and the boundaries of the phase transition and stability can not be depicted analytically, which can only be described numerically by using Eqs (7)(8)(9)(10)(11), however, in particular cases, the analytical conditions of phase transition and stability can been acquired.  9), the stability condition can be obtained as i.e., c = 2π for D = 2 and c = 8π(2π) 1/2 /(5 5/4 ω 1/2 ) = 2.68π for D = 3. Thus, in this case, the phase transition only depends on SOC and RC, while the stability criteria only depends on geometric dimensionality and external trap potential.
Case II: For k 2 0 2 Ω = or Ω = k 0 = 0, the phase transition only depends on g/g 12 (see Eq. (8)). The BEC is in phase I for g 12 < g and in phase II for g 12 > g. Thus the intra-species (inter-species) repulsion promotes the system to convert from phase II (I) into phase I (II). Therefore, the phase transition has nothing to do with SOC, RC, geometric dimensionality and external trap potential. In this case (i.e., Ω = k 0 = 0), by using Eqs (8)(9)(10)(11), the stability condition can be obtained as 12 Thus, the collapse occurs if Eq. (14) is not satisfied, i.e., as g 12 = 0, the collapsing critical intra-species attraction g cri = −c = −2.68π for D = 3 and g cri = −c = −2π for D = 2 37 . This agrees with the previous theoretical 17 and experimental 17,27,28 results for single component BEC (g 12 = 0).
Case III: For a given inter-atomic repulsive interaction, or weak attractive interaction, i.e., g > −c and g 12 > −c, the BEC is always stable, and the condition of phase transition can also been obtained analytically. From the stationary equation of Eq. (8), one can acquire R g g k (( )/(2 )) /(2 ) − Ω for a crossover from s = 0 into s ≠ 0, i.e., from phase I into Phase II. By inserting it into the stationary equation of Eq. (11), the condition of phase transition is obtained as (2 )/ Evidently, the phase structure of BEC can be elaborated by adjusting SOC, RC, interaction and the system's geometric dimensionality. When f > 0, the system enters into phase I, on the contrary, the BEC is in phase II for f < 0, which is also clearly shown in Fig. 2. For a fixed g 12 /g, as the increase (decrease) of RC (SOC), the system transfers from phase II into phase I. When g 12 /g changes, the critical value k ( / ) 0 2 c Ω of phase transition from phase II to phase I is gradually changes, as well as Ω < k ( / ) 2 0 2 c for g 12 /g < 1 and k ( / ) 2 0 2 c Ω > for g 12 /g > 1. Moreover, the phase transition boundary of 2D system is closer to Ω = k / 2 0 2 than that in 3D case, i.e., the phase transition is more sensitive to Ω k / 0 2 for a lower dimension system. Thus, the phase transition can more easily take place in 2D than 3D system. To more clearly understand the nature of this phase transition, E k / 2 0 2 ∂ ∂ versus k 0 is demonstrated in Fig. 3. The second derivative of energy has a jump at k 0 = k 0c , where the system exhibits a phase transition from phase I to phase II, and the jump point k 0c shifts right as the increase of Ω and shifts left as the increase of g 12 /g. These indicate that this phase transition has a second-order nature. The critical values k 0c are satisfied with the condition of phase transition f = 0 for interatomic repulsive interaction or weak attractive interaction, and also in agreement with Fig. 2.    Case IV: For general case in 3D system, i.e., both repulsive and attractive interaction, although the ground state can not been obtained analytically, the complete phase diagram and stability diagram can be presented numerically by variational Eqs (7)(8)(9)(10)(11), which is explicitly demonstrated in g − g 12 plane (Fig. 1) and Ω-k 0 plane (Fig. 4). For strong inter-species attraction ( Fig. 4(a,b)), there is that the collapse occurs as k / c 0 2 1 η Ω > , otherwise the stable phase II exists, where η c1 is determined by the coupling effects of the interatomic interaction, the harmonic trap potential and the system's geometric dimensionality. In this case, the stabilizing mechanisms are illustrated that can lead to the BEC polarized, then stabilize the BEC with strong inter-species attractive interaction, and make the BEC enter into phase II. In other way, the weak RC and strong SOC (i.e., k / c 0 2 1 η Ω < ) can produce an effective inter-species repulsive interaction, which can balance the mean-filed inter-species attractive interaction to prevent the collapse and result in stable phase II. Further, η c1 in 3D system is smaller than that in 2D system, which indicates that the BEC is more unstable in 3D than 2D system. These are also demonstrated in Fig. 1(a2,b2,a3,b3), where the stability and phase boundaries are depicted in g − g 12 plane.
For strong intra-species attraction (Fig. 4(c,d)), the Ω-k 0 plane is divided into three regions, the collapse region , the stable phase II region for η η < Ω < k / c c 1 0 2 2 , and stable phase I for can generate a spin polarized state, which can stabilize the collapse produced by strong intra-species attractive interaction, and make the BEC convert into phase II, whereas k / c 0 2 2 η Ω > can introduce a spin overlapped state, which can also stabilize the collapse produced by strong intra-species attractive interaction, and create a stable phase I. In other words, the strong RC and weak SOC (i.e., η ) produces an effective intra-species repulsive interaction, which can balance the mean-field intra-species attractive interaction then prevent the collapse. Compared with the 3D system, in 2D system, the region of collapse (stable BEC) decreases (increases), and the stability boundary is closer to Ω = k / 2 0 2 , namely, the BEC is more stable in 2D system. These are also evidently depicted in g-g 12 plane (see Fig. 1(a1,b1,a3,b3)). In brief, the strong attractive collapse may be stabilized by elaborating the SOC, RC, harmonic trap potential, geometric dimensionality, and atomic interaction. Moreover, the collapse can be more easily stabilized in 2D than 3D system for all cases.
Collapse dynamics. The breathing behavior of the BEC can be obtained by Eqs (7)(8)(9). Obviously, the breathing behavior is closely related to SOC, RC, and interatomic interaction. While the BEC width declines to zero in However, due to the interference of the two species under the competitive effects between SOC and interatomic interaction, a new unique quantum phase, i.e., the stripe phase, has been predicted theoretically [12][13][14][15] and observed experimentally 6 , which display a periodic density modulation spatially. The region for existing the stripe phase can be qualitatively obtained by similar variational methods where the trial wave function can be assumed as a superposition of two plane wave with momentum ±k x corresponding to Eq. (3), i.e., the stripe phase only exists in the region of g 12 < g and Ω k / 2 0 2  . As shown in Fig. 8, the dynamic evolutions for the density profiles of 3D wave packets in different phase regions are demonstrated. In the parameter region of the stripe phase (see Fig. 8(a)), the evolutions of the density profiles demonstrates as the spatially periodic density modulation along x direction because the SOC is only added in x direction, i.e., there indeed is a stable stripe phase in this region. Otherwise, there is not a stable stripe phase out of this region (see Fig. 8(b)). Thus, the occurrence of the stripe phase can modify the phase diagram of zero-momentum and plane wave phases, whereas their stability criteria can not be changed.

Discussion
In conclusion, the combination of analytical and numerical methods depicts the phase and stability diagrams, while illustrates the stability mechanism of the trapped D-dimensional binary BEC with SOC. The dependence of phase transition and the stability of the system on SOC, RC, interatomic interaction, geometric dimensionality and external trap potential is presented in full parameter space. It is shown that the SOC-induced modification of dispersion relations can generate a spin polarized state or an overlapped state, which can compensate for the strong interatomic attraction, and stabilize the collapse, then changes the stability criteria. In particular, the condition of phase transition and stability criteria are gained, which represents that the phase transition and stability not only rely on interatomic interaction, SOC and RC, but also depend on the system's geometric dimensionality. A full and deep understanding of the dependence of phase transition and stability mechanism on geometric dimensionality and external trap potential is demonstrated. It is also conformed by the collapse dynamics, which indicates that from 2D to 3D system, the mean-field attraction for inducing the collapse is reduced and the collapse speed is enhanced, then the collapse can be more easily stabilized in 2D system. That is, the collapse dynamics may be elaborated experimentally by adjusting the SOC, RC, geometric dimensionality, harmonic trap and atomic interaction.

Methods
Here we have investigated the ground state properties of spin-orbit coupled Bose-Einstein condensate in a harmonic trap potential by using variational approach from the mean field Gross-Pitaevskii equations. The stability phase diagram and stability mechanism are analytically obtained, while collapse dynamics, collapse time and collapse speed are also presented, which both confirmed by the direct numerical simulations of Gross-Pitaevskii equations by means of the fourth-order Runge-Kutta scheme.