Exploring Bistability in Rotating Bose-Einstein
Condensates by a Quotient Transformation Invariant
Continuation Method
Yueh-Cheng Kuo∗ Wen-Wei Lin† Shih-Feng Shieh‡ Weichung Wang§
June 4, 2010
Abstract
We study the discrete nonlinear Schr¨odinger (DNLS) equations that model rotating Bose-Einstein Condensates (BEC) both analytically and numerically. Due to the difficulties associated with transformation invariant solutions, standard continuation methods may not properly follow the solution curves of the DNLS equations. We propose a quotient transformation invariant con-tinuation method to circumvent this obstacle. We also analyze the bifurca-tion properties of the primal stalk solubifurca-tion curve corresponding to the DNLS equations for an isotropic trap. In numerical computation, the existence of a bistable region corresponding to the bound states with a 0- or 1-vortex is shown. This finding not only agrees with the physics of the experimental phenomena, but also explains why a 0- or 1-vortex may be observed within a certain region that has an angular velocity. Numerical evidence shows that trap potentials have little effect on the width of the bistable regions. In contrast, intra-component scattering length significantly affects the bistable region.
Keywords. Rotating Bose-Einstein condensates, nonlinear Schr¨odinger equa-tions, quotient transformation invariant continuation method, bistability. ∗
Department of Applied Mathematics, National University of Kaohsiung, Kaohsiung 811, Tai-wan ([email protected]).
†
Department of Applied Mathematics, National Chiao Tung University, Hsinchu 300, Taiwan ([email protected]).
‡
Department of Mathematics, National Taiwan Normal University, Taipei 116, Taiwan ([email protected]).
§
Department of Mathematics, National Taiwan University, Taipei 106, Taiwan ([email protected]).
1
Introduction
Bose-Einstein condensates (BEC) have been produced and extensively studied in laboratories in the last 15 years. Studies include the exploration of quantized vor-tex states and their connections with superfluidity [1, 20]. In many circumstances, a laser beam rotating with an angular velocity ω is applied to the magnetic trap holding the atoms to create a harmonic and anisotropic potential. Vortex nucle-ation [20] and vortex arrays [21] have also been observed.
In the rotating frame, a BEC with an external driven field at temperatures T much smaller than the critical temperature Tc [23] can be well described by
time-independent nonlinear Schr¨odinger (NLS) equations [2, 4, 13]. In particular,
−1 2∇ 2φ(x) + V (x)φ(x) + α|φ|2φ(x) + ωι∂ θφ(x) = λφ(x), (1.1a) for x ∈ Ω ⊆ R2 with Z Ω |φ(x)|2dx = 1, (1.1b)
where Ω is a bounded, smooth domain and ι = √−1. In (1.1a), V (x) ≥ 0 is the magnetic trapping potential, ∂θ = x∂y − y∂x is the z-component of the
an-gular momentum, ω is the anan-gular velocity of the rotating laser beam, α denotes the intra-component scattering length and λ is the chemical potential or eigen-value. The ground state solution of a BEC can be found by minimizing the energy functional E(φ) under the conditions given by (1.1b) [2, 4]. That is,
Minimize φ E(φ) subject to Z Ω |φ(x)|2dx = 1 and φ(x) = 0 for x ∈ ∂Ω, (1.2a) where E(φ) = Z Ω 1 2|∇φ| 2+ V |φ|2+α 2|φ| 4+ ωιφ∗∂ θφ (1.2b)
and φ∗ denotes the complex conjugate of φ. In (1.2b), the value of R
Ωφ ∗∂
θφ is
purely imaginary. The NLS equations (1.1) are also regarded as the Euler-Lagrange equation of the optimization problem given by (1.2).
To numerically solve a rotating BEC, a continuous normalized gradient flow (CNGF) method with a backward Euler finite difference discretization has been proposed in [6] for computing the ground, symmetric and central vortex states. The paper analyzes the existence/nonexistence problem for ground states and
finds that symmetric and central vortex states are independent of the angular rotational momentum. Furthermore, a continuation method has been proposed in [7] for computing the energy levels of a rotating BEC. Also, non-rotating BECs have been studied numerically and analytically in [4, 5, 8, 9, 16].
While several numerical studies on rotating BECs and the corresponding dis-crete nonlinear Schr¨odinger equations have been conducted, computational chal-lenges still persist. We focus on the following computational issue. Due to the transformation invariant property (in particular, rotational invariance) of the so-lution set of the NLS equations (1.1) or its corresponding discrete version, the solution set forms a two-dimensional manifold. Therefore, each solution point un-dergoes the singularity in a standard continuation scheme, making it difficult to trace the solution set. In addition, the bistability, i.e. the saddle points of (1.2), further complicates solutions based on gradient methods.
Our approach for overcoming these numerical difficulties is to quotient the solution set to obtain a one-dimensional solution curve. Our findings are briefly summarized as follows.
• We propose the quotient transformation invariant continuation method that follows the one-dimensional solution curves. These curves are obtained by quotienting the two-dimensional solution’s manifold. We then elaborate on the details of the method, including the prediction-correction path procedure and the detection of bifurcation points. We also propose an easy way to check whether the computed solutions are stable.
• We analyze the bifurcation properties of the primal stalk solution curve of the model equations analytically.
• We assert the existence of bistable solutions to the equations that model a rotating BEC via numerical methods. This bistability is imperceptible in time-dependent numerical methods and experiments. The bistability agrees with the experimental phenomena that either a 0- or 1-vortex is observed in a certain region that has an angular velocity [20]. Finally, our numerical exper-iments suggest that a potential trap induces linear changes to the bistable regions, and the intra-component scattering length affects the bistable re-gions exponentially.
This paper is organized as follows. In Section 2 we develop the quotient trans-formation invariant continuation method for solving the NLS equations numeri-cally. In Section 3 we show bifurcation for a rotating BEC. In Section 4 we present the numerical results for stable and bound state solutions of a rotating BEC by using the QTICM. Finally, a conclusion is given in Section 5.
1.1 Discretization of the NLS equations
To investigate the NLS equations (1.1) and the associated optimization problem, we can discretize the differential equations, for example, by a finite difference method. In particular, we suppose that the operator −12∇2+V (x) and the angular
momentum ∂θin (1.1a) are discretized by the central difference approximation with
a grid size h. Let A ∈ RN ×N be the discretization matrix of the operator −12∇2+
V (x) combined with the Dirichlet boundary condition. We can see that A is an irreducible and symmetric positive definite matrix with non-positive off-diagonal entries (i.e. an irreducible M-matrix). Let S ∈ RN ×N be the discretization matrix of ∂θ where S is a skew-symmetric matrix.
We denote the approximation of the wave function φ in (1.1a) as uc= u1+ιu2,
where u1, u2 ∈ RN. Then, the discrete nonlinear Schr¨odinger (DNLS) equations
corresponding to (1.1) can be formulated as
Auc+ α¯uc◦ u c2 + ωιSuc= λuc, (1.3a)
uHc uc=
1
h2. (1.3b)
Here u c2 = uc◦ ucand ◦ denotes the Hadamard product. Similarly, the
optimiza-tion problem (1.2a) becomes the finite-dimensional optimizaoptimiza-tion problem: Minimize
u1,u2∈RN
E(u1, u2)
subject to u>1u1+ u>2u2 = h12,
(1.4a)
with the energy functional
E(u1, u2) = h2 uHc Auc+ α 2u 2H c u c2 + ωιuHc Suc = h2 " 2 X i=1 u>i Aui+ α 2(u 2 1 + u 2 2 ) >(u 2 1 + u 2 2 ) + 2ωu > 2Su1 # . (1.4b)
Throughout this paper, we use the bold face letters or symbols to denote a matrix or a vector. For u = (u1, . . . , un)>, v = (v1, . . . , vn)> ∈ Cn, u ◦ v =
(u1v1, . . . , unvn)> denotes the Hadamard product of u and v, u r = u ◦ · · · ◦
u denotes the r-time Hadamard product of u, |u| = (|u1|, . . . , |un|)> and kuk2
denote the vector for the absolute value of each component and the 2-norm of u, respectively, and [[u]] := diag(u) denotes the diagonal matrix of u. For A ∈ Rm×n, N (A) and R(A) denote the null space and column space of A, respectively. For A ∈ Rn×n, A > 0 (≥ 0) denotes a symmetric positive (semi-positive) definite matrix.
2
Quotient Transformation Invariant Continuation Method
In this section, we develop a quotient transformation invariant continuation method (QTICM) to solve the DNLS equations (1.3). First, we derive the parameter de-pendent polynomial systems that are solved at each iteration of the continuation method in Section 2.1. A strategy for circumventing the difficulties associated with transformation invariant solutions is presented in Section 2.2; this strategy is discussed in detail in Section 2.3. Methods for detecting the bifurcation points and stability of the solutions are proposed in Sections 2.4 and 2.5, respectively.
2.1 A parameter dependent polynomial system
A continuation method usually follows solution curves by solving sequences of parametrized polynomial systems. The polynomial systems considered here are described as follows.
Letting u1 and u2 be the real and imaginary parts of uc, respectively, we can
rewrite uc= u1+ ιu2 ∈ CN and the DNLS equations (1.3) as:
Au1+ α([[u 12]] + [[u2 2]])u1− ωSu2= λu1, (2.1a)
Au2+ α([[u 12]] + [[u2 2]])u2+ ωSu1= λu2, (2.1b)
u>1u1+ u>2u2 =
1
h2. (2.1c)
Now, we denote p as the positive continuation parameter, such that p is incorpo-rated into the system (2.3) by letting
α = α0+ µ0p,
ω = ω0+ ν0p,
(2.2)
for the known constants α0, µ0, ω0and ν0. Furthermore, by letting u = (u>1, u>2)>∈
R2N and z = (u>, λ)> ∈ R2N +1, we can rewrite the DNLS equations (2.1) as a parameter dependent polynomial system:
G(z, p) = 0, (2.3) where G ≡ (G1, G2, g) : R2N +1× R → R2N +1 is given by G1(z, p) = Au1+ α([[u 2 1 ]] + [[u 2
2 ]])u1− ωSu2− λu1, (2.4a)
G2(z, p) = Au2+ α([[u 12]] + [[u2 2]])u2+ ωSu1− λu2, (2.4b)
g(z) = 1 2 u>1u1+ u>2u2− 1 h2 . (2.4c)
Consequently, the Jacobian of G is given by DG = [Gz, Gp] ∈ R(2N +1)×(2N +2). (2.5) Thus, Gz(z, p) = [Gu|Gλ] = " H(z, p) −u −u> 0 # ∈ R(2N +1)×(2N +1), (2.6a) Gp(z, p) = µ0 −Su2 Su1 0 + ν0
([[u 12 ]] + [[u 22]])u1
([[u 12 ]] + [[u 22]])u2
0 ∈ R 2N +1, (2.6b) H(z, p) = " b A + 2α[[u 12]] 2α[[u1◦ u2]] − ωS 2α[[u1◦ u2]] + ωS A + 2α[[ub 2 2 ]] # (2.7)
is a 2N × 2N symmetric matrix, and
b
A = A + α([[u 12 ]] + [[u 22 ]]) − λIN. (2.8)
2.2 The transformation invariant solutions and quotient solutions
Now, we discuss why a standard continuation method is inappropriate for our problem and how we may circumvent the difficulty. Let the solution set M of (2.3) be defined as follows:
M = {(z>, p)>| G(z, p) = 0}. (2.9)
We define a function u(ϑ) : [0, 2π] → R2N by
u(ϑ) = "
cos ϑu1+ sin ϑu2
− sin ϑu1+ cos ϑu2
#
. (2.10)
It is clear that if (u>1, u>2, λ, p)>∈ M and both λ and p are fixed, then
G(u(ϑ), λ, p) = 0 (2.11)
for all ϑ ∈ [0, 2π]. That is, the solution set M is actually a two-dimensional man-ifold on R2N +2 containing the transformation invariant solutions u(ϑ)’s. In other words, for a particular solution (u(0)>, λ, p)>∈ M, all of the solutions in the form (u(ϑ)>, λ, p)> have the same energy and density functions with the same shape. Namely, E(u(0)) = E(u(ϑ)) and |u(0)| 2 = |u(ϑ)| 2 for all ϑ ∈ [0, 2π]. This
two-dimensional solution manifold thus may prevent a standard continuation method from successfully following a desired solution curve by varying the parameter p.
(a) '( )ϑ u ( )ϑ u prediction direction (b)
Figure 2.1: Conceptual solution pipelines that show: (a) why a standard continua-tion method may not lead to meaningful numerical solucontinua-tions, and (b) our approach for computing the quotient transformation invariant solutions.
This observation can be understood conceptually with Figure 2.1. We image the shape of the solution is a pipeline in R2N +2; the solution curve (u(ϑ)>, λ, p)> for a certain λ and p is given by the black line in the middle of the pipeline of Figure 2.1-(b). A standard continuation method may irregularly circle around the pipeline solution set as shown in Figure 2.1-(a).
Now, let
G(y(s)) = 0 (2.12)
be the arc length parameter, where s ∈ R. By differentiating equation (2.12) with respect to s, we obtain:
DG(y(s)) ˙y(s) = 0, (2.13)
where DG(y(s)) is defined in (2.5) and ˙y(s) = ( ˙z(s)>, ˙p(s))> denotes a tangent vector to C at y(s). Equation (2.13) suggests that the tangent vector ˙y(s) to C at y(s) is the nontrivial solution of the (2N + 1) × (2N + 2) homogeneous system DG(y(s))w = 0, when DG(y(s)) has a full row rank. However, the Jacobian matrix DG(y(s)) considered here does not have a full row rank, as the dimension of the solution set M is 2. In this case, standard continuous methods may not follow the solution set successfully.
To avoid circling around the solution manifold, our approach is to consider a quotient solution set by applying additional hyperplane constraints on the system. First, we parametrize the solution curve (2.3) by introducing arc length parameter s, so that the equivalent solutions u(ϑ) with ϑ ∈ [0, 2π] are quotient. Furthermore, the quotient solution curve can be parametrized as:
Now, we shall discuss how we can trace the quotient solution curve C by using a continuation method with additional hyperplane constraints. To do so, we can compute the tangent vector of u(ϑ) for a certain ϑ. Without losing generality, we may choose ϑ = 0 to obtain
∂u ∂ϑ(0) = (u > 2, −u > 1) > , (2.15)
and then develop a quotient transformation invariant continuation method by incorporating the following prediction and correction schemes.
1. In the prediction step, a prediction direction ˙yi = ( ˙z>i , ˙pi)> should
sat-isfy (2.13) and the vector ˙zi should be orthogonal to the tangent vector ∂u
∂ϑ(0) >, 0>
.
2. In the correction step, we compute the correction vector that satisfies an additional hyperplane constraint whose normal vector is ∂u∂ϑ(0)>, 0, 0>
. By doing so, the correction solution is located in the quotient curve C (2.14). Figure 2.1-(b) demonstrates this idea visually.
We have introduced the main ideas of the quotient transformation invariant continuation method. In the next subsection, we will discuss in detail how this QTICM follows the parametrized quotient solution curve with modified prediction and correction steps.
2.3 Prediction and correction
Let yi= (z>i , pi)>∈ R2N +2 be a point that has been accepted as an approximate
solution on C. A “good” prediction direction ˙yi = ( ˙z>i , ˙pi)> should satisfy (2.13)
and the vector ˙z should be orthogonal to ∂u∂ϑ(0)>, 0>
. Consequently, we apply an additional hyperplane constraint to equation (2.13):
a>ϑ˙z = 0, (2.16) where aϑ= ∂u∂ϑ(0)>, 0
>
= (u>2, −u>1, 0)>. It follows that the prediction direction ˙
y ∈ R2N +2 should satisfy the bordered linear system Gz Gp a>ϑ 0 c>i ci y˙i= 0 0 1 , (2.17)
where (c>i , ci)> ∈ R2N +2 is a suitable constant vector. Note that equation (2.17)
can be geometrically interpreted as follows. We first use the Euler predictor yi+1,1 = yi+ hi( ˙yi/k ˙yik2)
to predict a new point yi+1,1, where hi > 0 is a suitable step length and ˙yi is the
tangent vector at yi that is obtained by solving the bordered linear system (2.17).
In the correction step, the next acceptable solution yi+1 must pass through
the hyperplane whose normal vector is (a>ϑ, 0)> with the Euler predictor yi+1,1=
(zi+1,1, pi+1,1). In other words, the solution curve C is now determined by the
underlying system equations: G(z, p) = 0, a>ϑz = a>ϑzi+1,1, ˙z>i z + ˙pip = ˙y>i yi+1,1.
Starting from the predictor yi+1,1, Newton’s method is typically used to compute
the correctors. By setting yi+1,l+1 = yi+1,l + δl for l = 1, 2, . . ., we solve the
bordered linear system using an iterative method: Gz(yi+1,l) Gp(yi+1,l) a>ϑ 0 ˙z>i p˙i δl= −G(yi+1,l) −ρϑ,l −ρl , (2.18)
with ρϑ,l = a>ϑ(zi+1,l− zi+1,1) and ρl= ˙y>i (yi+1,l− yi+1,1). If the sequence {yi+1,l}
converges numerically for l = l∞, we accept yi+1= yi+1,l∞as a new approximation
to the solution curve C.
While both linear systems (2.17) and (2.18) are overdetermined and have unique solutions, an efficient way to solve these two systems is to rewrite them in the form: " B f g> γ # " x σ # = " q ρ # , x(2N + 2) = 0, (2.19) where B = " Gz aϑ a>ϑ 0 # ∈ R(2N +2)×(2N +2) (2.20)
is symmetric, f , g, q ∈ R2N +2 and x(2N + 2) denotes the (2N + 2)-th entry of x. The linear system (2.19) can be easily solved by the well-known block elimination algorithm (see [10, 12, 15]).
We conclude this discussion of prediction and correction with the following two remarks.
(i) The main step in the prediction or correction scheme is to solve linear systems in the form of Bξ = g.
(ii) Since the linear systems (2.17) and (2.18) are consistent, equation (2.19) is also consistent. Thus, the solution x of first equation of (2.19) satisfies x(2N + 2) = 0.
2.4 Detecting bifurcation points
We have discussed how we may follow the solution curve by the QTICM. Now we focus on detecting the bifurcation points, or the singular points, along the solution curves.
Let C be the solution curve defined in (2.14), where y(s) ∈ C and
J(s) = " Gz(y(s)) Gp(y(s)) a>ϑ 0 # ∈ R(2N +2)×(2N +2), (2.21a) J(s) = " Gz(y(s)) a>ϑ # ∈ R(2N +2)×(2N +1). (2.21b)
As described in [3, 11, 15], a point y(s) ∈ C is said to be a regular point if rank(J(s)) = 2N + 1 (i.e., dim N (J(s)) = 1) and a singular point if rank(J(s)) ≤ 2N (i.e., dim N (J(s)) ≥ 2). For a regular point y(s), the tangent vector ˙y(s) is uniquely determined by the linear system (2.17).
We propose to use the inverse power method to detect the singularity of the symmetric matrix B(s) described in (2.20) at the point y(s) ∈ C. Furthermore, the tangent vectors at the singular point y(s∗) ∈ C can also be found using the methods reported in [15, pp. 88-99] and [17]. This proposed scheme can be justified by Theorem 2.1 for the following mild assumption. Suppose we consider only the case
[Gp(y(s))>, 0]>∈ R(J(s)) for each singular point y(s) ∈ C, (2.22)
i.e., dim N (J(s)) = dim N (J(s)) − 1 and dim N (J(s)) ≥ 2, when detecting the bifurcation points. For this situation, the tangent vector at a singular point has a nonzero component at ˙p(s) and is expected to appear in the solution curve C of (2.12). Furthermore, as shown in Theorem 2.1, detecting the singular points of the solution curve C is equivalent to detecting the singularity of the symmetric matrix B in (2.20).
Theorem 2.1. If the condition (2.22) holds, then the following statements are equivalent: (i) rank(J(s)) ≤ 2N +1, (ii) N (J(s)) 6= {0}, and (iii) B(s) is singular. Proof. See Appendix.
2.5 Testing for stability
Along the solution curves determined by the QTICM, several types of solutions can be computed. Except for the ground and excited state solutions, we can define the so-called stable state solutions and then explain how to detect whether or not the computed solution u∗ of the DNLS equations (2.3) is stable.
It is well-known that the ground state solution of the set of DNLS equations (2.3) is found by minimizing the energy functional E(u1, u2) in (1.4b) on the
sphere: S = (u>1, u>2)> ∈ R2N|u>1u1+ u>2u2 = 1 h2 . (2.23)
On the other hand, an excited state solution is a saddle point of the energy func-tional E(u1, u2) on the sphere S. We can further define a stable state solution as
follows: note that the definition suggests that the ground state solution is also a stable state solution.
Definition 2.1. The solution u = (u>1, u>2)>∈ S of the DNLS equations (2.3) is a stable state solution if it is a local minimum of the energy functional E(u) on the sphere S.
Verifying local minimum of E(u) on S can be achieved by applying standard optimization results to the optimization problem (1.4a). We first define the La-grangian function of the optimization problem (1.4a):
L(u, λ) = E(u) − λg(u), (2.24)
where g(u) = 12 u>1u1+ u>2u2− 1/h2 is given by (2.4c). By combining the
fol-lowing two theorems, we find a practical way to identify stable sate solutions. Theorem 2.2. (Nocedal-Wright [22]) Suppose that the Karush-Kuhn-Tucker (KKT) conditions
∇uL(u∗, λ∗) = 0 (2.25a)
and
u∗ ∈ S (2.25b)
are satisfied for a certain u∗ ∈ R2N and λ∗. Suppose also that
w>∇2uuL(u∗, λ∗)w > 0, for all w>u∗= 0, w 6= 0. (2.26)
Theorem 2.3. (i) The KKT condition (2.25) is equivalent to equation (2.3). (ii) the Hessian condition (2.26) is equivalent to
Hu⊥(y) := Z>(u)H(y)Z(u) > 0, (2.27)
where Z(u) is a matrix whose columns form a basis of u⊥= {x ∈ R2N|x>u = 0} and H(y) = H(u, λ, p) = ∇2uuL(u, λ) is defined by (2.7) and (2.8).
Proof. A proof by definition is straightforward for this theorem.
Since Theorem 2.3 shows that the conditions (2.25) and (2.26) are equivalent to the equations (2.3) and (2.27), respectively, we can make the following obser-vation: Let y∗ = (u∗, λ∗, p∗) ∈ R2N +2 be a solution of DNLS equations (2.3); then the detection of the stability of E(u) in (1.4b) over S at y∗ is equivalent to the detection of the positivity of the projected Hessian matrix Hu⊥(y∗). In the
following theorem, we give a sufficient condition for determining if the projected Hessian matrix Hu⊥(y∗) is positive definite.
Theorem 2.4. Let y∗ = (u∗, λ∗, p∗) ∈ R2N +2 be a solution of DNLS equations (2.3). If a = u∗>H(y∗)u∗ 6= 0, then i " H(y∗) u∗ u∗> 0 #! = i(Hu⊥(y∗)) + (1, 1, 0),
where i(·) = (i+(·), i−(·), i0(·)) and i+(·), i−(·) and i0(·) are the number of positive,
negative, and zero eigenvalues of a symmetric matrix, respectively. Proof. See Appendix.
In short, by checking whether or not the projected Hessian matrix Hu⊥(y∗) is
positive definite, we may determine the stability of the corresponding solution u∗.
2.6 A short summary
To solve the DNLS equations (1.3), we have presented the main ideas of the QTICM from analytical and geometrical viewpoints. We have also discussed how the predictors, correctors, and bifurcation points can be computed. In addition, we have shown how the stability of the solutions can be verified.
We conclude the development of the QTICM by a remark regarding the choice of discretization schemes. While the central finite differences is used to demon-strate the idea of QTICM, we can apply other spatial discretizations in the QTICM to gain efficiency or accuracy improvements. Some possible discretizations include
(higher order) finite differences, finite volume, finite element, and pseduospectral methods. After choosing a discretization method, we need to discretize the NLS equations (1.1) to obtain the corresponding DNLS equations that are similar to (1.3). We can separate the real and imaginary parts of the DNLS equations in a similar way as shown in (2.1). The additional hyperplane constraint can then be incorporated by taking (2.15) as the normal vector of the hyperplane. As long as the computed predictors and correctors along the solution curves are accurate enough, the QTICM with the particular discretization can trace the solution curves as intended.
3
Bifurcation of a Rotating BEC with Isotropic
Po-tential Traps
We proposed the QTICM for (1.3) in Section 2. Hereafter, we focus on the DNLS equations with a harmonic trap potential in the two dimensions. In particular, we assume that:
V (x, y) = γxx2+ γyy2,
where γx and γy are the trap frequencies.
In this section, we study the bifurcation of a rotating BEC with isotropic traps in a discoidal domain, i.e. the trap frequencies are equal (γx = γy) and
Ω = {x ∈ R2| kxk2 ≤ 1}. Now, we consider the quotient solution curve C of (2.3)
by letting p = ω and fixing α = α0. For any α > 0, the ground state solution φg(x)
is unique [19] and is real radial-symmetric for ω = 0, i.e. ∂θφg(x) = 0. In fact, the
NLS equations (1.1) always have a symmetric solution φg(x) that is independent of the angular velocity. It is well-known [6, 24, 25] that when the angular velocity increases to a critical value, the bifurcations occurs in the ground state.
Similar problems also arise in the discretized system when equation (1.3) is considered for a radially symmetric trap in the discoidal domain and ω = 0. The ground state solution ug of the optimization problem
Minimize u∈RN u >Au + α 2u 2>u 2 subject to u>u = h12, (3.1)
is radial-symmetric. We make the following assumption: (H) The ground state solution of (3.1) ug ∈ N (S).
Note that the matrix S gives the discretization of the operator ∂θand the
assump-tion (H) is true for the discretizaassump-tion scheme using polar coordinate (see e.g. [18]). Using the assumption (H) and letting λg be the eigenvalue of the DNLS equations (2.1) corresponding to the eigenvector (ground state) ug/√2, ug/√2 for ω = 0,
it can been seen that the primal stalk of the solution curve of the DNLS equations (2.3) is Cp = ( 1 √ 2u g> ,√1 2u g> , λg, ω > | 0 ≤ ω ≤ ∞ ) ⊂ C. (3.2)
In this section, we prove that the curve C defined in (2.14) will have a bifurca-tion point at a finite ω = ω∗. To this end, we first define the matrices H1(z) and
H2 ∈ R2N ×2N by H1(z) + ωH2≡ " b A + 2α[[u 12 ]] 2α[[u1◦ u2]] 2α[[u1◦ u2]] A + 2α[[ub 2 2 ]] # + ω " 0 −S S 0 # = " b A + 2α[[u 12 ]] 2α[[u1◦ u2]] − ωS 2α[[u1◦ u2]] + ωS A + 2α[[ub 2 2 ]] # = H(y), (3.3)
where bA and H(y) are defined in (2.8) and (2.7), respectively.
Theorem 3.1. Suppose that 0 < α < ∞ and p = ω in (2.2), and the assumption (H) holds. Then, the primal stalk described by (3.2) has at least n = rank (S) − 3 bifurcation points at a finite valued ω = ωi∗, i = 1, . . . , n. That is, the matrix Bz(zg, ω) = " Gz aϑ a>ϑ 0 # in (2.20) is singular on Cp at ω = ω∗i, where zg = (√1 2u g>,√1 2u g>, λg)>.
Proof. From (2.6a), (2.15) and (3.3), we see that
B(zg, ω) = H Gλ aϑ G>λ 0 0 a>ϑ 0 0 ∈ R (2N +2)×(2N +2), (3.4) where Gλ = (−√12ug > , −√1 2u g> )> and aϑ = (√12ug > , −√1 2u g> )>. For ω = 0, we have H0≡ H(zg, 0) = H1(zg) = " b A + α[[ug]]2 α[[ug]]2 α[[ug]]2 A + α[[ub g]]2 # ∈ R2N ×2N (3.5a) with b A = A + α[[ug]]2− λgIN. (3.5b)
Using the fact that ug ∈ RN is the ground state of (3.1), λg is the smallest
eigenvalue of the matrix A + α[[ug]]2. On the other hand, A + α[[ug]]2 is an irreducible M-matrix and ug > 0 is the unique eigenvector corresponding to the smallest eigenvalue λg. Hence, bA in (3.5b) is positive semi-definite with a simple
zero eigenvalue. From (3.5a), it follows that λ(H0) = λ( bA + 2α[[ug]]2) ∪ λ( bA).
Here, λ(·) denotes the set of eigenvalues of a matrix. However, bA + 2α[[ug]]2 is positive definite since the ground state ug is a positive vector, implying that H
0
is positive semi-definite with one-dimensional null space. From (3.4) and since minkxk=1,x∈RN(x>, 0, 0)Bz(zg, 0)(x>, 0, 0)> = 0, using Courant-Fischer Theorem
[14, p.179] we find that the smallest third eigenvalue of Bz(zg, 0) is non-negative.
Hence, Bz(zg, 0) has 3 negative eigenvalues at most.
Next, we claim that the smallest rank (S) eigenvalues of B(zg, ω) will become negative when ω is sufficiently large. Since S is a real N × N skew-symmetric ma-trix, the eigenvalues of S can be written as ±ιµ1, . . . , ±ιµm, 0, . . . , 0, where µi > 0.
Note that m = rank (S)/2. Therefore, the eigenvalues of ωH2=
"
0 −ω ω 0
# ⊗S ∈
R2N ×2N are ±ωµ1, . . . , ±ωµm, ±ωµ1, . . . , ±ωµm, 0, . . . , 0. From (3.3) and (3.4), it
follows that lim ω→∞ 1 ωB(z g, ω) = H2 0 0 0> 0 0 0> 0 0 . (3.6)
This implies that the smallest 2m = rank (S) eigenvalues of the matrix B(zg, ω) are negative for values of ω that are sufficiently large. Also, using the fact that there are no more than 3 non-positive eigenvalues of B(zg, 0), it follows that there exist at least n = 2m − 3 bifurcation points, ωi∗ ∈ (0, ∞) i = 1, . . . , 2m − 3, such that the matrix B(zg, ωi∗) is singular. The proof is complete.
4
Bistability of a Rotating BEC
In Section 3, we analyzed bifurcation phenomena for the DNLS equations with isotropic traps (γx = γy). We found that the curve of the ground state solution
has a bifurcation point. However, it remains unclear whether bistable states exist for certain ω values. In this section we perform numerical studies by considering various settings for γx and γy using our MATLAB version of the QTICM. Finding
an S-shape-like bistability in the DNLS equations is not only a novel result to the best of the authors’ knowledge, but it would also provide a clue as to why both 0- and 1-vortex configurations can be found experimentally for a certain angular velocity, as reported in [20].
To conduct the numerical experiments, we assume that ω is the continuation parameter p. That is, we set µ0 = 0, ω0 = 0, and ν0 = 1 in (2.2). The mesh
Newton’s correction is chosen to be Tol = 10−8. To compute the initial solution of the continuation method, we set ω = 0 and use the method discussed in [9].
Our computational results are summarized in the following three subsections. In Section 4.1, we assert the existence of a bistable region by examining the bi-furcation diagram, solution profiles, and the corresponding solution stability. In particular, we find an interval [ω1, ω2] where the DNLS equations (2.3) have two
types of stable state solutions containing 0- and 1-vortex solutions. Since the trap potential and intra-component scattering length are also tunable factors in BEC experiments, we study how the trap potentials and intra-component scat-tering length affect the bistable region and the corresponding energy curves in Sections 4.2 and 4.3, respectively.
4.1 Bistability and bifurcation diagram
In this subsection, we assume that α = 100, V (x) = 1.5x2+ y2, and Ω = [−5, 5] ×
[−6, 6]. Figure 4.2 shows the bifurcation diagram of the DNLS equations (2.3) for ω ∈ (0, 0.8). The stable and unstable state solution curves are plotted with blue lines and red lines, respectively, and the corresponding nodal domains of density |u1+ ιu2|2 are attached next to the solution curves. Note that the solution curves
are defined in (2.14), the solution stability is defined in Definition 2.1, and the stability criteria are described in Theorem 2.4.
Figure 4.2-(a) shows how the solution curves bifurcate along the continuation parameter ω ∈ (0, 0.8). Starting at ω = 0, we trace the stable solution curve C0S with a 0-vortex until the first bifurcation point, which is located at ω = 0.607. The corresponding (excited) solutions of the bifurcation branch C1Ur and C1U
l contain
vortices entering the domain from right and left hand side, respectively. The solution curve bifurcates again at ω = 0.397, where a single vortex is found at the center of the domain. We can continually increase ω to obtain a stable solution curve C1S with a 1-vortex until the next bifurcation point occurs. Note that two turning points at ω = 0.520 and 0.728 are also marked.
Figure 4.2-(a) also shows that there is bistability in the DNLS equations. We denote the bistable region in terms of ω as [ω1, ω2] = [0.397, 0.607]. In this region,
the set of DNLS equations contains two stable solutions. We can denote the 0-, 0/1-and 1-vortex regions as [0, 0.397], [0.397, 0.607], 0/1-and [0.607, 0.728], respectively. These numerical results agree with the experiments reported in [20]. Madison et al. report that either the 0- or the 1-vortex configuration is found after the 0-vortex region as the angular velocity increases. After the 0/1-vortex configuration, the 1-vortex region is found.
0.397 0.520 0.607 0.728 ω 2N 1rU 1U 1S 0S (a) 3.8 4.8 4.6 4.4 4.2 4 0.397 0.607 0.728 ω ( ) E u (b)
Figure 4.2: (a) A bifurcation diagram of the solution curve C with V (x) = 1.5x2+y2 and α = 100. The blue lines indicate stable solutions. The red lines indicate unstable solutions. The corresponding nodal domains, with a density of |u1+ιu2|2,
are attached next to the solution curves. (b) Energy curves of the DNLS equations corresponding to the solution curves shown in part (a).
h ω1 CR ω2 CR Energy (ω = 0) CR
h0 0.3967274381981 - 0.6065765780083 - 4.369505165262
-h1 0.3978874922653 - 0.6059499956020 - 4.369781584647
-h2 0.3984879007864 1.62 0.6057003973605 2.27 4.369781584646 2.15
h3 0.3987605482082 1.95 0.6055944117548 2.11 4.369950507681 1.92
Table 4.1: The points ω1 and ω2 corresponding to the bistable region and the
energy for ω = 0 are computed by using different mesh size hi = 15(23)i for i =
0, 1, 2, 3. The convergence rates are presented in the columns labeled by CR.
during experiments. The number of vortices observed in the steady state depends on the initial conditions. Various initial BEC states may result in either a 0-or 1-v0-ortex. It is w0-orth mentioning that iterative numerical methods like [4, 5] can be used to solve the time-dependent problem corresponding to the DNLS equations. However, in the bistable region, these methods converge to one of the stable solutions, depending on the initial guess. Consequently, bistability can be imperceptible in these time-dependent numerical methods.
Figure 4.2-(b) demonstrates the corresponding energy curves of the bifurcation diagram shown in Figure 4.2-(a) for ω ∈ (0, 0.8). The nodal domains of the density |u1+ ιu2|2 are also included.
To assess the aforementioned numerical findings, we study how the mesh size affects the computed bifurcation points and energies. In particular, by letting
the mesh size hi = 15(23)i for i = 0, 1, 2, 3, we compute the two bifurcation points
(ω1 and ω2) associated with the bistable region, the energies for ω = 0, and the
corresponding convergence rates. As shown in Table 4.1, the QTICM converges to the energy and the bifurcation points around second order, which is parallel to the convergence rate of the central finite difference scheme used in this study. Note
that the convergence rates are computed by log(ηi−2−ηi−1
ηi−1−ηi )/log(
3
2), where ηi is the
computed values of the energy or ω corresponding to hi.
The computational results have demonstrated that the QTICM can success-fully trace the solution curves for both stable (local minimum points) and excited (saddle points) solutions. The numerical experiment results based on the method also lead to the discovery of the bistable region.
4.2 Effect of trap potential
Now, we study how the trap potential affects the bistable region. In particular, we consider the trap Vγx(x) = γxx
2 + γ
yy2 for the trap frequencies γx ∈ [1, 2]
be α = 100 and the computational domain to be Ω = [−6, 6] × [−6, 6]. We use [ω1(γx), ω2(γx)] to denote the bistable region, which is a function of γx.
Figure 4.3-(a) demonstrates the values of ω1(γx) and ω2(γx) versus γx∈ [1, 2].
The figure also plots ωc to indicate the frequency at which the energy curves
corresponding to the 0- and 1-vortices cross. The figure shows that both ω1(γx)
and ω2(γx) increase almost linearly as γx increases. In addition, the widths of
the bistable regions, i.e., |ω2(γx) − ω1(γx)|, only slightly vary as γx varies. In
Figure 4.3-(b), the energy curves for the stable state solutions are plotted for γx= 1.0, 1.2, 1.4, 1.6, 1.8, and 2. The energy curves gradually change from linear
to quadratic-like functions.
From Figure 4.3 we conclude that the trap frequency marginally affects the width of the bistable regions, while the bistable region moves to the right as γx
increases.
4.3 Effect of the intra-component scattering length
Finally, we study how the intra-component scattering length affects the bistable region. We assume the intra-component scattering length α ∈ [10, 200], the po-tential trap V (x) = 1.5x2+ y2, and computational domain Ω = [−5, 5] × [−6, 6]. We use [ω1(α), ω2(α)] to denote the bistable region, which is a function of α.
Figure 4.4-(a) demonstrates the values of ω1(α) and ω2(α) versus α ∈ [10, 200].
The figure also plots ωc(α) to indicate the frequency at which the energy curves
corresponding to 0- and 1-vortices intersect. The curves of ω1(α) and ω2(α)
de-crease exponentially. Unlike the behavior reported in Section 4.2, the widths of the bistable regions grow larger, but gradually approach a constant value as α increases. In Figure 4.4-(b), the energy curves for the stable state solutions are plotted for α = 10, 40, 80, 120, 160, and 200.
Figure 4.4 shows that the intra-component scattering length significantly af-fects the width of the bistable regions, while the bistable region moves to the left as α increases.
5
Conclusions
We have developed a quotient transformation invariant continuation method for the numerical computation of the ground and excited states of a rotating BEC. Our method circumvents the difficulties associated with transformation invariant solutions. We analytically studied the bifurcation properties of the primal stalk solution curve corresponding to the DNLS equations for an isotropic trap. In ad-dition, numerical experiments have shown that our method is reliable and finds
1.0 1.2 1.4 1.6 1.8 2.0 0.7 0.6 0.5 0.4 0.3 ω x γ 2( )x ω γ 1( )x ω γ ( ) c x ω γ (a) 0 0.2 0.4ω 0.6 0.8 0 0.2 0.4ω 0.6 0.8 0 0.2 0.4ω 0.6 0.8 0 0.2 0.4ω 0.6 0.8 0 0.2 0.4 0.6 0.8 ω 0 0.2 0.4ω 0.6 0.8 ( ) E u 4.0 3.9 3.8 3.7 ( ) E u 4.2 4.1 4.0 3.9 ( ) E u 4.4 4.3 4.2 4.1 4.0 ( ) E u 4.5 4.4 4.3 4.2 4.1 ( ) E u 4.6 4.5 4.4 4.3 4.2 ( ) E u 4.9 4.7 4.5 4.3 1.0 x γ = γ =x 1.2 1.4 x γ = γ =x 1.6 1.8 x γ = γ =x 2.0 (b)
Figure 4.3: (a) The change to the bistable region [ω1(γx), ω2(γx)] in terms of
γx ∈ [1, 2]. Here, we assume that the potential is given by V (x) = γxx2 + y2
and α = 100. The frequency corresponding to the intersection of 0- and 1-vortex energy curves is denoted by ωc(γx). (b) Energy curves containing 0- or 1-vortex
2( ) ω α 1( ) ω α 10 40 80 α120 160 200 1.0 0.9 0.8 0.7 0.6 0.5 0.4 0.3 0.2 ω ( ) c ω α (a) 0 0.2 0.4 0.6 0.8 1.0ω 0 0.2 0.4 0.6 0.8 1.0ω 0 0.2 0.4 0.6 0.8 1.0ω 0 0.2 0.4 0.6 0.8 1.0ω 0 0.2 0.4 0.6 0.8 1.0ω 0 0.2 0.4 0.6 0.8 1.0ω 1.8 1.7 1.6 1.5 ( ) E u 4.0 3.9 3.8 3.7 3.6 ( ) E u 5.6 5.4 5.2 ( ) E u 3.0 2.9 2.8 2.7 2.6 4.9 4.7 4.5 6.2 6.0 5.8 ( ) E u ( ) E u ( ) E u 10 α = α =40 80 α = α =120 160 α = α =200 (b)
Figure 4.4: (a) The change to the bistable region [ω1(α), ω2(α)] in terms of α ∈
[10, 200]. Here, we assume that the potential is given by V (x) = 1.5x2+ y2. The frequency corresponding to the intersection of 0- and 1-vortex energy curves is denoted by ωc(α). (b) The energy curves containing a 0- or 1-vortex for α = 10,
0.762
0.520 0.607
0.837
ω
2N
0.634
0.770
0.882
0.905 ω = 0.986 ω =0.898
0.923 ω =Figure 5.5: A snapshot of bifurcation diagram for ω ranges from 0.500 to 0.986.
We set V (x) = 1.5x2+ y2 and α = 100.
the bistable solutions corresponding to stable states with a 0- or 1-vortex. This result agrees with experimental observations and explains these observations nu-merically. We have also reported how the trap frequency and the intra-component scattering length affect the bistability.
Behaviors like vortex patterns and stability of rotating Bose-Einstein conden-sates in rapid rotating region have not been understood completely. As reported in [20], multiple vortices are observed as the value of ω is increased, which com-plicates the dynamics of the system. We use the QTICM to trace some solution curves as ω approaches 1. Figure 5.5 shows some of the connected solution curves that ω ranges from 0.500 to 0.986. The figure suggests that more and more bifur-cation branches can be found by the QTICM as ω increases. However, it is not clear such a path following scheme may find all the solution curves of the DNLS equations and then illustrate the complete solution behaviors. In other words, the QTICM may act as a tool to find some of the connected solutions of the DNLS equations, it remains a challenge to completely understand the dynamic system as ω increases to 1.
6
Appendix
6.1 Proof of Theorem 2.1
It is straightforward to see that statements (i) and (ii) are equivalent. Thus, we show that N (J(s)) 6= {0} is equivalent to the fact that B(s) is singular.
Necessity condition. If N (J(s)) 6= {0} then there exists a nonzero vector x ∈ R2N +1 such that J(s)x = 0. It is easy to verify that B(s)(x>, 0)> = 0 and, therefore, B(s) is singular.
Sufficiency condition. If B(s) is singular, then there exists a nonzero vector (x>, ξ)>∈ R2N +2 such that B(s) " x ξ # = 0, (6.1) which implies Gzx + ξaϑ= 0. (6.2)
Differentiating (2.11) with respect to ϑ, we have
Gu
∂u
∂ϑ(0) = 0. (6.3)
Since aϑ= (∂u∂ϑ(0)>, 0)>, using (2.6a), (6.3) becomes
Gzaϑ= 0. (6.4)
Multiplying a>ϑ from the left to (6.2), we have a>ϑGzx + (a>ϑaϑ)ξ = 0.
Since Gz is symmetric, from (6.4), we see that a>ϑGzx = 0. Hence ξ = 0. From
(6.1), it can be seen that J(s)x = 0. Namely, statement (ii) holds.
6.2 Proof of Theorem 2.4
In order to prove the Theorem 2.4, we first show some preliminary propositions and theorems.
Definition 6.1. A is congruent to B if there is a non-singular matrix U such that A = U>BU.
Proposition 6.1. If A is n × n symmetric matrix, kuk = 1 and a = u>Au 6= 0, then " A u u> 0 # is congruent to " P>AP − a1u˜˜u> 0 0> a # , where P = In− uu>, a = u>Au and ˜u = P>Au + u.
Proof. Let U1 =
"
P u u> 0
#
, where P = In− uu>. Since kuk = 1, then the
matrix P is an orthogonal projection of Rn on u⊥, i.e., P is a symmetric matrix with P2 = P. It is easy to check that U
1 is a symmetric orthogonal matrix. Then
we compute U>1 " A u u> 0 # U1 = " P>AP P>Au + u u>AP + u> u>Au # . (6.5) Let U2 = " In 0 −1 a u˜ > 1 #
where ˜u = P>Au + u. Since a = u>Au 6= 0, from (6.5) we have U>2U>1 " A u u> 0 # U1U2 = " P>AP − 1au˜˜u> 0 0> a # .
The proof is completed.
Since the eigenvalues of the symmetric matrix A are real, we adopt the con-vention that they are labeled according to non-decreasing size:
λmin= λ1 ≤ λ2≤ · · · ≤ λn−1≤ λn= λmax.
Definition 6.2. Let A be an n × n real symmetric matrix. The inertia of A is the ordered triple
i(A) = (i+(A), i−(A), i0(A))
where i+(A) is the number of positive eigenvalues of A, i−(A) is the number
of negative eigenvalues of A and i0(A) is the number of zero eigenvalues of A
(multiplicity is counted for each).
The following propositions can then be verified.
Proposition 6.2. Let A, B be n×n real symmetric matrices. Then A is congruent to B if and only if i(A) = i(B).
Proposition 6.3. Suppose A is real symmetric. Let ˜u = P>Au + u and let P = In− uu> be an orthogonal projection of Rn on u⊥. Then
a. i(P>AP − τ ˜u˜u>) = i(Z>AZ) + (0, 1, 0) if τ > 0, b. i(P>AP − τ ˜u˜u>) = i(Z>AZ) + (1, 0, 0) if τ < 0,
Proof. Let U1 = [u, u2, · · · , un] = [u, Z] be an orthogonal matrix with first column
is u, then
PU1= (I − uu>)U1 = [0, Z], (6.6)
and since u>P>Au = 0, then we have
U>1u = U˜ >1(u + PAu) = e1+ [0, w>]>, (6.7)
where w = Z>PAu ∈ Rn−1. From (6.6) and (6.7), we compute U>1(P>AP − τ ˜u˜u>)U1 = U1>P>APU1− τ U>1u˜˜u>U1
= " 0> Z> # A[0, Z] − τ " 1 w # [1, w>]. (6.8) Let U2 = " 1 −w> 0 In−1 #
and from (6.8) we have
U>2U>1(P>AP − τ ˜u˜u>)U1U2 = " −τ 0> 0 Z>AZ # (6.9) and U>2U>1(P>AP)U1U2 = " 0 0> 0 Z>AZ # . (6.10)
From (6.9), (6.10) and Proposition 6.3, it is easy to verify
i(P>AP − τ ˜u˜u>) = i(Z>AZ) + (0, 1, 0), if τ > 0, i(P>AP − τ ˜u˜u>) = i(Z>AZ) + (1, 0, 0), if τ < 0.
Proof of Theorem 2.4. Let v = u∗/ku∗k2, then it is easily seen that
i " H(y∗) u∗ u∗> 0 #! = i " H(y∗) v v> 0 #! . (6.11)
First, we consider the case a > 0. From Proposition 6.1, 6.3 and 6.3, we have
i " H(y∗) v v> 0 #! = i " P>H(y∗)P −ku∗k22 a v˜˜v > 0 0> kua∗k2 2 #! = i P>H(y∗)P −ku ∗k2 2 a v˜˜v > + (1, 0, 0)
where P = I2N− vv> and ˜v = P>H(y∗)v + v. From (6.11) and (6.12), we obtain i " H(y∗) u∗ u∗> 0 #! = i(Hu⊥(y∗)) + (1, 1, 0).
Similarly, in the case a < 0, we also have
i " H(y∗) u∗ u∗> 0 #! = i(Hu⊥(y∗)) + (1, 1, 0).
The proof is complete.
Acknowledgments
The authors are grateful to the anonymous referees for their valuable suggestions and comments. This work is partially supported by the National Science Council of Taiwan, the National Center for Theoretical Sciences, and the Taida Institute of Mathematical Sciences.
References
[1] J. R. Abo-Shaeer, C. Rarnan, J. M. Vogets, and W. Ketterle. Observation of vortex lattices in Bose-Einstein condensates. Science, 292(5516):476, 2001. [2] A. Aftalion and Q. Du. Vortices in a rotating Bose-Einstein condensate:
Critical angular velocities and energy diagrams in the Thomas-Fermi regime. Phys. Rev. A, 64:063603, 2001.
[3] E. L. Allgother and K. Gerog. Numerical path following. Handbook of Ne-merical Analysis, Edited by P. G. , 5, 1996.
[4] W. Bao. Ground states and dynamics of multi-component Bose-Einstein con-densates. SIAM Multiscale Modeling and Simulation, 2(2):210–236, 2004. [5] W. Bao and Q. Du. Computing the ground state solution of Bose-Einstein
con-densates by a normalized gradient flow. SIAM J. Sci. Comput., 25(5):1674– 1697, 2004.
[6] W. Bao, H. Wang, and P. Markowich. Ground, symmetric and central vertix states in rotatiog Bose-Einstein condensates. Comm. Math. Sci., 3(1):57–88, 2005.
[7] S.-L. Chang and C.-S. Chien. Adaptive continuation algorithms for computing energy levels of rotating Bose-Einstein condensate. Comput. Phys. Comm., 177:707–719, 2007.
[8] 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. J. Comp. Phys., 210:439–458, 2005.
[9] S.-M. Chang, W.-W. Lin, and S.-F. Shieh. Gauss-Seidel-type Methods for Energy States of A Multi-Component Bose-Einstein Condensate. J. Comp. Phy., 202:367–390, 2005.
[10] W. Govaerts. Stable solvers and block elimination for bordered systems. SIAM J. Matrix Anal. Appl., 12:469–483, 1991.
[11] W. Govaerts. Numerical Methods for Bifurcations of Dynamical Equilibria Computation,. SIAM, Philadelphia, 2000.
[12] W. Govaerts and J.D. Pryce. Mixed block elimination for linear systems with wider borders. IMA J. Numer. Anal., 13:161–180, 1993.
[13] D. S. Hall, M. R. Matthews, J. R. Ensher, C. E. Wieman, and E. A. Cor-nell. Dynamics of component separation in a binary mixture of Bose-Einstein condensates. Phys. Rev. Lett., 81:1539–1542, 1998.
[14] R. A. Horn and C. R. Van Johnson. Matrix Analysis. Cambridge University Press, 1990.
[15] H. B. Keller. Lectures on Numerical Methods in Bifurcation Problems. Springer-Verlag, Berlin, 1987.
[16] Y.-C. Kuo, W.-W. Lin, and S.-F. Shieh. Bifurcation Analysis of a Two-Component Bose-Einstein Condensate. Physica D, 211:311–346, November 2005.
[17] Y.-C. Kuo, W.-W. Lin, S.-F. Shieh, and W. Wang. A hyperplane-constrained continuation method for near singularity in coupled nonlinear Schr¨odinger equations. Applied Numerical Mathematics, 60(5):513–526, 2010.
[18] M.-C. Lai. A note on finite difference discretizations for Poisson equation on a disk. Numerical Methods for Partial Differential Equations, 17(3):199–203, 2001.
[19] 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.
[20] K. W. Madison, F. Chevy, W. Wohlleben, and J. Dalibard. Vortex formation in a stirred bose-einstein condensate. Phys. Rev. Lett., 84(5):806–809, Jan 2000.
[21] K. W. Madison, F. Chevy, W. Wohlleben, and J. Dalibard. Vortices in a stirred Bose-Einstein condensate. Journal of Modern Optics, 47(14-15):2715– 2724, 2000.
[22] J. Nocedal and S. J. Wright. Numerical Optimization. Springer-Verlag, New York, 2 edition, 2006.
[23] L. P. Pitaevskii and S. Stringari. Bose-Einstein condensation. Clarendon Press, 2003.
[24] R. Seiringer. Gross-pitaevskii theory of the rotating Bose gas. Communica-tions in Mathematical Physics, 229(3):491–509, September 2002.
[25] R. Seiringer. Ground state asymptotice of a dilute rotation gas. J. Phy. A, 36:9755–9778, 2003.
![Figure 4.3: (a) The change to the bistable region [ω 1 (γ x ), ω 2 (γ x )] in terms of](https://thumb-ap.123doks.com/thumbv2/9libinfo/7412333.103988/20.892.231.656.257.837/figure-change-bistable-region-ω-γ-ω-terms.webp)
![Figure 4.4: (a) The change to the bistable region [ω 1 (α), ω 2 (α)] in terms of α ∈](https://thumb-ap.123doks.com/thumbv2/9libinfo/7412333.103988/21.892.231.653.235.855/figure-change-bistable-region-ω-α-ω-terms.webp)
