Phase Diagram of Solitons in the Polar Phase of a Spin-1 Bose-Einstein Condensate

We theoretically study the structure of a stationary soliton in the polar phase of spin-1 Bose--Einstein condensate in the presence of quadratic Zeeman effect at zero temperature. The phase diagram of such solitons is mapped out by finding the states of minimal soliton energy in the defining range of polar phase. The states are assorted into normal, anti-ferromagnetic, broken-axisymmetry, and ferromagnetic phases according to the number and spin densities in the core. The order of phase transitions between different solitons and the critical behaviour of relevant continuous transitions are proved within the mean-field theory.

Generation of topological defects caused by spontaneous symmetry breaking (SSB) is a universal phenomenon constantly addressed in cosmology, high energy and condensed matter physics [1][2][3].Depending on the type of SSB transition and the degrees of freedom of order parameters topological defects are created in various forms, such as vortices/strings [4,5], domain walls/solitons [6], hedgehogs/monopoles, or their composites.In the past decades, the multi-component or spinful superfluids of liquid 3 He [7,8] and gaseous Bose-Einstein condensates (BECs) [9,10], which induce SSB in numerous ways because of their multiple degrees of freedom of order parameters, have served as a testing ground for the creation of novel topological defects.In these multi-component systems, the core of a defect is not necessarily singular.A simple example is the so-called coreless vortex in which the vortex core of a superfluid component is occupied by other components [11], in contrast to the vortex in a single-component superfluid, whose core is singular due to the nonexistence of superfluid order.Another remarkable example is the core structure of vortices in superfluid 3 He-B [12], where vortices undergo a phase transition by changing the core structure depending on the pressure and temperature.
Recently, Kang et al. [13] observed the wall-vortex composite defects spontaneously generated in a quasi two-dimensional (2D) spinor BEC of sodium quenched from the antiferromagnetic phase to the polar (P) phase [14].In the early stage of evolution, the spontaneous breaking of discrete symmetry causes the formation of domain walls as dark solitons in the |m = 0 Zeeman component.In contrast to the dark soliton in a scalar BEC [15], whose lifetime is generally short, the wall structure exists for a long time with its core occupied by the |m = ±1 Zeeman components, forming a composite defect of wall and half-quantum vortices [13].This suggests that, owing to its solitary nature, soliton not only acquires a role in one-dimensional (1D) systems [6,[16][17][18][19] but can also serve as the building block of some exotic composite defects in the multi-component systems in a higher dimension.Despite that the experiment in Ref. [13] has revealed certain dynamical features of the solitons in a spinor BEC, a theoretical investigation of such solitons in the presence of quadratic Zeeman effect is lacking [20][21][22][23][24].
In this Letter, we investigate the core structure of solitons in the P phase of spin-1 BECs, using the Gross-Pitaevskii formalism.The phase diagram, featuring a variety of stationary states distinguished by the number and spin densities in the soliton core, is obtained and illustrated in Fig. 1 (a).Upon identifying the phase boundaries, we determine the order of phase transition between two bordering phases and further analyze the critical behaviour of those continuous phase transitions by applying a mean-field-theory based Ginzburg-Landaulike approach.Our findings unfold the less-explored aspects of soilton physics in the multi-component superfluid systems.
Formulae-In the mean-field theory [9], a spin-1 BEC can be well described by a vectorial macroscopic wave function, Ψ(r, t) = (Ψ 1 , Ψ 0 , Ψ −1 ) T at very low temperatures.A stationary solution of Ψ corresponds to the local minimum of the thermodynamic energy functional with the energy density given by where M is the atomic mass, p and q the linear and quadratic Zeeman coefficients respectively, µ > 0 the chemical potential, n = the number density, and f = (f x , f y , f z ) T is the spin density vector with and The constants c d and c s denote separately the strengths of density-density and spin-spin interactions mediated by s-wave collisions.In this paper, we assume p = 0 [13], and thus the P phase is now delimited by 2 ≥ q/µ > 0 and c s /c d > −q/2µ > −1 [10].Without loss of generality, the ground state of spin-1 BEC is represented by Ψ P = (0, √ n b , 0), where n b is the bulk density.Note that the arbitrariness of the global phase ensures that Ψ P can be real-valued.The chemical potential µ = c d n b and the density healing length ξ = / √ M µ characterize the energy and length scales of P phase.
Consider a flat domain wall normal to the x-axis in a uniform spinor condensate.The translational invariance in the y-and z-directions renders the problem into a 1D case.Denoting the wave function as Ψ soliton (x) = (ψ 1 , ψ 0 , ψ −1 ) T , and assuming that the soliton is centered at x = 0, then in analogy to the dark soliton in the scalar BEC, ψ 0 obeys the boundary condition ψ 0 (x = ±∞) = ± √ n b , imposing a π-phase jump across the soliton core.We also assume that Ψ soliton produces no currents, i.e., j m = |ψ m | 2 ∇arg(ψ m ) = 0 for the stationary solution [25].Although ψ 0 vanishes at the core center, the ψ ±1 components can readily occupy the core region, forming the so-called bright-dark-bright soliton [22,26,27].Note that E[ Ψ] is invariant under a rotational transformation of f about the z-axis.We thus choose (f x , f y ) = (f ⊥ , 0) for simplicity.Consequently, the fields ψ m (x) can be real-valued.
The soliton solution with lowest energy is determined by minimizing the soliton energy, This quantity corresponds to the excess energy per unit area of the soliton, which plays the role of the tension coefficient of domain wall [28].Requiring that Ψ soliton minimizes the soliton energy, leads to the following coupled GP equations, (3) Using imaginary-time propagation method, the desired solutions are obtained by numerically solving Eq. (3) subjected to the Neumann boundary condition ∂ x ψ m (x = ±L/2) = 0, where the system size L is sufficiently large compared to the width of the core [29].
Phase Diagram-Figure 1 (a) shows the phase diagram of soliton obtained numerically by minimizing the soliton energy.Here, by rescaling energy and length by µ and ξ, the number of independent parameters reduces to two, namely, q/µ and c s /c d .We conclude that the core structure of these soliton states can be distinguished into four types: normal (N), antiferromagnetic (AF), brokenaxisymmetry (BA) and ferromagnetic (F), as shown in Fig. 2. The N-core is prescribed by ψ ±1 (x) = 0, or equivalently, ψ 0 is precisely the dark-soliton solution in the usual scalar BECs [9], such that Ψ soliton of this kind features the "normal (N) state" at its core center as shown in Fig. 2 (a).The AF-core illustrated in Fig. 2 (b) is just the case observed in Ref. [13], which is occupied by the "AF state" with ψ 1 = −ψ −1 , i.e., ψ 1 and ψ −1 are out of phase.Correspondingly, the BA-core is occupied by the "BA state" with ψ 1 = ψ −1 , i.e., ψ 1 and ψ −1 are in phase.We notice that, in contrast to the N-and AF-cores whose spin density vectors f vanish everywhere, the spins of the BA-core are all aligned in the transverse direction with a node at the core center, in other words, f z = 0, as shown in Fig. 2 (c).Finally, the F-core illustrated in Fig. 2 (d) can be considered as a generalization of BA-core, which is occupied by the "F state" with f z = 0 in the core.The F-core features nonzero spin distribution in both transverse and axial directions.
The phase diagram can be understood qualitatively by considering two length scales associated with the thickness of solitons.For simplicity, let us first consider the case of magnetization-free AF-N transition for c s > 0.Then, the relevant length scales in this case are related to the chemical potential and quadratic Zeeman term by ξ = / √ M µ and ξ q = / √ M q, respectively [13].For the AF-core, we have ξ q > ξ in the regime of small q/µ, and the ψ ±1 components are largely accommodated in the core region.With increasing q/µ, ξ q becomes comparably smaller than ξ, and both ψ ±1 components in the AF-core diminish until they totally vanish before transitioning to the N-core.As a result, we conclude that core size∼ max(ξ, ξ q ).
For the ferromagnetic interaction (c s < 0), the above argument on the length scales becomes intricate as the nonzero f of the F-and BA-cores would introduce additional length scale related to the spin-spin interaction.To probe into the essence of the phase transition between different types of soliton, it is necessary to appeal to other more rigorous approaches, and to this end, we introduce the real scalar fields which are defined as Accordingly, the number and spin densities are respectively expressed as We note that there is a jump of Φ 0 ± across the horizontal axis, c s = 0, implying that both of the AF-F and AF-BA phase transitions are discontinuous or first order.On the other hand, Φ 0 ± varies continuously along the horizontal-axis: the AF-N, F-BA and BA-N phase transitions are continuous or second order.We investigate the order of phase transition by computing the derivative of the soliton energy with respect to q/µ and c s /c d , and the conclusion is consistent with the above observation [30].
Perturbative approach-To describe the phase boundaries of the N-core region in the phase diagram quantitatively, we introduce the equations of Φ ± , obtained from Eq. (3), In the close vicinity to the left of the N-core region in the phase diagram, Φ ± are perturbatively small in magnitude and we can assume ψ 0 = √ n b tanh(x/ξ).Similar to the case of coreless vortices in the segregated binary condensates [31], we neglect the cubic terms in the righthand-side of Eq. ( 5), and the core-occupying components can be approached by the exact solution of Pöschl-Teller equation [32,33], where E ± = −q − (1 ± 1)(c s /c d )µ are the eigenvalues determined by the potentials The above perturbative treatment suggests that Φ ± respectively correspond to the bound-state solutions of lowest energy eigenvalues ±,0 associated with V ± in Eq. ( 6), i.e., where [34].The phase boundaries of the AF-N and BA-N transitions are determined by E ± = ±,0 and thus we obtain equations for the phase boundaries, .
(8) Remarkably, the above equations coincide closely with our numerical results in Fig. 1.In particular, the condition c s /c d > −1/2 is coincidentally satisfied since the BA-N boundary terminates at (q/µ, c s /c d ) = (1, −1/2) on the boundary of the bulk BA phase.
Flat-core limit-To gain an insight into the smallq regime, where ψ ±1 are barely suppressed in the core, in the phase diagram, we introduce a variational analysis assuming the flat-core limit n = const.In the limit q → 0 + , the soliton solution can be well approximated by the bright-dark-bright soliton solution in multicomponent condensates [22,26,27] Here φ v (0 ≤ φ v ≤ π) and ξ v are variational parameters, and the total density is normalized to the constant bulk density n b .Substituting Ψ v into Eq.( 2), the soliton energy is explicitly expressed by where (10) with respect to ξ v and φ v , we obtain the following energetically favorable states: (i) AF-core with (φ v , ξ v ) = (−π/4, ±ξ µ/2q) , which is a local minimum (saddle point) for c s > 0 (c s < 0).Note that these states are all doubly degenerate [35] .
The above variational analysis theoretically justifies the small-q regime of the phase diagram, and, in particular, accounts for the existence of the F-core soliton.The analysis also describes q/µ− and c s /c d -dependence of the soliton thickness ξ v ; the thickness could depend on the spin-spin interaction through c s /c d in addition to the quadratic Zeeman effect for c s < 0 as mentioned before.
Critical behaviour of the observables-The solitons could undergo a continuous phase transition from AF-to N-core, or successive continuous phase transitions from F-to BA-and then to N-core. Figure 1 (b) and (c) reveal that these phase transitions are also accompanied by the vanishing of Φ 0 + or Φ 0 − at a critical point q = q c , which suggests that Φ 0 ± can serve as the effective order parameters for the continuous transition.In what follows, we employ a mean-field analysis similar to the Ginzburg-Landau theory [28] to probe the critical behaviour of the continuous phase transition between different core structures.
We expand the soliton energy in terms of the effective order parameter ϕ near the critical point (q ≤ q c ) as . However, such power-law behavior of spin densities does not arise in the AF-N transition, as f vanishes for both AF-and N-cores.
The analysis for the F-BA transition is delicate.As shown in the inset of Fig. 3 (c), both Φ 0 ± 2 vary virtually linearly in the critical region.We notice that while Φ 0 − 2 vanishes at q c , Φ 0 + 2 decreases to a nonzero value at q c .Thus we assume Φ 0 + 2 = B 0 + B 1 |q − q c | for q ≤ q c , where B 0 , B 1 > 0, and the central density is given by n 0 = B 0 + (β + B 1 ) |q − q c |.Here the coefficients β, B 0 and B 1 can be determined from the linear fitting of n 0 .Furthermore, it is straightforward to show that f max ), which nicely fit to the numerical simulations of f max ⊥ and f max z , as demonstrated in Fig. 3 (c).We notice that, in the close vicinity of q c , the power-law behavior f max z ∝ |q − q c | is resumed.
Concluding remarks-We have theoretically inquired into the core structure of stationary solitons in the P phase of spin-1 BECs.In the presence of quadratic Zeeman effect, four different types of structure are identified according to the number, spin densities in the core.Our theoretical analyses support well the numerically obtained phase diagram of solitons.The critical behaviour of continuous phase transition between different types of solitons is predicted by introducing the Ginzburg-Landau-like mean-field theory.
A variety of types of solitons found in this study can be in principle generated via quantum quenches within the reach of the state-of-the-art techniques [13,36] or phase imprinting method [37].Our study thus offers possibilities to explore different aspects of soliton physics, e.g., magnetic phase transition and spontaneous symmetry breaking of soliton cores.One of the follow-up studies would be the phase diagram of solitons in other experimentally accessible phases of spinor BECs.Dynamical properties of solitons, such as the instability or collisions between solitons, are also important issues to be looked into while exploring the non-equilibrium dynamics in multi-component superfluids.
Acknowledgement I.-K.Liu and S.-C.Gou acknowledge the financial support from the Taiwan MOST 103-2112-M-018-002-MY3 grant.H. Takeuchi is supported by the Japan Society for the Promotion of Science (KAK-ENHI Grant No. JP17K05549), and in part by the OCU "Think globally, act locally" Research Grant for Young Scientists 2019 through the hometown donation fund of Osaka City.We appreciate the fruitful discussions with Yong-il Shin, Yu-Ju Lin and Thomas Bland.

A. Numerical analysis
For a homogeneous spin-1 condensate in the P phase, we choose the chemical potential µ = c d n b (n b is a constant density) [10] as the characteristic energy.Accordingly, the characteristic length and time scales are defined as ξ = / √ M µ (the healing length) and τ = /µ, respectively.By rescaling the energy functional, Eq. ( 1), with respect to these variables, one obtains the dimensionless energy functional in the absence of linear Zeeman term, The scaled wave function is related to the original one by Ψ m = Ψ m / √ n b .Accordingly, the dimensionless total number density is given by n = , and the scaled spin-density vector are given by is . Now we consider a soliton solution by chossing xaxis as the coordinate normal to the wall with the wave function, Ψ solition = (ψ 1 (x ), ψ 0 (x ), ψ −1 (x )) T .Consequently, the dimensionless soliton energy is where E bulk = L/2ξ = L /2.The corresponding dimensionless GP equations are derived by using Hartree variational principle, δα /δψ m = δ[ Ψ solition ]/δψ m = i∂ t ψ m , that are explicitly expressed as, A3) where t = t/τ is the dimensionless time.To find the soliton solution Ψ soliton minimizing the soliton energy, we employ the imaginary-time propagation method to solve for the lowest state of Eq. (A3).We select initial conditions of the form, where A 1 and A −1 are complex numbers, and, owing to the rotational symmetry about z axis in spin space, one can choose θ = 0 without loss of generality.Furthermore, in order to carry out the imaginary-time propagation with a faster convergence, we impose the following conditions (1) and A 1 and A −1 are in phase; (4) A 1 = A −1 = 0 such we can arrive at the designated AF-, BA-, F-or N-core more efficiently.Note that if A 1 and A −1 are neither in phase nor out of phase, it turns out that the imaginary-time evolution of the initial states will end up in a soliton-free bulk solution in the P phase because such states are not spin-conserving.
In the simulations, we consider a 1D mesh where the grid size is dx = 0.25.The length of the mesh ranges from L = 80 to 720 to guarantee the boundary conditions ψ ±1 (x = ±L/2) = 0, ψ 0 (x = ±L/2) = ±1 or ∓1, as well as the Neumann boundary condition, ∂ x ψ m (x = ±L/2) = 0 are fulfilled.Generally speaking, L = 80 is sufficiently large for q ≥ 0.0375 as the largest length scale of the core size is ξ v = 1/ √ 2q ∼ 6.3 for q = 0.0125 in the flat-core approximation.The convergence of the imaginary-time propagation is controled by the magnitude of the local error produced during the imaginary-time evolution governed by Eq. ( A3) based on the criterion max(δµ ) < 2.5 × 10 −13 .
The order of phase transition between two different types of soliton cores is further examined by computing the derivatives ∂α /∂q and ∂α /∂c s numerically.In Fig. 4 (a), a smooth variation of ∂α /∂q is clearly observed, whereas ∂α /∂c s exhibits a discontinuity at c s /c d = 0 as depicted in Fig. 4 (b).With a fixed density, the results suggest that the AF-N, F-BA and BA-N transitions are continuous or second order, and the AF-F and AF-BA transitions are discontinuous or first order.
Replacing the Pöschl-Teller potential in Eq. (B1) with the effective potentials V ± (x) = −U ± sech 2 (x/ξ), we have such that the corresponding λ parameters are given by Consequently, the bound-state solutions exist only when the following conditions are met: c s /c d > −1/2 for g > 0, and (1 ± 1)(c s /c d ) > −8 for λ ± > 1. Denote ±,n as the n-th energy eigenvalues for the potential V ± (x) respectively.For V − , it follows that λ − = 2 > 1, and according to Eq. (B2), there is only one bound-state Φ − with energy level −,0 = −µ/2 exists.Thus Eq. ( 6) has non-trivial solution if and only if E − = −,0 , implying that Next, we solve the eigenvalue equation of Φ + in the regime, −q/2µ < c s /c d < 0, where the BA-core is favorable.This yields 1 < λ + < 2, implying that only one bound-state solution exists, which has the energy level Likewise, Eq. ( 6) has non-trivial solution if and only if E + = +,0 , which then gives the equation of the BA-N phase boundary Equations ( B5) and (B7) represent the boundaries of the AF-N and BA-N phase transitions respectively, which are mutually exclusive.This rules out the possibility of the F-N transition as none of Φ ± vanish in the F-core phase.Since Eq. (B7) are well-founded in the regime −q/2µ < c s /c d < 0 and intersects with Eq. (B5) at (q/µ, c s /c d ) = (1/2, 0), it follows that the whole border separating the N-core phase from other possible phases can be mapped out by piecewisely joining the lower part (c s < 0) of Eq. (B7) to the upper part (c s > 0) of Eq. (B5) at the point (q/µ, c s /c d ) = (1/2, 0) , as shown in Fig. 1.

C. The Ginzburg-Landau-like approach
In the framework of mean-field theory, the secondorder phase transition, which is characterized by the vanishing of order parameter at the critical point, can be described by the Ginzburg-Landau (GL) theory.In view of the continuous nature of the AF-N,F-BA and BA-N phase transitions, by analogy with the GL theory, the soliton energy is expanded in terms of the order parameter in the critical region for q ≤ q c by α where ϕ is the order parameter, and specifically, ϕ = Φ 0 + for the BA-N and ϕ = Φ 0 − for the AF-N and F-BA transitions.The coefficients α i (i = 0, 1, 2) in Eq. (C1) are all positive and determined by Φ ± and Φ 0 .For the BA-N and AF-N transitions, α 0 can be exactly determined and is explicitly given by α 0 = 4µn b ξ/3.Minimization of α with respect to ϕ yields with β = α 1 /α 2 .Substituting Eq. (C2) into Eq.(C1), we get In Fig. 5 (a) and (b), the calculated α c are well fitted by the quadratic function of Eq. (C3) for AF-N and BA-N transitions, with α 1 = 0.56 and 0.57 respectively.

AF-N transition
In this case, ϕ = Φ 0 − , and the two parameters β = 2n b /µ and q c = µ/2 are exactly determined based on the results obtained by Pöschl-Teller approach.Since |ϕ| in other words, the central density linearly decreases to zero in the critical region for q ≤ q c , which is clearly illustrated in Fig. 3 (a).

BA-N transition
In this case, ϕ = Φ 0 + .However, as α i become c sdependent due to the spin-spin interaction, β and q c can only be determined numerically.For example, given c s = −0.0125cd , we find β = 1.97n b /µ and q c = 0.51µ.Like the AF-N phase transition, the central density exhibits linear decline near the critical point, n 0 = |ϕ| 2 = Φ 0 + 2 ∝ |q − q c | for q ≤ q c as shown in Fig. 3 (b).In what follows, we shall show that, in addition to the central density, the maximal spin density also exhibit power law behavior near the critical point.Recall that for the BA-core, f z (x) = 0 and f ⊥ (x) = 2Φ 0 Φ + .Defining max [f ⊥ (x)] = f max ⊥ and from Fig. 2, we see that f max ⊥ always occurs around the core.Assuming that f max ⊥ occurs at x = x f in the proximity of the core center x = 0, and expanding f max ⊥ to the first order, we obtain The above behavior is numerically verified and is shown in Fig. 3 (b) in good agreement.

F-BA transition
In this case, ϕ = Φ 0 − , and the two parameters β 1 and q c are to be determined numerically.As shown in the inset of Fig. 3 (c), in the F-BA transition, |Φ 0 − | 2 declines linearly in the critical region and becomes zero at q c as expected, whereas Φ 0 + linearly decreases to a nonzero value at q c .Thus, in addition to the decline of the order parameter, Φ 0 The coefficients β, B 0 and B 1 are determined by the linear fit to n 0 .Given c s = −0.0125cd , we have q c = 0.25µ, β = 1.04n b /µ, B 0 = 0.5 and B 1 = 0.99.Likewise, we define max [f z (x)] = f max z , and from Fig. 2, it follows that f max z = f z (x = 0).Consequently, in the critical region for q ≤ q c , we get and from Eq. (C4), it follows that

FIG. 1 .
FIG. 1.(a) The typical phase diagram of stationary soliton states in the defining range of polar phase (p = 0).(b) and (c) Variations of the order parameters Φ 0 ± = Φ±(x = 0).The phase boundary drawn in black solid line [AF-F and AF-BA in (a)] indicates a phase transition of first order, across which the order parameter changes discontinuously.On the other hand, the phase boundary drawn in broken line [AF-N, F-BA and BA-N in (a)] indicates a phase transition of second order, across which the order parameter changes continuously.

2 ∝
2).Here ϕ = Φ 0 − for the AF-N and F-BA transition, and ϕ = Φ 0 + for the BA-N transition.Minimization of α c with respect to ϕ yields |ϕ| 2 = β |q − q c | with β = α 1 /α 2 , which manifests itself as a linear ramping down of the central particle density,n 0 = n(x = 0) = |ϕ| 2 = Φ 0 ± |q − q c |,in the critical regions of the AF-N and BA-N transitions, respectively, as demonstrated in Fig. 3 (a)-(b).The parameter β can be determined by linear fitting of n 0 .Likewise, the maximal spin density around the core demonstrates a power law in the critical region.As is derived in the Sec.C ofsupplemental material for the BA-N transition, f max ⊥

FIG. 5 .− 2
FIG. 5.The αc (black dots) as a function of q/µ is shown for (a) AF-N; (b) BA-N and (c) F-AB transitions.The grey solid lines in (a) and (b) indicate the curve fitting using Eq.(C3).In (d), the αc as a function of Φ 0 − 2 is illustrated, which