Singularity in the matrix of the coupled Gross-Pitaevskii equations and the related state-transitions in three-species condensates

An approach is proposed to solve the coupled Gross-Pitaevskii equations (CGP) of the 3-species BEC in an analytical way under the Thomas-Fermi approximation (TFA). It was found that, when the strength of a kind of interaction increases and crosses over a critical value, a specific type of state-transition will occur and will cause a jump in the total energy. Due to the jump, the energy of the lowest symmetric state becomes considerably higher. This leaves a particular opportunity for the lowest asymmetric state to replace the symmetric states as the ground state. It was further found that the critical values are related to the singularity of either the matrix or a sub-matrix of the CGP. These critical values are not arising from the TFA but inherent in the CGP, and they can be analytically expressed. Furthermore, a model (in which two kinds of atoms separated from each other asymmetrically) has been proposed for the evaluation of the energy of the lowest asymmetric state. With this model the emergence of the asymmetric ground state is numerically confirmed under the TFA. The theoretical formalism of this paper is quite general and can be generalized for BEC with more than three species.

Furthermore, based on the analytical formalism and the singularity of the equations, effort is made to divide the whole parameter-space into zones, each supports a specific spatial configuration. This provides a primary frame for plotting the phase-diagrams in the future. Besides, a model for calculating the total energies of asymmetric states has also been proposed. The possibility of the emergence of asymmetric g.s. has been primarily evaluated.
The theoretical formalism of this paper is quite general and can be generalized for the K-BEC with K larger than 3.

Hamiltonian and the coupled Gross-pitaevskii equations
We assume that the 3-BEC contains N S S-atoms with mass m S and interacting via (1) 2 . H B and H C are similarly defined. We assume that no spatial excitations are involved in the g.s. Thus, each kind of atoms are fully condensed into a state which is most advantageous for binding. Accordingly, the total many-body wave function of the g.s. can be written as where ε A is the chemical potential. Via cyclic permutations of the three indexes (A, B, C) together with (u 1 , u 2 , u 3 ).
From eq. (3) we obtain the other two for u 2 and u 3 . It is emphasized that the normalization ∫ = u dr 1 l 2 (l = 1, 2 and 3) is required.

Formal solutions under the Thomas-Fermi approximation
Since N A , N B and N C are considered to be large, TFA has been adopted. The applicability of this approximation has been evaluated via a numerical approach given in refs 22 and 35 and will be discussed later. Under the TFA, the CGP become , and ′ d l l is the algebraic cominor of α l′l . (ii) Form II: When one and only one of the wave functions is zero inside a domain (say, u n /r = 0), the other two must have the unique form as where l, m, and n are in a cyclic permutation of 1-2-3 (the same in the follows), When the Form II has u n /r = 0, a more precise notation Form II n is adopted for the detail. (iii) Form I: When one and only one of the wave functions is nonzero in a domain (say, u l /r ≠ 0), it must have the unique form as l ll l 2 2 2 Obviously, u l /r in this form must descend with r. For the case u l /r ≠ 0, the more precise notation Form I l is adopted. (iv) Form 0: In this form all the three wave functions are zero.
If a wave function (say, u l /r) is nonzero in a domain but becomes zero when r = r o , then a downward form-transition (say, from Form III to II) will occur at r o . Whereas if u l /r is zero in a domain but becomes nonzero when r = r o , then a upward form-transition (say, from Form II to III) will occur at r o . r o is named a form-transition-point, and it appears as the boundary separating the two connected domains. In this way the formal solutions serve as the building blocks, and they will link up continuously to form an entire solution. They must be continuous at the form-transition-points because the wave functions satisfy exactly the same set of nonlinear equations at those points. However, their derivatives are in general not continuous at the boundaries.

An approach for obtaining analytical solutions of the CGP
In this section we consider the case that all the parameters are given and the values of the three {ε l } have been presumed. In this case all the formal solutions are known. We will propose an approach to link up the formal solutions to form a chain as a candidate of an entire solution. To this aim we first introduce a number of features related to the linking.
(i) For Form I to III, at least one of the wave function is descending with r.
The proof of this feature is referred to ref. 36, where it is proved that at least one of Y l (or Y l n ( ) for a given n) is positive. This feature implies that, when r increases, the occurrence of a downward form-transition is inevitable, unless a upward form-transition takes place prior to the downward transition. In any cases a formal solution must transform to another form somewhere (except Form 0).
(ii) For a formal solution existing in a domain, the right boundary of the domain and the successor (the successive formal solution) in the next domain have been prescribed when the three {ε l } have been presumed.
To prove this feature, as an example, we assume that u 1 /r and u 3 /r are nonzero in a domain while u 2 /r is zero. This assumption implies that we have assumed − ≥ X Y r 0 when r is given inside the domain (refer to eq. (10)). We define = r X Y / (2) or ≤ 0). Similarly, we (2) or ≤ 0), and = r X Y / 2 2 2 2 or ∞ (if both X 2 and Y 2 are negative or otherwise). Then, the smallest one among r 1 , r 2 , and r 3 is just the right boundary of the domain. Say, if r 1 is the smallest, then u 1 /r → 0 when r → r 1 , and the successor will have the Form I 3 . If r 2 is the smallest, then u 2 /r will emerge at r 2 , and the successor will have the Form III, and so on. Since r 1 , r 2 , and r 3 are prescribed, the right boundary and the successor are prescribed (iii) Once the formal solution in the first domain (starting from r = 0) is prescribed, the formal solutions will link up one-by-one to form a chain in a unique way. There are seven types of formal solutions (say, in Form II 2 , or in Form I 3 , and so on). Each type can appear in a chain at most once.
Obviously, since the successor in each step of linking is uniquely prescribed, the whole chain is prescribed. Since the right boundary of a type is prescribed, the type can not appear twice. (iv) Once a formal solution in a chain is in Form 0, the chain will end. This is because no wave functions can emerge from an empty domain. Otherwise, if u 1 /r emerges alone, it must have the form as eq. (12). This form prohibits the uprising of u 1 /r. If u 1 /r and u 2 /r emerge at the same place, Y 1 (3) and Y 2 (3) must both be negative. This violates the feature (i). If all the {u l /r} emerge at the same place, all the three {Y l } must be negative. This violates also the feature (i).
Based on the above features, we propose an approach as follows: First, we design a chain for a type of entire solutions denoted as, for an example, II 2 -III-II 1 -I 3 (it implies that the first domain has a Form II 2 , the next domain has a Form III, the third domain has a Form II 1 , while the last domain has a Form I 3 ). The prescription on the linking appears as a number of requirements (inequalities) imposing on the W-strengths and the presumed {ε l }. When all the α ′ { } ll and the {ε l } are given inside a specific scope, all the requirements can be met and the designed chain as a candidate can be achieved. At this stage the normalization has not yet been considered. When the three are further introduced, not only the scope but the values of the set {ε l } can be fixed. Then, the candidate will be a realistic entire solution of the CGP. In general, the three equations can uniquely determine the three unknowns {ε l }, unless the design itself is not reasonable. Thus, when the parameters are given in a reasonable scope, we can uniquely find out a realistic entire solution, which is a chain of formal solutions with a specific linking. A detailed practice of this approach for miscible states is given in ref. 36.  , the matrix-of-equations becomes singular and accordingly = D 0. When c AB is close to and crosses over this critical value (from (c) to (d)), {Y l } will suddenly change their signs. It implies a down-falling wave function suddenly becomes up-rising, and accordingly the whole pattern is changed greatly. This is definitely associated with a state-transition in which all the A-atoms suddenly jump from a core to a shell, while all the B-atoms jump in a reverse way as clearly shown in (c) and (d). Accompanying the great change, a remarkable increase in the total energy is expected (this expectation is confirmed below).

State-transition and the singularity of the matrix
From the equality = D 0, it is straight forward to obtain Similarly, p CA and p AB can be defined by permuting the indexes.
Note that: (i) When p BC p CA < 0, c AB crit (3) does not exist (i.e., the matrix will not become singular). Therefore, even a Form III is contained in a chain, the variation of c AB does not assure the occurrence of the state-transition. Only if the other five strengths are so chosen that p BC p CA ≥ 0, the critical point could exist and the transition could occur. (ii) c AB crit (3) deviates remarkably from the well known critical value = c c c for 2-BEC. Thus, the state-transition is remarkably affected by the influence of the third kind of atoms. However, if both c BC → 0 and c CA → 0 (i.e., the influence is removed), one can prove from eq. .
(iii) c AB crit (3) depends on the other five strengths but not on the particle numbers, trap frequencies, and masses. This feature is the same as what has found in 2-BEC. Thus, in an experiment, the variation of the parameters other than the strengths will not change the critical values. (iv) The state-transitions caused by the variation of other strengths can be similarly deduced. For an example, for the intra-species interaction of the A-atoms, when c A increases and arrives at a critical value the matrix will become singular and the transition will occur.
In summary, for an entire solution contains a Form III, when the variation of the strengths leads to a cross-over of the singular point of the matrix, a state-transition will occur. Since the singularity of the matrix is inherent in the CGP but not a product of the TFA, thus the occurrence of the state-transitions at the critical values holds beyond the TFA. In fact, in the earliest study of the 2-BEC, the instability in the neighborhood of the critical value = c c c AB crit A B (2) (the singular point of the two-rank matrix) has been pointed out 8 .
(2) State-transition occurring at the singular point of a sub-matrix-of-equations When an entire solution contains a Form II, another type of state-transition might occur. In Fig. 2, the entire solution is III-II 3 -I 2 in (a) to (c), and is III-II 3 -I 1 in (d) to (f), where the second domain has the Form II 3  , which is the singular point of the sub-matrix of the equations for u 1 /r and u 2 /r only. When c AB is close to this value ((c) and (d)) the second domain becomes very narrow, and the two wave functions become very steep. During the cross-over, Y 1 (3) and Y 2 (3) change their signs and a transition occurs as shown in (c) and (d). Nonetheless, different from the one found in Fig. 1, only a part of the A-and B-atoms are actively taking part in this transition, namely, a part of A-atoms rush out from the core and form a shell, while a part of outward B-atoms rush from the shell into the core. Thus, the corresponding change in spatial configuration is relatively milder. The change appears essentially in the second and the third domains where the C-atoms are absent. Accordingly, the critical value is not at all affected by the C-atoms and is identical to the value of  Fig. 2, c AB crit (3) does not exist in this case due to p BC p CA < 0.

Total energy of symmetric states and the great jump
When the total energy of the lowest symmetric state E tot is higher than the total energy of the lowest asymmetric state E tot asym , the g.s. will be asymmetric. Thus, > E E tot t ot asym is the discriminant to judge whether the g.s. is asymmetric.
When the wave functions are known we can obtain the total energy as (the kinetic energy has been omitted) . Note that, in Fig. 3(a), the crossing over c AB crit (2) does not cause an effect because the associated transition could occur only if the chain contains the building block II 3 , this building block is absent in Fig. 1. While in Fig. 3 (3) does not exist because p BC p CA < 0 as mentioned. Figure 3 confirms that the state-transition has caused a great change in E tot . For the transition shown in Fig. 1(c,d), when the B-atoms rush in, E 2 will increase (because a more compact distribution leads to the increase of the factor ∫ u r r dr ( / ) 2 4 2 ) while P 2 will remarkably decrease. The decrease over takes the increase. We found that, for each B-atom, (E 2 + P 2 )/N B decreases from 1.80 to 1.41. On the other hand, for each A-atom, (E 1 + P 1 )/N A increases from 1.08 to 1.50. Since

Asymmetric states and their total energy
We know from the study of the 2-BEC 9, 11, 12 that, when V AB is sufficiently strong, the A-and B-atoms might give up the symmetry of the trap for lowering the g.s. energy. Therefore, we propose a model where only the distributions of the A-and B-atoms are asymmetric, while the C-atoms are symmetric. Let O denotes the center of the trap. Let a sphere with radius R AB centered at O be divided into two parts by a plane perpendicular to the Z-axis. The plane intersects the Z-axis at z = z 0 (−R AB < z 0 < R AB ). Let the A-atoms be distributed in the lower part of the sphere, and the B-atoms in the upper part. Let the C-atoms be symmetrically distributed in another sphere with radius R C and centered also at O. Then, we assume is for the normalization.
When the values of R AB , R C , and z 0 are assumed, from eqs (16, 17, 18 and 19), the total energy for the asymmetric state (with the kinetic energies neglected), denoted as E tot asym , can be obtained. The parameters R AB , R C , and z 0 are considered as variable. Eventually, they fixed at the values that lead to the minimum of E tot asym . The E tot asym obtained via such a variational procedure is in general higher than its actual value. Thus, in any cases, if we found < E E tot asym tot , the asymmetric state will inevitably replace the symmetric g.s.
The comparison of the two energies are shown in Fig. 3, where (a) and (b) are associated with Figs 1 and 2, respectively. Figure 3a demonstrates clearly that E tot asym is remarkably lower than E tot when > c c AB AB crit (3) . Thus the jump provides a good opportunity for the lowest asymmetric state to replace the symmetric g.s. Whereas when < c c AB AB crit (3) , although E tot asym is remarkably higher than E tot as shown in the figure and E tot will decrease further with the decrease of c AB , we can only say that the g.s. is very probable to be symmetric. This is a point to be further studied.

Division of the parameter-space
If the whole parameter-space Σ, in which a point is associated with a set of parameters, can be divided into zones each supports a specific configuration, various phase-diagrams could be plotted. Thereby the essential features of the system and the effects of the parameters can be visualized. Due to having so many parameters, the phase-diagrams of a 3-BEC would be very complicated. At this moment we are not able to plot them. The following is a primary attempt along this line.
There are four well defined and important surfaces in Σ. They are expressed via the equations = D 0 and = d 0 1 to 3). In other words, each surface is an aggregation of a kind of singular points. We have proved under the TFA that a crossing over these surfaces may cause a state-transition and accordingly an increase of E tot . When the TFA is removed, in a domain of r in which all the {u i /r} are nonzero, the exact CGP can be written as  where ″ u i is the second-order derivative of u i against r. The appearance of the common factor D 1/ at the right side implies that the left-side (namely, the wave functions) is extremely sensitive against the parameters when they are given in the neighborhood of the surface = D 0. This is an important feature of the CGP. When a point in Σ crosses over = D 0, the factor D 1/ changes from ∞ to ±∞. Therefore, the entire solutions (if it contains a Form III) will undergo a dramatic change, and the occurrence of the state-transition (found before under the TFA, refer to Fig. 1) is inevitable Thus, this kind of transition is inherent in the CGP. For the kind of entire solutions containing a Form III, once the variation of the parameters leads to a crossing over the surface = D 0, the transition (denoted as trans-III) happens definitely.
Similarly, in a domain of r in which u n /r = 0, the exact CGP can be written in a form in which both (u l /r) 2 and (u m /r) 2 are proportional to a common factor d 1/ nn . Thus, for the type of entire solutions containing a Form II n , the crossing over the surface = d 0 nn will also lead to a great change in u l /r and u m /r, and accordingly another kind of state-transition (denoted as trans-IIn) occurs as shown in Fig. 2.
Let us define a subspace Σ III as follows. When a set of parameters leads to an entire solution containing a Form III, then the associated point belongs to Σ III , otherwise belongs to its complement. Let the part of the surface = D 0 located inside Σ III be denoted as σ III . Then, σ III appears as a boundary, the crossing over this boundary leads to the trans-III. Similarly, let Σ II 3 denotes the subspace containing the points each leads to an entire solution containing a Form II 3 . Let the part of the surface = d 0 nn located inside the subspace Σ II 3 be denoted as σ II 3 . Then, the crossing over σ II 3 leads to the transition trans-II 3 . We can further define σ II 1 and σ II 2 in a similar way. These four surfaces (σ III and the three σ II i ) together form the boundaries and provide a primitive division of Σ. At the two sides of each boundary, the entire solutions are greatly different.
Nonetheless, these boundaries are not the actual boundaries for the phase-diagrams of the g.s. The latter can be made certain only if exact calculations on both the lowest symmetric and asymmetric states have been performed. However, since the crossing over the above boundaries leads to an increase of E tot and the increase may be large (as shown in Fig. 3). Thus the increase provides an excellent opportunity for the lowest asymmetric state to replace the lowest symmetric state and become the g.s. Therefore, we believe that the exact boundaries for the phase diagrams would partially overlap the boundaries from singularity.

Final remarks
(1) A general approach is proposed to solve the CGP for 3-BEC in an analytical way. TFA has been adopted.
The essence of this approach is to find out the building blocks, i.e., the formal solutions, and the rules for their linking. The entire solutions of the CGP appear as a chain of them. This approach is applicable for obtaining solutions with their chains in various types, and can be generalized to K-BEC with K larger than three. For examples, in a domain where all the K {u l /r} are nonzero, the formal solution has exactly the same expressions as shown in eqs (5, 6 and 7) except that the related matrixes are K-rank. (2) The main result of this paper is the finding of the state-transitions caused by the singularity of the (sub) matrix-of-equation and the associated increase of E tot during the transition. The singularity is not a by-product of the TFA, but an important feature inherent in the CGP. Note that the critical behavior of the multiband superconductors was found to be substantially affected by the interband coupling [31][32][33] . Similarly, the critical point for the state-transition found in this paper differs remarkably from the one of the 2-BEC (refer to Fig. 1) due to the inter-species coupling. Note that the 3-BEC contains three subsystems, each contains two species. Similar to the hidden criticality found also in multiband superconductivity 32 , the critical points of the three subsystems appear as the hidden critical points of the 3-BEC. Under specific conditions state-transitions will also occur at these hidden points (refer to Fig. 2).
(3) A model for asymmetric states has been proposed. Via numerical calculations on some examples, it is demonstrated that the lowest asymmetric state replaces the symmetric states and become the g.s. when the strength of an inter-species interaction arrives at and exceeds its critical value. (4) The whole parameter-space is primitively divided into zones separated by four surfaces as boundaries, each is an aggregation of a kind of singular points. The spatial configurations at the two sides of a boundary are greatly different due to the state-transition occurring during the crossing over the boundaries. The transition is accompanied with an energy increase, the amount of increase might be very large. Thus the state-transition provides an excellent opportunity for the emergence of the asymmetric g.s. Therefore, it is expected that the exact boundaries designating the zones of asymmetric g.s. would overlap partially with the boundaries arising from the singularity. This remains to be confirmed.