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.

Introduction

Accompanying the progress in techniques, the research into the 2-species Bose-Einstein condensates (2-BEC) is gaining an increasing attention in recent years in both the experimental aspect[1-7] and theoretical aspect[8-26]. Many distinguished features have been found, say, the existence of various phases, the critical value in the inter-species interaction and the related instability, the emergence of asymmetric ground state (g.s.)[9, 11, 12], the appearance of vortex, and so on. The study for BEC with three species (3-BEC) has also started[27-30]. It is very interesting to see how the phenomena found in 2-BEC would recover in 3-BEC and whether new phenomena would emerge. Note that, for multiband superconductivity, the interband couplings among a set of different band condensates are important to the critical behavior of the system. Critical temperatures are thereby substantially affected (not determined alone by the Cooper-pair amplitude of a single band)[31-33]. Thus, it is reasonable to expect that the critical phenomena found in 2-BEC (say, a state-transition has been found to take place when the strength of the inter-species interaction arrives at a critical value[34]) might also be affected and new critical phenomena might emerge. Since the BEC with more than two species are experimentally achievable, it is meaningful to perform theoretical research at this stage.

This paper is dedicated to a primary theoretical study on the 3-BEC based on the coupled Gross-Pitaevskii equations (CGP). Although exact numerical solutions of the CGP are very valuable, it is not easy to extract the underlying physics simply via numerical results. In order to gain more insight into the physin class="Chemical">cs, it is helpful to obtain approximate analytical solutions. Therefore, the Thomas-Fermi approximation (TFA) has been adopted. Under the TFA, we provide an approach for obtaining analytical solutions. Thereby the wave functions and the total energies can be obtained in an analytical form, and these quantities can relate directly to the parameters involved. This facilitates greatly related physical analysis. We found that the singularity of the (sub)matrix-of-equations is crucial to the behavior of the BEC. Specific state-transitions will be induced when the parameters vary and cross over a singular point of the matrix. This will be studied in detail below.

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-atoms with mass m and interacting via , (S = A, B and C). The particle numbers are assumed to be huge (say, ≥10000). The inter-species interactions are with the strength c . These atoms are confined by the isotropic harmonic traps We introduce a mass m and a frequency ω. Then, ħω and are used as units for energy and length. The spin-degrees of freedom are assumed to be frozen. The total Hamiltonian iswhere γ  = (m /m )(ω /ω)2. H and H 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 aswhere u 1, u 2 and u 3 are for the A-, B-, and C-atoms, respectively.

In the set of the CGP, the one for u 1 iswhere ε 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 (l = 1, 2 and 3) is required.

Formal solutions under the Thomas-Fermi approximation

Since N , N and N 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 n class="Chemical">TFA, the CGP becomewhere α 11 = N c /(4πγ ), α 22 = N c /(4πγ ), α 33 = N c /(4πγ ), α 12 = N c /(4πγ ), α 21 = N c /(4πγ ), α 13 = N c /(4πγ ), α 31 = N c /(4πγ ), α 23 = N c /(4πγ ), α 32 = N c /(4πγ ), they are called the weighted strengths (W-strengths) and they are related as α 12 α 23 α 31 = α 21 α 32 α 13. ε 1 = ε /γ , ε 2 = ε /γ , ε 3 = ε /γ , they are the weighted energies for a single particle. Recall that there are originally 15 parameters (N , m , ω , c , c ). Their combined effects are fully represented by the nine α (only eight of them are independent). Thus, based on the W-strengths, related analysis could be simpler. In this paper all the interactions are considered as repulsive. Accordingly, all the W-strengths are positive. Furthermore, it is safe to assume that all the three u /r ≥ 0 (because they do not contain nodes).

The set of W-strengths forms a matrix (i.e., the matrix-of-equations) with matrix elements . The determinant of is denoted by . The set of equation (4) has four forms of formal solutions, each would hold in a specific domain of r:

Form III: When all the three wave functions are nonzero in a domain, they must have the unique form aswhere is a determinant obtained by changing the l column of from (α 1, α 2, α 3) to (ε 1, ε 2, ε 3). is also a determinant obtained by changing the l column of to (1/2,1/2,1/2). Once all the parameters are given, the three Y are known because they depend only on α . However, the three X have not yet been known because they depend also on (ε 1, ε 2, ε 3). When Y is positive (negative), u /r goes down (up) with r. This point is notable because the main feature of the formal solution depends on the signs of {Y }.

The set {X } and the set {ε } are related as where , and is the algebraic cominor of α .

Form II: When one and only one of the wave functions is zero inside a domain (say, u /r = 0), the other two must have the unique form aswhere l, m, and n are in a cyclic permutation of 1-2-3 (the same in the follows),Once the parameters are given, the six (n′ ≠ n) are known, while the six have not yet. When is positive (negative), goes down (up) with r. When the Form II has u /r = 0, a more precise notation Form II is adopted for the detail.

Form I: When one and only one of the wave functions is nonzero in a domain (say, u /r ≠ 0), it must have the unique form asObviously, u /r in this form must descend with r. For the case u /r ≠ 0, the more precise notation Form I is adopted.

Form 0: In this form all the three wave functions are zero.

If a wave function (say, u /r) is nonzero in a domain but becomes zero when r = r , then a downward form-transition (say, from Form III to II) will occur at r . Whereas if u /r is zero in a domain but becomes nonzero when r = r , then a upward form-transition (say, from Form II to III) will occur at r . r 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 {ε } 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.

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 (or 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).

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 {ε } 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 and when r is given inside the domain (refer to eq. (10)). We define or ∞ (if or ≤ 0). Similarly, we define or ∞ (if or ≤ 0), and 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 I3. 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

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 II2, or in Form I3, 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.

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, and must both be negative. This violates the feature (i). If all the {u /r} emerge at the same place, all the three {Y } 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, II2-III-II1-I3 (it implies that the first domain has a Form II2, the next domain has a Form III, the third domain has a Form II1, while the last domain has a Form I3). The prescription on the linking appears as a number of requirements (inequalities) imposing on the W-strengths and the presumed {ε }. When all the and the {ε } 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 equations are further introduced, not only the scope but the values of the set {ε } 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 {ε }, 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.

State-transition and the singularity of the matrix

Based on the above approach, numerical calculations for two types of examples are performed. Related wave functions are plotted.

State-transition occurring at the singular point of the matrix-of-equations

Figure 1(b–e) is for II2-III-II1-I3, while (a) and (f) are for III-II1-I3. From (a) to (f) c is increasing while the other parameters remain unchanged. Thus these patterns demonstrate the effect of c .

Figure 1

Wave functions u 1/r (in solid line), u 2/r (in dash line), and u 3/r (in dash-dot line) against r. ħω and are used as units for energy and length. c is given at six values marked in the panels. Other parameters are fixed and are chosen quite arbitrary but having p and p both being negative. They are N  = 30000, N  = 11000, N  = 29000, c  = 4 × 10−4 (in ħωλ 3, the same for other strengths), c  = 1.4 × 10−3, c  = 1.2 × 10−3, c  = 3.8 × 10−4, c  = 4.2 × 10−4, γ  = γ  = γ  = 1. In 1b, the Forms of the solutions are marked in the associated domains separated by the vertical dotted lines.

For 1b as an example, the second domain is in Form III. For this form both {X } and {Y } are proportional to . Therefore, when , the wave functions will become extremely steep and the second domain will become extremely narrow as shown in (c). It turns out that, when , the matrix-of-equations becomes singular and accordingly . When c is close to and crosses over this critical value (from (c) to (d)), {Y } 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).

From the equality , it is straight forward to obtainwhereSimilarly, p and p can be defined by permuting the indexes.

Note that: the matrix will become singular and the transition will occur.

When p p  < 0, 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 does not assure the occurrence of the state-transition. Only if the other five strengths are so chosen that p p ≥ 0, the critical point could exist and the transition could occur.

deviates remarkably from the well known critical value for 2-BEC. Thus, the state-transition is remarkably affected by the influence of the third kind of atoms. However, if both c  → 0 and c  → 0 (i.e., the influence is removed), one can prove from eq. (13) that .

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.

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 increases and arrives at a critical value

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 n class="Chemical">TFA. In fact, in the earliest study of the 2-BEC, the instability in the neighborhood of the critical value (the singular point of the two-rank matrix) has been pointed out[8].

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-II3-I2 in (a) to (c), and is III-II3-I1 in (d) to (f), where the second domain has the Form II3 (namely, only the A- and B-atoms are contained in this domain). For the Form II3 the critical value of c is , which is the singular point of the sub-matrix of the equations for u 1/r and u 2/r only. When c 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, and 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 2-BEC. Incidentally, although the Form III is contained in Fig. 2, does not exist in this case due to p p  < 0.

Figure 2

The same as in Fig. 1 but the parameters are so given that p and p are in opposite signs. The details of parameters are N  = N  = N  = 30000, c  = 4 × 10−4 (in ħωλ 3, the same in the follows), c  = 2 × 10−3, c  = 1 × 10−3, c  = 11.5 × 10−4, c  = 10.5 × 10−4, γ  = γ  = γ  = 1.

Total energy of symmetric states and the great jump

When the total energy of the lowest symmetric state E is higher than the n class="Chemical">total energy of the lowest asymmetric state , the g.s. will be asymmetric. Thus, 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)where i = 1, 2 and 3. They are associated with A–, B–, and C– atoms, respectively. , , , and so on. Let N = N  + N  + N . Examples of E /N versus c are plotted via the solid lines shown in Fig. 3. The other parameters in (a) and (b) are the same as in Figs 1 and 2, respectively. A distinguished feature is the appearance of the great jump at (a) and (b). Note that, in Fig. 3(a), the crossing over does not cause an effect because the associated transition could occur only if the chain contains the building block II3, this building block is absent in Fig. 1. While in Fig. 3(b) does not exist because p p  < 0 as mentioned.

Figure 3

E /N (in solid line for the lowest symmetric state) and (in horizontal dash line for the lowest asymmetric state) versus c . Other parameters in (a,b) are the same as in Figs 1 and 2, respectively. The unit ħω is used for energy, and 10−4 ħωλ 3 for c . Note that the ranges of c in (a,b) are different.

Figure 3 confirms that the state-transition has caused a great change in E . 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 ) while P 2 will remarkably decrease. The decrease over takes the increase. We found that, for each B-atom, (E 2 + P 2)/N decreases from 1.80 to 1.41. On the other hand, for each A-atom, (E 1 + P 1)/N increases from 1.08 to 1.50. Since in this example, totally, E increases remarkably. This examples demonstrates that, although the critical value for the transition depends only on the strengths, the magnitude of the energy gap depends also on other parameters. The magnitude can be very large or quite small (say, in the above examples, the magnitude can be tuned by varying N and/or N ). Since all the A- and B- atoms are involved in the transition, the excitation is collective in nature.

Asymmetric states and their total energy

We know from the study of the 2-BEC[9, 11, 12] that, when V 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 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  < z 0 < R ). 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 and centered also at O. Then, we assumeif r ≤ R and z ≤ z 0. Otherwise, it is zero. Where is for the normalization.if r ≤ R and z > z 0. Otherwise, it is zero. Where andif r ≤ R . Otherwise, it is zero. Where When the values of R , R , 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 , can be obtained. The parameters R , R , and z 0 are considered as variable. Eventually, they fixed at the values that lead to the minimum of . The obtained via such a variational procedure is in general higher than its actual value. Thus, in any cases, if we found , 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 is remarkably lower than E when . Thus the jump provides a good opportunity for the lowest asymmetric state to replace the symmetric g.s. Whereas when , although is remarkably higher than E as shown in the figure and E will decrease further with the decrease of c , 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 and (i = 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 . When the n class="Chemical">TFA is removed, in a domain of r in which all the {u /r} are nonzero, the exact CGP can be written aswhere is the second-order derivative of u against r. The appearance of the common factor 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 . This is an important feature of the CGP. When a point in Σ crosses over , the factor 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 , the transition (denoted as trans-III) happens definitely.

Similarly, in a domain of r in which u /r = 0, the exact CGP can be written in a form in which both (u /r)2 and (u /r)2 are proportional to a common factor . Thus, for the type of entire solutions containing a Form II, the crossing over the surface will also lead to a great change in u /r and u /r, and accordingly another kind of state-transition (denoted as trans-IIn) occurs as shown in Fig. 2.

Let us define a subspace Σ as follows. When a set of parameters leads to an entire solution containing a Form III, then the associated point belongs to Σ, otherwise belongs to its complement. Let the part of the surface located inside Σ be denoted as σ . Then, σ appears as a boundary, the crossing over this boundary leads to the trans-III. Similarly, let denotes the subspace containing the points each leads to an entire solution containing a Form II3. Let the part of the surface located inside the subspace be denoted as . Then, the crossing over leads to the transition trans-II3. We can further define and in a similar way. These four surfaces (σ and the three ) 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 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

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 /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.

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 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 inten class="Chemical">rband coupling[31-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).

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.

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.