GSI method :

(i) Given A ∈ R^{N ×N}, α1 ≥ 0, α2 ≥ 0, ρ1, ρ2, β > 0 and u⁽⁰⁾_j > 0 with ku^jk2 = 1, j = 1, 2, let n = 0;

(ii) Repeat n: until convergence,

Use the Shift-Invert Arnoldi algorithm [16, 18] or the Jacobi-Davidson [17] to solve the minimal eigenvalue λ⁽ⁿ⁺¹⁾₁ with eigenvector u⁽ⁿ⁺¹⁾₁ and λ⁽ⁿ⁺¹⁾₂ with eigenvector u⁽ⁿ⁺¹⁾₂ of A⁽ⁿ⁺¹⁾₁ and A⁽ⁿ⁺¹⁾₂ , respectively, where

A⁽ⁿ⁺¹⁾₁ := A + α1[[u⁽ⁿ⁾₁ ² ]] + βρ1[[u⁽ⁿ⁾₂ ²]], (6.1) A⁽ⁿ⁺¹⁾₂ := A + α2[[u⁽ⁿ⁾₂ ² ]] + βρ2[[u⁽ⁿ⁺¹⁾₁ ²]]. (6.2) (iii) Compute the residual,

res⁽ⁿ⁺¹⁾_j = A⁽ⁿ⁺¹⁾_j u⁽ⁿ⁺¹⁾_j − λ⁽ⁿ⁺¹⁾j u⁽ⁿ⁺¹⁾_j , j = 1, 2. (6.3)

(iv) Ifkres⁽ⁿ⁺¹⁾j k² < Tol, j = 1, 2, then stop; else n← n + 1, go to Repeat.

An essential equivalent condition of convergence of GSI method has been proven in [8].

Theorem 6.1. [8] The GSI method converges to (u^∗₁, λ^∗₁, u^∗₂, λ^∗₂) locally starting with an initial in T^o ×T^o if and only if u^∗ = (u^∗₁, u^∗₂) is a strictly local minimum of E(u) as in (1.11b).

Theorem 4.1 in Section 4 shows that a supercritical pitchfork bifurcation of E(u) (1.11b) occurs at some finite value of β = β₁^∗ (or ˜β₁^∗). From Theorem 6.1 we see that the GSI method can be used to compute two strictly local minima of E(u) for β > β₁^∗ (or ˜β₁^∗). In [6] a continuation BSOR-Lanczos-Galerkin (BSOR-LG) method has been developed for computing other positive bound states of a multi-component BEC. The solution curve of NAEP (1.10) with respect to β can then be followed by the proposed

method [6] for β ≥ β2^∗ (or ˜β₁^∗), where β₂^∗ (or ˜β₁^∗) is the second singular bifurcation point of C^α2,ρ2.

We compute the solution curve of ground/positive bound states of NAEP (1.10) on a rectangle domain [−1, 1] × [−2/3, 2/3] with grid size h = 1/30, with (α¹, α2, ρ1, ρ2) = (0.1, 0.1, 1, 1) and (α1, α2, ρ1, ρ2) = (0.1, 0.12, 1, 1.1), respectively, by using the GSI method and the BSOR-LG continuation method. In Figure 6.1 and Figure 6.2, we plot the bifurcation diagrams for two cases of the solution curves versus β ∈ (0, 50). The nodal domains of positive bound states are attached near the solution curve, where the left figure is the level set of u^∗₁ and the right figure is the level set of u^∗₂.

For (α1, α2, ρ1, ρ2) = (0.1, 0.1, 1, 1), in Figure 6.1 we see that the NAEP (1.10) has only identical ground state solutions for 0≤ β ≤ β1^∗, and the solution curve undergoes a supercritical pitchfork bifurcation at β = β₁^∗. The ground states begin to separate at β > β₁^∗. From Theorem 4.1 and Theorem 5.2 it is expected that the ground states of (1.10) should be little by little mutually separated with continually increasing β. The structure of the phase separation will reach a stage of segregated nodal domains, when β ≈ 10⁵. Next, we come back to the bifurcation point β₁^∗ on the primal stalk and trace the solution curveC^α1,ρ1 with identical bound states for β > β₁^∗. By path following, we encounter a sequence of bifurcation points {β1^∗, β₂^∗, β₃^∗, β₄^∗, β₅^∗} on the primal stalk. For each bifurcation branch at β_k^∗, k = 2, . . . , 5, if we trace the solution curve C^α1,ρ1 with β > β_k^∗ a new structure of positive bound states is found.

In Figure 6.2 we plot the solution curve of NAEP (1.10) with (α1, α2, ρ1, ρ2) = (0.1, 0.12, 1, 1.1). In Table 4.2 we see that the vectors ξ_i,1 in Lemma 4.1, i = 1, 2, 4, 5 at the 1-st, the 2-nd, the 4-th and the 5-th bifurcation points are anti-symmetric and the vector ξ_3,1 at the 3-rd bifurcation point is symmetric when (α1, α2, ρ1, ρ2) = (0.1, 0.1, 1, 1). Therefore, from Remark 4.2 (iii) and Theorem 4.1 one can shown that the negative gradient flow of the FOP (1.11a) undergoes supercritical pitchfork

bifur-cation at ˜β₁^∗, ˜β₂^∗, ˜β₄^∗, ˜β₅^∗ with (α1, α2, ρ1, ρ2) = (0.1, 0.12, 1, 1.1). The most interesting phenomenon here is that the solution curve of the NAEP (1.10) for the nonidentical case breaks at the 3-rd bifurcation point, i.e., the 3-rd supercritical pitchfork bifurca-tion at β = β₃^∗with α1 = α2, ρ1 = ρ2disappears and becomes a saddle-node bifurcation at β = ˜β₃^∗ when the intra- and inter-component scattering lengths are chosen slightly different with (α1, α2, ρ1, ρ2) = (0.1, 0.12, 1, 1.1) (See Figure 6.2).

7 Conclusions

In this paper, we mainly prove that the negative gradient flow of the FOP (1.11) under-goes supercritical pitchfork bifurcations for ground/positive bound state solutions at some finite values of repulsive scattering lengths. Furthermore, symmetric pitchfork bi-furcations for ground/positive bound states occur when two intra- and inter-component scattering lengths are equal. We also show that ground/positive bound states may re-pel each other and form segregated nodal domain as the repulsive scattering length goes to infinity.

In the future, we are interested in proving the existence of a multi-pitchfork bifur-cation of a multi(m)-component BEC with m ≥ 3, at some finite value of β^∗, and the symmetry of ground state solutions with respect to Ω, at the finite value β^∗ when a m-component BEC has equal intra- and inter-component scattering lengths.

A Appendix

Proof of Theorem 4.2. As in the proof of Theorem 4.1, hereafter, we set βnew = βold− β˜₁^∗. It suffices to show that a(0) = D⁴f (0, 0)q1(0)⁴ and b(0) = D³f (0, 0)q1(0)². To

this end, from (4.45) and (4.49a) we first note that Differentiating θ twice and three times, respectively, with respect to v₁ and noting that q^>_1,i(0)q1,i(0) = pi for some p1, p2 > 0 satisfying p1+ p2 = 1, we have

Multiplying (A.4) by q1(0) from the right, using (A.2) and using the fact ξ^>_1,1u^∗₁ = 0, we have

D²T(v^∗(0))q1(0)²

= 0^>_{N −1},−ξ1,1,N(0)²

u^∗_1,N(0) − p1

u^∗_1,N(0), 0^>_{N −1},−ξ1,2,N(0)²

u^∗_2,N(0) − p2

u^∗_2,N(0)

!>

= ν2(0). (A.5)

Analogously, using (A.3b), we also have

D³T(v^∗(0))q₁(0)³ = 0^>_{N −1},3ξ1,1,N(0)³

u^∗_1,N(0)² + 3p1ξ1,1,N(0)

u^∗_1,N(0)² , 0^>_{N −1},3ξ1,2,N(0)³

u^∗_2,N(0)² +3p2ξ1,2,N(0) u^∗_2,N(0)²

!_>

.

Now, differentiating f (z, β) and gj four times with respect to z at z = 0 and β = 0, and from (4.43) we get

D⁴f (0, 0)q1(0)⁴ =D⁴_uE(u^∗(0))ξ₁(0)⁴ + 6D³_uE(u^∗(0))ξ₁(0)²(D²T(v^∗(0))q1(0)²) + 3D_u²E(u^∗(0))(D²T(v^∗(0))q1(0)²)²

+ 4D_u²E(u^∗(0))ξ₁(0)(D³T(v^∗(0))q1(0)³)

+ DE(u^∗(0))(D⁴T(v^∗(0))q1(0)⁴) (A.6) and

D⁴_ugj(u^∗(0))ξ₁(0)⁴+ 6D_u³gj(u^∗(0))ξ₁(0)²(D²T(v^∗(0))q₁(0)²)

+ 3D_u²gj(u^∗(0))(D²T(v^∗(0))q₁(0)²)²+ 4D²_ugj(u^∗(0))ξ₁(0)(D³T(v^∗(0))q₁(0)³) + Dugj(u^∗(0))(D⁴T(v^∗(0))q₁(0)⁴) = 0. (A.7)

Subtracting (A.6) from ˆλ1× (A.7)j=1+ ˆλ2× (A.7)j=2 and using (4.37) we have D⁴f (0, 0)q1(0)⁴

= D⁴_uE(u^∗(0))ξ₁(0)⁴+ 6D_u³E(u^∗(0))(D²T(v^∗(0))q1(0)²)ξ₁(0)² + 3h

D²_uE(u^∗(0))−

ˆλ1D²_ug1(u^∗(0)) + ˆλ2D_u²g2(u^∗(0))i

(D²T(v^∗(0))q1(0)²)² + 4h

D²_uE(u^∗(0))−

ˆλ1D²_ug1(u^∗(0)) + ˆλ2D_u²g2(u^∗(0))i

ξ₁(0)(D³T(v^∗(0))q1(0)³).

(A.8) With a = (a^>₁, a^>₂)^> we compute

D⁴_uE(u)a⁴ =ρ2(6α1a^>₁[[a1◦ a1]]a1+ 2ρ1β˜₁^∗a^>₁[[a2◦ a2]]a1) + ρ1(6α2a^>₂[[a2 ◦ a2]]a2+ 2ρ2β˜₁^∗a^>₂[[a1◦ a1]]a2)

+ 4ρ1ρ2β(a^>₁[[a2◦ a1]]a2+ a^>₁[[a1◦ a2]]a2). (A.9) Then, the first term of (A.8) is positive by evaluating (A.9) at β = ˜β₁^∗, (a^>₁, a^>₂)^> = (ξ^>_1,1(0), ξ^>_1,2(0))^> and u^∗(0) = (u^∗₁(0), u^∗₂(0)), i.e.,

D⁴_uE(u^∗(0))ξ₁(0)⁴ =6ρ2α1ξ^>_1,1(0)[[ξ_1,1(0) ² ]]ξ_1,1(0) + 6ρ1α2ξ^>_1,2(0)[[ξ_1,2(0) ² ]]ξ_1,2(0) + 12ρ1ρ2β˜₁^∗ξ^>_1,1(0)[[ξ_1,2(0) ²]]ξ_1,1(0) > 0. (A.10) Substituting (A.10), (A.1), (A.5) , (4.39b) into (A.8), we obtain

D⁴f (0, 0)q1(0)⁴ =6ρ2α1ξ^>_1,1(0)[[ξ_1,1(0) ²]]ξ_1,1(0) + 6ρ1α2ξ^>_1,2(0)[[ξ_1,2(0) ²]]ξ_1,2(0)

+12ρ1ρ2β˜₁^∗ξ^>_1,1(0)[[ξ_1,2(0) ² ]]ξ_1,1(0) + 6ν1(0)^>ν2(0) + 3ν2(0)^>J(0)ν2(0)

≡a(0). (A.11)

Similar to part (iii) in the proof of Theorem 4.1, one can also compute D³f (0, 0)q₁(0)²p = D_u³E(u^∗(0))ξ₁(0)²(DT(v^∗(0))p)

+ 2(D_u²E(u^∗(0))− (ˆλ1D_u²g₁(u^∗(0)) + ˆλ₂D_u²g₂(u^∗(0))))ξ₁(0)(D²T(v^∗(0))q₁(0)p) +

D_u²E(u^∗(0))− (ˆλ1D_u²g₁(u^∗(0)) + ˆλ₂D²_ug₂(u^∗(0)))

(DT(v^∗(0))p)(D²T(v^∗(0))q₁(0)²) (A.12)

for any p∈ R^{N −1}. Substituting (A.1), (4.39b), (A.5) into (A.12), we have

D³f (0, 0)q1(0)² = DT(v^∗(0))^>ν1(0) + DT(v^∗(0))^>J(0)ν2(0) = b(0). (A.13) From (A.11) and (A.13), we complete the proof.

References

[1] P. Ao and S. T. Chui. Binary Bose-Einstein condensate mixtures in weakly and strongly segregated phases. Phys. Rev. A, 58:4836–4840, 1998.

[2] W. Z. Bao. Ground states and dynamics of multi-component Bose-Einstein con-densates. SIAM Multiscale Modeling and Simulation, to appear.

[3] W. Z. Bao and Q. Du. Computing the ground state solution of Bose-Einstein condensates by a normalized gradient flow. SIAM J. Sci. Comput., to appear.

[4] M. S. Bazaraa, H. D. Sherali, and C. M. Shetty. Nonlinear programming: theory and algorithms. Wiley, New York, 2 edition, 1993.

[5] A. Berman and R. J. Plemmons. Nonnegative matrices in the mathematical sci-ences. Academic Press, New York, 1979.

[6] S. M. Chang, Y. C. Kuo, W. W. Lin, and S. F. Shieh. A continuation BSOR-Lanczos-Galerkin method for positive bound states of a multi-component Bose-Einstein condensate. NCTS preprints in Math#2002-24, Tsing-Hua Unv., Tai-wain, 2003.

[7] S. M. Chang, C. S. Lin, T. C. Lin, and W. W. Lin. Segregated nodal domains of two-dimensional multispecies Bose-Einstein Condensates. Physica D, to appear.

[8] S. M. Chang, W. W. Lin, and S. F. Shieh. Gauss-Seidel-type methods for en-ergy states of a multi-component Bose-Einstein condensate. NCTS preprints in Math#2002-25, Tsing-Hua Unv., Taiwain, /J. Comp. Phy., to appear.

[9] B. D. Esry and C. H. Greene. Spontaneous spatial symmetry breaking in two-component Bose-Einstein condensates. Phys. Rev. A, 59:1457–1460, 1999.

[10] B. D. Esry, C. H. Greene, J. P. Burke Jr, and J. L. Bohn. Hartree-Fock theory for double condensates. Phys. Rev. Lett., 78:3594–3597, 1997.

[11] S. Gutpa, Z. Hadzibabic, M. W. Zwierlein, C. A. Stan, K. Dieckmann, C. H.

Schunck, E. G. M. van Kempen, B. J. Verhaar, and W. Ketterle. Radio-frequency spectroscopy of Ultracold Fermions. Science, 300:1723–1726, 2003.

[12] D. S. Hall, M. R. Matthews, J. R. Ensher, C. E. Wieman, and E. A. Cornell. Dy-namics of component separation in a binary mixture of Bose-Einstein condensates.

Phys. Rev. Lett., 81:1539–1542, 1998.

[13] E. H. Lieb, R. Seiringer, and J. Yngvason. Bosons in a trap: a rigorous derivation of Gross-Pitaevskii energy functional. Phys. Rev. A, 61:043602, 2000.

[14] C. J. Myatt, E. A. Burt, R. W. Ghrist, E. A. Cornell, and C. E. Wieman. Produc-tion of two overlapping Bose-Einstein condensates by sympathetic cooling. Phys.

Rev. Lett., 78:586–589, 1997.

[15] C. Robinson. Dynamical systems: stability, symbolic dynamics, and chaos. CRC Press, 1995.

[16] A. Ruhe. Rational Krylov: A practical algorithm for large sparse nonsymmetric matrix pencils. SIAM J. Sci. Comput., 19(5):1535–1551, 1998.

[17] G. L. G. Sleijpen and H. A. van der Vorst. A Jacobi-Davidson iteration method for linear eigenvalue problems. SIAM J. Matrix Anal. Appl., 17:401–425, 1996.

[18] D. C. Sorensen. Implicit application of polynomial filters in a k-step Arnoldi method. SIAM J. Matrix Anal. Appl., 13:357–385, 1992.

[19] E. Timmermans. Phase separation of Bose-Einstein condensates. Phys. Rev. Lett., 81:5718–5721, 1998.

[20] M. Trippenbach, K. G´oral, K. Rz¸a˙zewski, B. Malomed, and Y. B. Band. Structure of binary Bose-Einstein condensates. J. Phys. B: At. Opt. Phys., 33:4017–4031, 2000.

[21] S. Wiggins. Introduction to applied nonlinear dynamical systems and chaos.

Springer-Verlag, 1990.

β₁^∗ β₂^∗ β₃^∗ β₄^∗ β β₅^∗

Figure 6.1: Bifurcation diagram of a two-component BEC on the rectangle domain ([−1, 1] × [−2/3, 2/3]) with grid size h = 1/30 and (α1, α2, ρ1, ρ2) = (0.1, 0.1, 1, 1).

The bifurcation points are, respectively, β₁^∗ = 3.1940, β₂^∗ = 6.1873, β₃^∗ = 7.8002, β₄^∗ = 10.956, β₅^∗ = 16.278. The level sets of u1 (left) and u2 (right) are attached near the solution curve.

β˜₁^∗ β˜₂^∗ β˜₃^∗ β˜₄^∗ β

β₅^∗ β ≈ 50

Figure 6.2: Bifurcation diagram of a two-component BEC on the rectangle domain ([−1, 1] × [−2/3, 2/3]) with grid size h = 1/30 and (α1, α2, ρ1, ρ2) = (0.1, 0.12, 1, 1.2).

The bifurcation points are, respectively, ˜β₁^∗ = 3.0525, ˜β₂^∗ = 5.9281, ˜β₃^∗ = 8.188, ˜β₄^∗ = 10.47, ˜β₅^∗ = 15.78. The level sets of u1 (left) and u2 (right) are attached near the solution curve.