Ground states and their characterization of spin-$F$ Bose-Einstein condensates

TL;DR

This paper introduces a novel optimization-based numerical method for computing ground states of spin-F Bose-Einstein condensates, including complex spin-3 cases, using manifold optimization and multigrid acceleration.

Contribution

It develops the first applicable algorithm for arbitrary spin BECs, employing Fourier pseudospectral discretization, manifold optimization, and multigrid techniques.

Findings

Efficient computation of ground states for spin-1, spin-2, and spin-3 BECs.

Demonstration of interesting physical phenomena in various dimensions.

Validation of the method's effectiveness through extensive numerical experiments.

Abstract

The computation of the ground states of spin-$F$ Bose-Einstein condensates (BECs) can be formulated as an energy minimization problem with two quadratic constraints. We discretize the energy functional and constraints using the Fourier pseudospectral schemes and view the discretized problem as an optimization problem on manifold. Three different types of retractions to the manifold are designed. They enable us to apply various optimization methods on manifold to solve the problem. Specifically, an adaptive regularized Newton method is used together with a cascadic multigrid technique to accelerate the convergence. According to our limited knowledege, our method is the first applicable algorithm for BECs with an arbitrary integer spin, including the complicated spin-3 BECs. Extensive numerical results on ground states of spin-1, spin-2 and spin-3 BECs with diverse interaction and optical lattice potential in one/two/three dimensions are reported to show the efficiency of our method and to demonstrate some interesting physical phenomena.

Tables

Table 1. Table 1: Numerical results of spin-1 BEC in 2D

	ARNT				RGBB				RTR
retr	f	nrmG	iter	time	f	nrmG	iter	time	f	nrmG	iter	time
$M=0.0$
R1	15.1032	7.8e-07	4 (38)	17.0	15.1032	8.3e-07	239	18.1	15.1032	6.2e-07	20 (17)	31.6
R2	15.1032	7.7e-07	4 (38)	15.3	15.1032	6.3e-07	257	15.1	15.1032	6.3e-07	20 (17)	31.0
R3	15.1032	4.4e-07	4 (35)	15.9	15.1032	3.3e-07	255	15.9	15.1032	6.2e-07	20 (17)	28.9
$M=0.2$
R1	15.1411	6.6e-07	4 (44)	19.9	15.1411	9.3e-07	258	21.5	15.1411	2.6e-07	21 (26)	54.9
R2	15.1411	8.1e-07	4 (42)	18.6	15.1411	1.0e-06	254	15.1	15.1411	2.7e-07	21 (26)	51.2
R3	15.1411	7.6e-07	4 (42)	18.6	15.1411	7.9e-07	261	18.2	15.1411	2.6e-07	21 (26)	50.8
$M=0.5$
R1	15.3436	6.7e-07	4 (64)	27.2	15.3436	8.1e-07	431	36.4	15.3436	2.1e-07	21 (31)	61.2
R2	15.3436	3.5e-07	4 (67)	24.9	15.3436	9.1e-07	421	26.4	15.3436	1.9e-07	21 (31)	61.4
R3	15.3436	5.3e-07	4 (64)	25.8	15.3436	9.0e-07	429	26.9	15.3436	1.9e-07	21 (31)	63.5
$M=0.9$
R1	15.9621	7.7e-07	4 (51)	23.8	15.9621	8.8e-07	323	28.4	15.9621	3.9e-07	21 (26)	54.3
R2	15.9621	7.5e-07	4 (52)	21.7	15.9621	1.0e-06	279	18.1	15.9621	3.9e-07	21 (26)	51.3
R3	15.9621	5.3e-07	4 (55)	23.5	15.9621	7.6e-07	483	31.8	15.9621	3.6e-07	21 (26)	51.6

Table 2. Table 2: Numerical results of spin-1 BEC in 3D

	ARNT				RGBB				RTR
retr	f	nrmG	iter	time	f	nrmG	iter	time	f	nrmG	iter	time
$M=0.0$
R1	55.4362	6.9e-07	4 (34)	1187.6	55.4362	1.0e-06	720	7197.9	55.4362	7.2e-07	17 (15)	1484.0
R2	55.4362	7.5e-07	4 (34)	1042.3	55.4362	9.5e-07	384	2872.6	55.4362	7.2e-07	17 (15)	1306.8
R3	55.4362	5.6e-07	4 (35)	1210.5	55.4362	7.0e-07	227	1169.3	55.4362	7.0e-07	17 (15)	1540.0
$M=0.2$
R1	55.4362	4.2e-07	4 (40)	1568.7	55.4363	9.0e-05	10000	47615.4	55.4362	7.3e-07	17 (15)	1523.8
R2	55.4362	5.8e-07	4 (33)	2041.6	55.4363	9.8e-05	10000	42311.5	55.4362	7.1e-07	17 (15)	1309.3
R3	55.4362	3.2e-07	4 (41)	2626.7	55.4363	9.1e-05	10000	48975.6	55.4362	7.1e-07	17 (15)	1538.8
$M=0.5$
R1	55.4362	5.6e-07	4 (41)	1395.6	55.4362	8.3e-07	201	1242.1	55.4362	7.0e-07	17 (15)	1514.3
R2	55.4362	3.2e-07	4 (41)	1392.8	55.4362	3.1e-07	310	1846.4	55.4362	7.0e-07	17 (15)	1309.3
R3	55.4362	5.7e-07	4 (36)	1394.5	55.4362	1.0e-06	295	1603.6	55.4362	7.2e-07	17 (15)	1532.3
$M=0.9$
R1	55.4362	5.4e-07	4 (39)	1421.1	55.4362	9.8e-07	444	4093.6	55.4362	6.8e-07	17 (15)	1527.9
R2	55.4362	2.9e-07	4 (41)	1323.6	55.4362	8.7e-07	410	2916.9	55.4362	7.2e-07	17 (15)	1297.3
R3	55.4362	6.8e-07	4 (34)	1350.7	55.4362	9.7e-07	236	1460.1	55.4362	6.9e-07	17 (15)	1548.9

Table 3. Table 3: Numerical results of spin-2 BEC in 2D

	ARNT				RGBB				RTR
retr	f	nrmG	iter	time	f	nrmG	iter	time	f	nrmG	iter	time
$M=0.0$
R1	14.3386	5.9e-07	4 (32)	8.3	14.3386	9.0e-07	280	12.1	14.3386	9.3e-07	17 (18)	18.3
R2	14.3386	5.9e-07	4 (32)	8.4	14.3386	8.0e-07	238	9.8	14.3386	9.3e-07	17 (18)	17.2
R3	14.3386	5.9e-07	4 (32)	8.2	14.3386	9.0e-07	225	8.3	14.3386	9.3e-07	17 (18)	17.2
$M=0.5$
R1	14.3730	3.8e-07	4 (59)	13.5	14.3730	1.3e-07	346	14.9	14.3730	4.7e-07	18 (22)	25.4
R2	14.3730	3.8e-07	4 (59)	13.8	14.3730	9.9e-07	323	13.9	14.3730	4.7e-07	18 (22)	24.3
R3	14.3730	3.8e-07	4 (59)	13.5	14.3730	9.6e-08	523	32.1	14.3730	4.7e-07	18 (22)	23.1
$M=1.5$
R1	14.6754	4.3e-07	4 (62)	13.5	14.6754	9.5e-07	462	19.0	14.6754	2.8e-07	18 (28)	30.5
R2	14.6754	4.3e-07	4 (62)	14.2	14.6754	9.9e-07	519	31.8	14.6754	2.8e-07	18 (28)	28.5
R3	14.6754	4.3e-07	4 (62)	13.8	14.6754	8.5e-07	402	19.4	14.6754	2.8e-07	18 (28)	28.5

Table 4. Table 4: Numerical results of spin-2 BEC in 3D

	ARNT				RGBB				RTR
retr	f	nrmG	iter	time	f	nrmG	iter	time	f	nrmG	iter	time
$M=0.0$
R1	46.2917	5.9e-07	4 (36)	2179.1	46.2917	4.7e-07	358	4192.9	46.2917	3.6e-07	18 (15)	2905.1
R2	46.2917	6.6e-07	4 (40)	2163.8	46.2917	9.7e-07	346	3667.8	46.2917	3.4e-07	18 (15)	2720.1
R3	46.2917	6.5e-07	4 (45)	2534.3	46.2917	7.7e-07	356	4219.3	46.2917	2.7e-07	18 (15)	2937.6
$M=0.5$
R1	45.7403	6.2e-07	5 (81)	4824.5	45.7403	9.3e-07	1129	9464.9	45.7403	5.5e-07	19 (29)	5432.8
R2	45.7403	8.4e-07	5 (70)	4176.5	45.7403	9.2e-07	1050	7630.2	45.7403	4.7e-07	19 (28)	5277.9
R3	45.7403	7.7e-07	5 (76)	4511.5	45.7403	9.6e-07	1296	12903.3	45.7403	5.3e-07	19 (28)	5483.4
$M=1.5$
R1	46.8619	7.2e-07	4 (56)	3032.4	46.8619	9.8e-07	391	3182.1	46.8619	3.3e-07	19 (22)	4371.3
R2	46.8619	6.6e-07	4 (56)	3038.1	46.8619	8.9e-07	380	2996.6	46.8619	3.4e-07	19 (22)	4142.4
R3	46.8619	5.8e-07	4 (55)	3151.1	46.8619	9.3e-07	515	6125.7	46.8619	3.5e-07	19 (22)	4370.7

Table 5. Table 5: Numerical results of spin-3 BEC in 2D

	ARNT				RGBB				RTR
retr	f	nrmG	iter	time	f	nrmG	iter	time	f	nrmG	iter	time
$M=0.0$
R1	11.8279	6.8e-07	4 (58)	29.8	11.8279	6.8e-07	634	50.3	11.8279	5.2e-07	28 (159)	359.0
R2	11.8279	7.1e-07	4 (49)	24.4	11.8279	8.6e-07	469	39.5	11.8279	1.4e-07	18 (58)	83.2
R3	11.8279	6.8e-07	4 (58)	25.9	11.8279	8.6e-07	496	32.9	11.8279	5.9e-07	32 (184)	476.4
$M=0.5$
R1	11.8334	4.9e-07	4 (75)	32.7	11.8334	7.8e-07	577	48.0	11.8334	1.8e-07	18 (79)	124.2
R2	11.8334	5.2e-07	4 (75)	35.4	11.8334	8.7e-07	472	37.7	11.8334	1.7e-07	18 (79)	127.5
R3	11.8334	5.9e-07	4 (73)	31.0	11.8334	9.5e-07	572	40.0	11.8334	1.9e-07	18 (42)	63.5
$M=1.5$
R1	11.8780	4.0e-07	5 (134)	51.6	11.8780	1.0e-06	2874	190.9	11.8780	3.0e-07	18 (51)	78.6
R2	11.8780	3.6e-07	6 (150)	79.6	11.8780	9.7e-07	857	74.9	11.8780	2.7e-07	18 (61)	86.9
R3	11.8780	6.2e-07	4 (95)	36.5	11.8780	9.9e-07	3123	221.6	11.8780	3.0e-07	18 (49)	69.3

Table 6. Table 6: Numerical results of spin-3 BEC in 3D

	ARNT				RGBB				RTR
retr	f	nrmG	iter	time	f	nrmG	iter	time	f	nrmG	iter	time
$M=0.0$
R1	42.9752	5.1e-07	5 (164)	4243.2	42.9752	9.6e-07	496	1328.1	42.9752	9.5e-08	18 (23)	794.0
R2	42.9752	5.9e-07	4 (92)	4405.6	42.9774	1.6e-04	10000	15255.9	42.9752	9.9e-08	18 (30)	1027.9
R3	42.9752	9.9e-07	40 (154)	16610.8	42.9767	2.2e-03	10000	14740.1	42.9752	4.7e-07	48 (98)	8453.7
$M=0.5$
R1	42.9822	8.8e-07	4 (115)	1050.1	42.9822	9.4e-07	864	1370.7	42.9822	7.7e-07	96 (159)	27219.4
R2	42.9822	8.8e-07	4 (115)	1050.4	42.9822	1.0e-06	1150	1956.0	42.9822	6.9e-07	21 (40)	1601.1
R3	42.9822	8.9e-07	4 (115)	1013.1	42.9822	1.0e-06	1328	2128.7	42.9822	6.7e-07	21 (40)	1539.8
$M=1.5$
R1	43.0399	2.8e-07	5 (151)	1959.2	43.0400	4.2e-03	10000	15199.3	43.0399	1.2e-07	82 (379)	55201.7
R2	43.0399	2.3e-07	5 (156)	2100.3	43.0399	2.1e-03	10000	15261.4	43.0417	1.6e-07	22 (36)	1479.3
R3	43.0399	2.5e-07	5 (151)	1683.8	43.0399	9.1e-07	1162	1802.2	43.0417	1.6e-07	22 (37)	1485.9

Table 7. Table 7: Ground state energies of spin-2 BEC

$M$	1D			2D			3D
	Case I	Case II	Case III	Case I	Case II	Case III	Case I	Case II	Case III
0.0	10.3700	25.6185	24.1144	9.5754	14.3386	13.9158	39.0045	46.9770	46.2917
0.5	10.3700	25.7372	22.9404	9.5754	14.3730	13.5746	39.0045	47.0301	45.7403
1.5	10.3700	26.8415	25.4640	9.5754	14.6754	14.2734	39.0045	47.5117	46.8619

Table 8. Table 8: Ground state energies of spin-3 BEC

$M$	1D		2D		3D
	Case I	Case II	Case I	Case II	Case I	Case II
0.0	17.1091	17.2527	11.7811	11.8279	42.9028	42.9752
0.5	17.1289	17.2693	11.7877	11.8334	42.9115	42.9822
1.5	17.2905	17.4034	11.8415	11.8780	42.9825	43.0399

Equations

$$E(\Phi (\cdot ))$$

$$F_x={\Phi }^{\ast }{\mathrm{f}}_x\Phi ,F_y={\Phi }^{\ast }{\mathrm{f}}_y\Phi ,F_z={\Phi }^{\ast }{\mathrm{f}}_z\Phi ,$$

$${\mathrm{f}}_x=\frac{1}{\sqrt{2}}\left(\begin{array}{ccc}0& 1& 0\\ 1& 0& 1\\ 0& 1& 0\end{array}\right),{\mathrm{f}}_y=\frac{i}{\sqrt{2}}\left(\begin{array}{ccc}0& -1& 0\\ 1& 0& -1\\ 0& 1& 0\end{array}\right),{\mathrm{f}}_z=\left(\begin{array}{ccc}1& 0& 0\\ 0& 0& 0\\ 0& 0& -1\end{array}\right),$$

$$\begin{array}{rl}& F_x=\frac{1}{\sqrt{2}}\left[{\overline{\varphi }}_1{\varphi }_0+{\overline{\varphi }}_0({\varphi }_1+{\varphi }_{-1})+{\overline{\varphi }}_{-1}{\varphi }_0\right],F_y=\frac{i}{\sqrt{2}}\left[-{\overline{\varphi }}_1{\varphi }_0+{\overline{\varphi }}_0({\varphi }_1-{\varphi }_{-1})+{\overline{\varphi }}_{-1}{\varphi }_0\right].\end{array}$$

$$E(\Phi (\cdot ))$$

$${\mathrm{f}}_x=\left(\begin{array}{ccccc}0& 1& 0& 0& 0\\ 1& 0& \sqrt{\frac{3}{2}}& 0& 0\\ 0& \sqrt{\frac{3}{2}}& 0& \sqrt{\frac{3}{2}}& 0\\ 0& 0& \sqrt{\frac{3}{2}}& 0& 1\\ 0& 0& 0& 1& 0\end{array}\right),{\mathrm{f}}_y=i\left(\begin{array}{ccccc}0& -1& 0& 0& 0\\ 1& 0& -\sqrt{\frac{3}{2}}& 0& 0\\ 0& \sqrt{\frac{3}{2}}& 0& -\sqrt{\frac{3}{2}}& 0\\ 0& 0& \sqrt{\frac{3}{2}}& 0& -1\\ 0& 0& 0& 1& 0\end{array}\right)$$

$${\mathrm{f}}_z=\text{diag}(2,1,0,-1,-2).$$

$$\begin{array}{rl}& F_x={\overline{\varphi }}_2{\varphi }_1+{\overline{\varphi }}_1{\varphi }_2+{\overline{\varphi }}_{-1}{\varphi }_{-2}+{\overline{\varphi }}_{-2}{\varphi }_{-1}+\frac{\sqrt{6}}{2}({\overline{\varphi }}_1{\varphi }_0+{\overline{\varphi }}_0{\varphi }_1+{\overline{\varphi }}_0{\varphi }_{-1}+{\overline{\varphi }}_{-1}{\varphi }_0),\\ & F_y=i\left[{\overline{\varphi }}_1{\varphi }_2-{\overline{\varphi }}_2{\varphi }_1+{\overline{\varphi }}_{-2}{\varphi }_{-1}-{\overline{\varphi }}_{-1}{\varphi }_{-2}+\frac{\sqrt{6}}{2}({\overline{\varphi }}_0{\varphi }_1-{\overline{\varphi }}_1{\varphi }_0+{\overline{\varphi }}_{-1}{\varphi }_0-{\overline{\varphi }}_0{\varphi }_{-1})\right],\\ & F_z=2|{\varphi }_2{|}^2+|{\varphi }_1{|}^2-|{\varphi }_{-1}{|}^2-2|{\varphi }_{-2}{|}^2.\end{array}$$

$$\mathbf{A}=\frac{1}{\sqrt{5}}\left(\begin{array}{ccccc}0& 0& 0& 0& 1\\ 0& 0& 0& -1& 0\\ 0& 0& 1& 0& 0\\ 0& -1& 0& 0& 0\\ 1& 0& 0& 0& 0\end{array}\right),$$

$$A_{00}=\frac{1}{\sqrt{5}}(2{\varphi }_2{\varphi }_{-2}-2{\varphi }_1{\varphi }_{-1}+{\varphi }_0^2).$$

$$E(\Phi (\cdot ))$$

$$\displaystyle\qquad\qquad+\frac{\beta_{3}}{2}\sum_{l=-2}^{2}|A_{2l}|^{2}\bigg{\}}\;d{\bf x},$$

$${\mathrm{f}}_x=\left(\begin{array}{ccccccc}0& \sqrt{3/2}& 0& 0& 0& 0& 0\\ \sqrt{3/2}& 0& \sqrt{5/2}& 0& 0& 0& 0\\ 0& \sqrt{5/2}& 0& \sqrt{3}& 0& 0& 0\\ 0& 0& \sqrt{3}& 0& \sqrt{3}& 0& 0\\ 0& 0& 0& \sqrt{3}& 0& \sqrt{5/2}& 0\\ 0& 0& 0& 0& \sqrt{5/2}& 0& \sqrt{3/2}\\ 0& 0& 0& 0& 0& \sqrt{3/2}& 0\end{array}\right),$$

$${\mathrm{f}}_y=i\left(\begin{array}{ccccccc}0& -\sqrt{3/2}& 0& 0& 0& 0& 0\\ \sqrt{3/2}& 0& -\sqrt{5/2}& 0& 0& 0& 0\\ 0& \sqrt{5/2}& 0& -\sqrt{3}& 0& 0& 0\\ 0& 0& \sqrt{3}& 0& -\sqrt{3}& 0& 0\\ 0& 0& 0& \sqrt{3}& 0& -\sqrt{5/2}& 0\\ 0& 0& 0& 0& \sqrt{5/2}& 0& -\sqrt{3/2}\\ 0& 0& 0& 0& 0& \sqrt{3/2}& 0\end{array}\right)$$

$${\mathrm{f}}_z=\text{diag}(3,2,1,0,-1,-2,-3).$$

$$\begin{array}{rl}& F_{+}=F_x+i{F}_y=\sqrt{6}{\overline{\psi }}_3{\psi }_2+\sqrt{10}{\overline{\psi }}_2{\psi }_1+2\sqrt{3}{\overline{\psi }}_1{\psi }_0+2\sqrt{3}{\overline{\psi }}_0{\psi }_{-1}+\sqrt{10}{\overline{\psi }}_{-1}{\psi }_{-2}+\sqrt{6}{\overline{\psi }}_{-2}{\psi }_{-3},\\ & F_z=3|{\psi }_3{|}^2+2|{\psi }_2{|}^2+|{\psi }_1{|}^2-|{\psi }_{-1}{|}^2-2|{\psi }_{-2}{|}^2-3|{\psi }_{-3}{|}^2.\end{array}$$

$$\mathbf{A}=\frac{1}{\sqrt{7}}\left(\begin{array}{ccccccc}0& 0& 0& 0& 0& 0& 1\\ 0& 0& 0& 0& 0& -1& 0\\ 0& 0& 0& 0& 1& 0& 0\\ 0& 0& 0& -1& 0& 0& 0\\ 0& 0& 1& 0& 0& 0& 0\\ 0& -1& 0& 0& 0& 0& 0\\ 1& 0& 0& 0& 0& 0& 0\end{array}\right),{\mathbf{A}}_0=\frac{1}{\sqrt{7}}\left(\begin{array}{ccccccc}0& 0& 0& 0& 0& 0& \frac{5}{2\sqrt{3}}\\ 0& 0& 0& 0& 0& 0& 0\\ 0& 0& 0& 0& \frac{-\sqrt{3}}{2}& 0& 0\\ 0& 0& 0& \sqrt{\frac{2}{3}}& 0& 0& 0\\ 0& 0& \frac{-\sqrt{3}}{2}& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& 0& 0\\ \frac{5}{2\sqrt{3}}& 0& 0& 0& 0& 0& 0\end{array}\right)$$

$${\vec{a}}_{\pm 1}=\frac{1}{\sqrt{7}}{\left(\frac{5}{2\sqrt{3}},-\frac{\sqrt{5}}{2},\frac{1}{\sqrt{6}},\frac{1}{\sqrt{6}},-\frac{\sqrt{5}}{2},\frac{5}{2\sqrt{3}}\right)}^T,{\vec{a}}_{\pm 2}=\frac{1}{\sqrt{7}}{\left(\sqrt{\frac{5}{6}},-\sqrt{\frac{5}{3}},\sqrt{2},-\sqrt{\frac{5}{3}},\sqrt{\frac{5}{6}}\right)}^T.$$

$$A_{00}=$$

$$A_{20}=$$

$$A_{2,\pm 1}=$$

$$A_{2,\pm 2}=$$

$$N(\Phi (\cdot )):=\Vert \Phi (\cdot ){\Vert }^2={\int }_{{\mathbb{R}}^d}\sum _{l=-F}^F|{\varphi }_l(\mathbf{x}){|}^{2}d\mathbf{x}=1,$$

$$M(\Phi (\cdot )):={\int }_{{\mathbb{R}}^d}\sum _{l=-F}^{F}l|{\varphi }_l(\mathbf{x}){|}^{2}d\mathbf{x}=M.$$

$$E_g:=E\left({\Phi }_g\right)=\underset{\Phi \in S_M}{\min }E\left(\Phi \right),$$

$$S_M=\left\{\Phi =({\varphi }_F,\dots ,{\varphi }_{-F}{)}^T\in {\mathcal{C}}^{2F+1}|\Vert \Phi \Vert =1,{\int }_{{\mathbb{R}}^d}\sum _{l=-F}^{F}l|{\varphi }_l(\mathbf{x}){|}^2=M,E(\Phi )<\infty \right\}.$$

$$E(\Phi )$$

$$\partial_{xx}^{f}\phi\Big{|}_{j}=-\frac{1}{n}\sum_{p=-n/2}^{n/2-1}\lambda_{p}^{2}\tilde{\phi}_{pl}e^{i\frac{2\pi jp}{n}},$$

$${\tilde{\varphi }}_{pl}=\sum _{j=0}^{n-1}{\varphi }_l(x_j)e^{-i\frac{2\pi jp}{n}},{\lambda }_p=\frac{2\pi p}{b-a},p=-\frac{n}{2},\cdots ,\frac{n}{2}-1.$$

$$E_h(X)$$

$$\rho =\mathrm{diag}(X{X}^{\ast }),{\mathrm{F}}_{\alpha }=\mathrm{diag}(X{f}_{\alpha }^T{X}^{\ast }),\alpha =x,y,z,$$

$${\mathrm{A}}_{00}=\mathrm{diag}(X\mathbf{A}{X}^T),{\mathrm{A}}_{2l}=\mathrm{diag}(X{\mathbf{A}}_l{X}^T),l=-2,\cdots ,2$$

$$N(\Phi (\cdot ))\approx h\sum _{j=0}^{n-1}\sum _{l=-3}^3|{\varphi }_{jl}{|}^2=\mathrm{tr}(X^{\ast }X)=1,$$

Full text

{AMS}

35Q55, 35A01, 81Q99

1 Introduction

Bose-Einstein condensate (BEC), first predicted by A. Einstein based on S. N. Bose’s work, refers to the state of matter in which part of the bosons occupy the same quantum state at extremely low temperature. The earliest experimental observations of BEC were announced in 1995 [5, 14, 16] and have attracted numerous researchers into the study of condensates of dilute gases ever since [4, 15, 17, 23, 25, 27, 28]. While in early experiments the spin degrees of freedom are frozen due to the magnetic trapping, the experimental realizations of spin-1 and spin-2 condensates have been achieved later by optical confinements [12, 19, 24, 29, 31] and revealed various exciting phenomena absent in single-component condensates.

Numerous theoretical studies of spinor condensates have been carried out after the experimental achievement [20, 22, 26, 30]. At zero temperature, a spin-$F$ ($F=1,2,\ldots$) BEC is described by a $2F+1$ vector wave function $\Phi=(\phi_{F},\cdots,\phi_{-F})^{T}\in\mathcal{C}^{2F+1}$ and a generalized coupled Gross-Pitaevskii equation (GPE). Three important invariants of it are the mass of the wave function, the magnetization and the energy per particle. A fundamental problem in BEC is to find the condensate stationary states, which is obtained by minimizing the Gross-Pitaevskii (GP) energy functional subject to the conservation of total mass and magnetization.

Different numerical methods have been proposed in the literature to compute the ground state of a spin-1 BEC [11, 7, 34, 35, 9]. Among them, a very popular method is the imaginary time method combined with a proper discretization scheme to evolve the resulted gradient flow equation under the normalization of the wave function [6, 8, 9, 11]. To apply the normalized gradient flow method to compute the ground state of a spin-$F$ BEC, $2F+1$ projection constants have to be determined in the normalization step, while only two normalization conditions (i.e., the two constraints) are given. In the literature, this method is applied to compute the ground state of a spin-1 BEC through the introduction of a random variable [34, 35] or a third normalization condition [9]. Recently, a projection gradient method [11, 32] has been proposed to compute ground states of spin-1 and spin-2 BEC, where a continuous normalized gradient flow (CNGF) was discretized by the Crank-Nicolson finite difference (CNFD) method with a proper and very special way to deal with the nonlinear terms. This scheme is proved to be mass- and magnetization-conservative and energy-diminishing in the discretized level. However, a fully nonlinear coupled system has to be solved at each time step.

Most of the existing numerical methods for computing the ground states of spinor BEC evolve from the gradient flow method, and thus converge at most linearly and/or require to solve a large scale linear system per iteration, which leads to quite expensive computational cost. Most of them are specially designed for spin-1 or spin-2 BEC, but the spin-3 cases are rarely discussed. Meanwhile, over the last decade, some advanced optimization methods have been developed for solving minimization problems on matrix manifolds, such as the Riemannian Newton methods and trust-region methods [2, 1] with superlinear or quadratic convergence rate. The aim of this paper is to explore a new way to compute the ground states of spinor BEC, and propose an efficient regularized Newton method for the general spin-$F$ cases. We first discretize the energy functional and the constraints with the Fourier pseudospectral schemes and thus approximate the original infinite dimensional problem by a finite dimensional minimization problem, of which the feasible region can be proven to be a Riemannian manifold. We give the formulas of Riemannian gradient and Hessian on this manifold, and then aim to apply an adaptive regularized Newton method to solve the Riemannian optimization problem. To improve the efficiency and stability, we adopt the cascadic multigrid technique and use a Riemannian gradient method with Barzilai-Borwein step size to compute initial points on each mesh. Three different retractions on the manifold are proposed for the implementation of Riemannian optimization algorithms. The first one is the classical projective retraction, and the second one comes from the normalized gradient flow [9]. The computation of them relies on finding a unique zero of a single-variable function, which can be done quite efficiently and accurately. The third retraction is proposed as an approximation of the first one, with a very brief closed-form formula. Extensive numerical experiments demonstrate that our approach can quickly compute an accurate approximation of the ground state, and is more stable than the classical Riemannian trust-region method. The algorithm remains effective even for the complicated spin-3 BEC in 3D with an optical lattice potential, for which there exists no applicable algorithm before.

The rest of this paper is organized as follows. Specific problem statements of spin-1, spin-2 and spin-3 BEC are given in Section 2. Discretizations of the energy functional and the constraints via the Fourier pseudospectral schemes are introduced in Section 3. In Section 4, we give some preliminaries on Riemannian optimization, and investigate the manifold structure of the feasible region. In Section 5, we present a modified version of the adaptive regularized Newton method for solving the discretized optimization problem. The three retractions are described in Section 6, and detailed numerical results are reported in Section 7 to illustrate the efficiency and accuracy of our algorithm. Finally, some conclusions are given in Section 8.

2 Problem Statement

The specific formulation of the minimization problem for computing the ground states of spin-1, spin-2 and spin-3 BEC is stated as follows:

Spin-1. For a spin-1 BEC, the GP energy functional for the spin-1 wave function is given by

[TABLE]

where ${\bf x}=x$ in 1D, ${\bf x}=(x,y)^{T}$ in 2D and ${\bf x}=(x,y,z)^{T}$ in 3D, $V({\bf x})$ is the external confining potential, $p$ and $q$ are the linear and quadratic Zeeman energy shifts, respectively. $\beta_{0}$ is the density dependent interaction strength between the particles and $\beta_{1}$ is the spin dependent interaction strength, $\mathbf{F}:=\mathbf{F}(\Phi)=(F_{x},F_{y},F_{z})^{T}\in\mathcal{R}^{3}$ is the spin vector given by

[TABLE]

where $\Phi^{*}=\overline{\Phi}^{T}$ is the conjugate transpose and $\mathrm{f}_{\alpha}$ ($\alpha=x,y,z$) are the 3-by-3 spin-1 matrices

[TABLE]

and $i=\sqrt{-1}$ is the imaginary unit. In detail, the components of spin vector $\mathbf{F}$ can be written explicitly as: $F_{z}=|\phi_{1}|^{2}-|\phi_{-1}|^{2}$,

[TABLE]

Spin-2. For a spin-2 BEC, the GP energy is given by

[TABLE]

where $\beta_{2}$ is the spin-singlet interaction strength and all the other parameters $p,q,\beta_{0},\beta_{1}$ are the same as those in the spin-1 case, $\mathbf{F}:=\mathbf{F}(\Phi)=(F_{x},F_{y},F_{z})^{T}\in\mathcal{R}^{3}$ is the spin vector defined by Eq. 2.2, with $\mathrm{f}_{\alpha}$ ($\alpha=x,y,z$) given by the 5-by-5 spin-2 matrices

[TABLE]

and

[TABLE]

Therefore, the spin vector $\mathbf{F}$ can be written explicitly

[TABLE]

Define the matrix

[TABLE]

then $A_{00}:=A_{00}(\Phi)=\Phi^{T}\mathbf{A}\Phi$ can be expressed as

[TABLE]

Spin-3. For a spin-3 BEC, the GP energy is given by

[TABLE]

where $\beta_{3}$ is the spin-quintet interaction strength, and all the other parameters $p,q,\beta_{0},\beta_{1},\beta_{2}$ are the same as those in the spin-1,2 cases. $\mathbf{F}:=\mathbf{F}(\Phi)=(F_{x},F_{y},F_{z})^{T}\in\mathcal{R}^{3}$ is the spin vector defined by (2.2), with $\mathrm{f}_{\alpha}$ ($\alpha=x,y,z$) given by the 7-by-7 spin-3 matrices

[TABLE]

[TABLE]

and

[TABLE]

The spin vector $\mathbf{F}$ can be written explicitly

[TABLE]

Define the matrices

[TABLE]

and $\mathbf{A}_{l}=(a_{l,jk})_{7\times 7}$ ($l=\pm 1,\pm 2$), where $a_{l,jk}$ is zero except for those $j+k=8-l$. For the simplicity of notations, we denote $\vec{a}_{l}=(a_{l,1(7-l)},a_{l,2(6-l)},\ldots,a_{l,(7-l)1})^{T}\in\mathcal{R}^{7-l}$ for $l=1,2$ and $\vec{a}_{l}=(a_{l,(1-l)7},a_{l,(2-l)6},\ldots,a_{l,7(1-l)})^{T}\in\mathcal{R}^{7+l}$ for $l=-1,-2$ with

[TABLE]

Then $A_{00}:=A_{00}(\Psi)=\Psi^{T}\mathbf{A}\Psi$ and $A_{2l}:=A_{2l}(\Psi)=\Psi^{T}\mathbf{A}_{l}\Psi$ can be expressed as

[TABLE]

For computing the ground state of a spin-$F$ BEC, the energy functional $E(\Phi(\cdot))$ is usually subject to the following two constraints, i.e. the mass (or normalization) as

[TABLE]

and the magnetization (with $M\in[-F,F]$) as

[TABLE]

The ground state $\Phi_{g}({\bf x})$ is obtained from the minimization of the energy functional subject to the conservation of total mass and magnetization:

Find $\left(\Phi_{g}\in S_{M}\right)$ such that

[TABLE]

where the nonconvex set $S_{M}$ is defined as

[TABLE]

For $M=\pm F$ in the spin-$F$ BEC, the constraints ensure only one component $\phi_{\pm F}$ is nonzero, and (2.21) reduces to the single component BEC ground state problems which have been considered. Therefore, we will assume $|M|<F$ for the spin-$F$ BEC ground states in the rest part of the paper.

2.1 Notoations

Given $X\in\mathcal{C}^{m\times n}$, $\bar{X},X^{T},X^{*}$ and $\Re(X)$ denote the complex conjugate, the transpose, the complex conjugate transpose, and the real part of $X$, respectively. The trace of $X$, i.e., the sum of the diagonal elements of $X\in\mathcal{C}^{n\times n}$, is denoted by $\mathrm{tr}(X)$. For a given vector $d\in\mathcal{C}^{n}$, the operator $\mathrm{diag}(d)$ returns a square matrix in $\mathcal{C}^{n\times n}$ with the elements of $d$ on the main diagonal, while $\mathrm{diag}(X)$ gives a column vector in $\mathcal{C}^{n}$ consisting of the main diagonal of $X$. The Euclidean inner product between two matrices $X,\,Y\in\mathcal{C}^{m\times n}$ is defined as $\langle X,Y\rangle:=\sum_{jk}X_{jk}\bar{Y}_{jk}=\mathrm{tr}(Y^{*}X)$.

3 Discretization Schemes

In this section, we introduce discretization of the energy functional (2.11) and constraints (2.19)- (2.20) in the constrained minimization problem (2.21) for the spin-3 case. It is similar and much easier to deal with the spin-1 and spin-2 cases. Due to the external trapping potential, the ground state of (2.21) decays exponentially as $|\mathbf{x}|\to\infty$. Thus we can truncate the energy functional and constraints from the whole space $\mathcal{R}^{d}$ to a bounded computational domain $U$ which is chosen large enough such that the truncation error is negligilbe with periodic boundary condition. Then we approximate spatial derivatives via the Fourier pseudospectral (FP) method and the integrals via the composite trapezoidal quadrature. For simplicity of notation, we only present the FP discretization in 1D. Extensions to 2D and 3D are straightforward for tensor grids and the details are omitted here for brevity.

For $d=1$, we take a bounded interval $U=(a,b)$. Let $h=(b-a)/n$ be the spatial mesh size with $n$ an even positive integer and denote $x_{j}=a+jh$ for $j=0,1,\cdots,n$. Let $\phi_{jl}$ be the numerical approximation of $\phi_{l}(x_{j})$ for $j=0,1,\cdots,n$ and $l=3,\cdots,-3$ satisfying $\phi_{0l}=\phi_{nl}$ and denote $X=(\sqrt{h}\phi_{jl})\in\mathcal{C}^{n\times 7}$ ($j=0,1,\cdots,n-1$, $l=3,\cdots,-3$).

[TABLE]

where the Fourier pseudospectral differential operator is given as

[TABLE]

with

[TABLE]

Introduce $\mathrm{V}=\mathrm{diag}(V(x_{0}),\cdots,V(x_{n-1}))$, $B=\mathrm{diag}(b)$ with $b=(b_{3},b_{2},\ldots,b_{-3})^{T}$ ($b_{l}=-pl+ql^{2}$, $l=3,\cdots,-3$), $\Lambda=\mathrm{diag}(\lambda_{-\frac{n}{2}}^{2},\cdots,\lambda_{\frac{n}{2}-1}^{2})$, and $C=(c_{jp})\in\mathcal{C}^{n\times n}$ with entries $c_{jp}=e^{-i\frac{2\pi jp}{n}}$ for $j=0,\cdots,n-1$ and $p=-\frac{n}{2},\cdots,\frac{n}{2}-1$. Plugging (3.25) and (3.26) into (3.24), and replacing $\phi_{l}(x_{j})$ by $\phi_{jl}$, we get the finite dimensional approximation to the energy functional defined as

[TABLE]

where $L=C^{*}\Lambda C$ is the matrix representation of the discrete negative Laplace operator and

[TABLE]

are column vectors. In fact, the first term in (3.27) can be computed efficiently at cost $O(n\ln n)$ through the discretized Fourier transform (DFT).

Similarly, let $D=\mathrm{diag}(3,\cdots,-3)$, the constraints (2.19)- (2.20) can be truncated and discretized as

[TABLE]

which immediately implies that the set $S_{M}$ can be discretized as

[TABLE]

Hence, the original problem (2.21) with $d=1$ can be approximated by the discretized minimization problem via the FP discretization:

[TABLE]

To solve the discrete minimization problem (3.33), it is often necessary to compute the gradient and Hessian matrix of the discrete energy $E_{h}(X)$. The second-order Taylor expansion of $E_{h}(X)$ can be expressed as

[TABLE]

where h.o.t. is short for the higher-order terms. By a simple calculation, we can get the gradient

[TABLE]

and the Hessian-vector product

[TABLE]

4 Manifold Structure

In the ground state of a spin-$F$ BEC, we have $M\leftrightarrow-M\Longleftrightarrow\phi_{l}\leftrightarrow\phi_{-l}.$ Thus we only discuss the cases where $M\geq 0$. Express $X$ as $X=X_{r}+iX_{i}$, where $X_{r},X_{i}\in\mathcal{R}^{n\times(2F+1)}$. Let $(X_{r};X_{i})=(u_{F},u_{F-1},\cdots,u_{0},\cdots,u_{-F+1},u_{-F})\in\mathcal{R}^{2n\times(2F+1)}$ and $u\in\mathcal{R}^{N}$ be the reconstructed column vector of this matrix, where $N=2n(2F+1)$. Introduce

[TABLE]

where $I_{2n}$ denotes the identity matrix of size $2n$. Then the constraints can be discretized as

[TABLE]

in which $0\leq M<F.$ Define $\widetilde{E}(u):=E_{h}(X)$, our model problem can be formulated as

[TABLE]

This is a nonconvex optimization problem with constraints. Observe that $\mathcal{M}$ is a level set of the function

[TABLE]

When $M\notin\mathcal{Z}$, $\nabla\mathbf{G}(u)=(u,\Gamma u)^{T}$ has full rank at every point $u\in\mathcal{M}$, thus according to Proposition 3.3.3 in [2], $\mathcal{M}$ is a closed embedded submanifold of $\mathcal{R}^{N}$ of dimension $N-2$. In the following discussion, we will let $u$ be an arbitrary point on the manifold $\mathcal{M}$ and assume $f$ is a smooth real-valued function in a neighborhood of $u$.

Given a curve $\gamma(t):\mathcal{R}\rightarrow\mathcal{M}\subset\mathcal{R}^{N}$ through $u$ at $t=0$, the associated tangent vector $\dot{\gamma}(0)$ can be represented by $\gamma^{\prime}(0)$ in the way that

[TABLE]

Since $\mathcal{M}$ is a level set of the constant-rank function G, the tangent space $T_{u}\mathcal{M}$ reads

[TABLE]

We naturally define the inner product $\langle\cdot,\cdot\rangle_{u}$ and the norm $\|\cdot\|_{u}$ on $T_{u}\mathcal{M}$ as

[TABLE]

Under such a metric, the Riemannian gradient grad$f(u)$, defined as the unique element of $T_{u}\mathcal{M}$ satisfying

[TABLE]

can be written as

[TABLE]

where $\mathrm{P}_{u}$ denotes the orthogonal projection from $\mathcal{R}^{N}$ onto $T_{u}\mathcal{M}$. From (4.42) we can easily derive the formula of $\mathrm{P}_{u}$:

Lemma 4.1** ($\mathrm{P}_{u}$).**

For an arbitrary point $w\in\mathcal{R}^{N}$, the orthogonal projection of it onto $T_{u}\mathcal{M}$ reads

[TABLE]

Proof 4.2**.**

From (4.42) and the definition of $\mathrm{P}_{u}$, we have

[TABLE]

and $w-\mathrm{P}_{u}w\in(T_{u}\mathcal{M})^{\perp}$, therefore there exists $\alpha_{w},\beta_{w}\in\mathcal{R}$ such that

[TABLE]

Noticing that $\mathrm{P}_{u}w\in T_{u}\mathcal{M}$ which implies

[TABLE]

We can obtain from (4.48) that

[TABLE]

In view of the fact that $u^{T}u=1$ and $u^{T}\Gamma u=M$, (4.50) can be simplified as

[TABLE]

It follows directly from Cauchy-Schwarz inequality that

[TABLE]

which ensures the linear system (4.51) has a unique solution:

[TABLE]

*Substituting (4.53) into (4.48) yields the formula (4.46). *

Let $\mathfrak{X}(\mathcal{M})$ be the set of smooth vector fields on $\mathcal{M}$. The Riemannian Hessian Hess$f(u)$ is a linear mapping from $T_{u}\mathcal{M}$ into itself defined as

[TABLE]

where $\widetilde{\nabla}$ denotes the Riemannian connection of $\mathcal{M}$. Since $\mathcal{M}$ is a Riemannian submanifold of $\mathcal{R}^{N}$, according to [2] its Riemannian connection reads

[TABLE]

Thus we have

[TABLE]

The formula of Hess$f(u)$ is given in following lemma:

Lemma 4.3** ($\mathrm{Hess}f(u)$).**

Given a tangent vector $\xi\in T_{u}\mathcal{M}$, and let $g$ and $H$ be the Euclidean gradient and Euclidean Hessian of $f$ respectively, then

[TABLE]

where

[TABLE]

*and $h_{e}:=H(u)\cdot\xi,\alpha_{u}:=u^{T}\Gamma^{2}u-M^{2},\beta_{u}:=u^{T}\Gamma h_{e}-\beta_{g}u^{T}\Gamma^{2}\xi$. *

Proof 4.4**.**

Recalling Lemma 4.1 and (4.45), we get

[TABLE]

and

[TABLE]

Since $(\nabla\alpha_{g}^{T}\xi)u,~{}(\nabla\beta_{g}^{T}\xi)\Gamma u\in(T_{u}\mathcal{M})^{\perp}$, we have

[TABLE]

*For $\xi\in T_{u}\mathcal{M}$, Lemma 4.1 and (4.60) lead to the formula (4.57). *

The first-order and second-order optimality conditions for optimization problems on Riemannian manifolds coincide with the conventional ones [33]. If $u^{\ast}$ is a local solution of (4.39), we have grad$\widetilde{E}(u^{\ast})=0$ and all the points $u$ at which grad$\widetilde{E}(u)=0$ are called stationary points of $\widetilde{E}$.

Line search optimization methods in $\mathcal{R}^{N}$ are based on the update formula

[TABLE]

where $\eta_{k}\in\mathcal{R}^{N}$ is the search direction and $t_{k}>0$ is the step size. Correspondingly, when (4.61) is generalized to a manifold, $\eta_{k}$ is selected as a tangent vector, and the line search procedure relies on the concept of retraction:

Definition 4.5** (retraction).**

A retraction on a manifold $\mathcal{M}$ is a smooth mapping $R$ from the tangent bundle $T\mathcal{M}$ onto $\mathcal{M}$ with the following properties. Let $R_{u}$ denote the restriction of $R$ to $T_{u}\mathcal{M}$.

• (i)

$R_{u}(0_{u})=u$, where $0_{u}$ denotes the zero element of $T_{u}\mathcal{M}$.

• (ii)

With the canonical identification $T_{0_{u}}T_{u}\mathcal{M}\simeq T_{u}\mathcal{M}$, $R_{u}$ satisfies

[TABLE]

where $\mathrm{id}_{T_{u}\mathcal{M}}$ denotes the identity mapping on $T_{u}\mathcal{M}$.

Remark 4.6**.**

When $M\in\mathcal{Z}$, $\mathcal{M}$ is not a well-defined manifold. However, by restricting the feasible region to

[TABLE]

*we can also define above structures and the formulas still work. This modification does not change our numerical experiments. *

5 A Modified Adaptive Regularized Newton Method

We aim to solve (4.39) with a modified version of the adaptive regularized Newton method (ARNT) developed in [21]. At the $k$-th iteration, ARNT uses a second-order Taylor model with a penalization term to approximate the original objective function but keeps the constraint $u\in\mathcal{M}$.

Specifically, the method replaces (4.39) with a sequence of quadratic subproblems:

[TABLE]

where $\langle\cdot,\cdot\rangle$ denotes the dot product in $\mathcal{R}^{N}$ and $H_{k}$ is the Euclidean Hessian of $\widetilde{E}$ at $u_{k}$. The subproblem (5.63) is solved approximately by applying a modified conjugate gradient (CG) method to the linear system

[TABLE]

The method terminates when either certain accuracy is reached or negative curvature is detected. It outputs two vectors $s_{k}$ and $d_{k}$, where $s_{k}$ is the solution computed by CG method and $d_{k}$ represents the negative curvature information. The new search direction $\xi_{k}$ is chosen as

[TABLE]

which is a descent direction (cf. Lemma 7, [21]).

After construction of $\xi_{k}$, a monotone Armijo-based curvilinear search is conducted to generate a trial point

[TABLE]

where $\varsigma$ is the smallest integer satisfying

[TABLE]

and $\rho,\delta\in(0,1),\alpha_{0}\in(0,1]$ are given constants.

In order to monitor the acceptance of the trial point $z_{k}$ and adjust the regularization parameter $\sigma_{k}$, the above procedure is embedded in a trust-region framework where $\sigma_{k}$ plays a similar role as the trust-region radius and is updated according to the ratio

[TABLE]

ARNT may exhibit a certain instability when directly applied to solve (4.39). To improve its performance, we combine it with the cascadic multigrid method in [13]. In detail, we first solve (4.39) on the coarsest mesh, and then use the obtained solution as the initial guess of the problem on a finer mesh, and repeat until reaching the finest mesh.

On each mesh, we use the Riemannian gradient method with a BB step size (RGBB) in [21] to compute an initial point for ARNT. At the $k$-th iteration, RGBB performs a nonmonotone Armijo-based curvilinear search along the steepest descent direction $\eta_{k}=-\mathrm{grad}\,\widetilde{E}(u_{k})$. Given $\rho,\varrho,\delta\in(0,1)$, it tries to find the smallest integer $\varsigma$ satisfying

[TABLE]

where the initial step size $\gamma_{k}$ is computed as $\gamma_{k}=|\langle s_{k-1},v_{k-1}\rangle_{u_{k}}|/\langle v_{k-1},v_{k-1}\rangle_{u_{k}}$, with $s_{k-1}=u_{k}-u_{k-1},\quad v_{k-1}=\mathrm{grad}\widetilde{E}(u_{k})-\mathrm{grad}\widetilde{E}(u_{k-1}).$ The value $C_{k+1}$ is calculated via $C_{k+1}=(\varrho Q_{k}C_{k}+\widetilde{E}(u_{k+1}))/Q_{k+1}$, with $C_{0}=\widetilde{E}(u_{0}),Q_{k+1}=\varrho Q_{k}+1$ and $Q_{0}=1$.

The modified adaptive regularized Newton method (still referred to as ARNT in this paper) is presented in Algorithm 1.

6 Retractions

The selection of retractions can affect the performance of Riemannian optimization algorithms. In this section, we try to find retractions of the form

[TABLE]

where $\psi$ is some “projection” from a neighborhood of $\mathcal{M}$ in $\mathcal{R}^{N}$ to $\mathcal{M}$.

For $w=(w_{F},w_{F-1},\ldots,w_{-F})\in\mathcal{R}^{N}$ ($w_{l}\in\mathcal{R}^{2n}$, $N=2n(2F+1)$, $l=F,\ldots,-F$), we define two functions $f_{1}(w),f_{2}(w)\in\{F,F-1,\ldots,-F\}$:

[TABLE]

Observe that at every point $u\in\mathcal{M}$, the constraints indicate

[TABLE]

Thus when $M\notin\mathcal{Z}$, we have $f_{1}(u)<M,f_{2}(u)>M$, which is equivalent to

[TABLE]

when $M\in\mathcal{Z}$, (6.73) also holds for $u\in\widetilde{\mathcal{M}}$. Define the open set $\Omega$ as

[TABLE]

then $\Omega\supset\mathcal{M}$ for $M\notin\mathcal{Z}$ and $\Omega\supset\widetilde{\mathcal{M}}$ for $M\in\mathcal{Z}$. For the simplicity of presentation, we will discuss three different retractions from $\Omega$ to $\mathcal{M}$ ($M\notin\mathcal{Z}$) and all the results hold also for $\widetilde{\mathcal{M}}$ ($M\in Z$).

6.1 Projective Retraction

The most intuitive retraction $\psi$ is given by the projection operator $\mathcal{P}_{\mathcal{M}}$, which is defined as

[TABLE]

According to [3], $\mathcal{P}_{\mathcal{M}}$ is a well-defined function (existence and uniqueness of the projection hold) in a neighborhood $\widetilde{\Omega}\subset\Omega$ of $\mathcal{M}$, and the mapping $R_{u}(\xi_{u}):=\mathcal{P}_{\mathcal{M}}(u+\xi_{u})$ is a well-defined retraction on $\mathcal{M}$, called the projective retraction in this paper. The explicit formula is given as follows.

Lemma 6.1**.**

For an arbitrary point $w\in\widetilde{\Omega}$, $\mathcal{P}_{\mathcal{M}}(w)$ reads

[TABLE]

where

[TABLE]

and $r$ is the unique zero of the function

[TABLE]

Proof 6.2**.**

Define the Lagrangian function of (6.75) as

[TABLE]

Let $z=\mathcal{P}_{\mathcal{M}}(w)$, then the condition $\nabla_{l}L(z,\mu,\lambda)=0$ gives

[TABLE]

Since $\|z-w_{l}\|_{2}^{2}=\|z\|_{2}^{2}+\|w\|_{2}^{2}-2\sum_{l}w_{l}^{T}z_{l}=1+\|w\|_{2}^{2}-2\sum_{l}w_{l}^{T}z_{l}$ ($z\in\mathcal{M}$), the projection $z$ maximizes $w_{l}^{T}z_{l}$ which leads to $1-\mu-l\lambda>0$ for $w_{l}\neq 0$. On the other hand, if $w_{l}=0$ then $z_{l}=0$, otherwise substituting $z_{l}$ with$-z_{l}$ yields a different projection of $w$, which contradicts the uniqueness.

Since $w_{f_{1}(w)}$ and $w_{f_{2}(w)}$ are nonzero, we have

[TABLE]

which is equivalent to $1-\mu-M\lambda>(f_{1}(w)-M)\lambda$ , $1-\mu-M\lambda>(f_{2}(w)-M)\lambda$ , and

[TABLE]

The inequalities in (6.81) indicate that $1-\mu-l\lambda>0$ for $l=f_{1}(w)+1,...,f_{2}(w)-1,$ and from (6.80) and $z\in\mathcal{M}$ we have

[TABLE]

In view of (6.82) - (6.83), denoting

[TABLE]

recalling the definition of function $h_{1}(\cdot)$ in (6.78), we have $s>0,\quad\frac{1}{f_{1}(w)-M}<r<\frac{1}{f_{2}(w)-M}$ and

[TABLE]

For any $t\in(\frac{1}{f_{1}(w)-M},\frac{1}{f_{2}(w)-M})$, we have

[TABLE]

In addition, noticing that

[TABLE]

$h_{1}(\cdot)$* has exactly one zero in $(\frac{1}{f_{1}(w)-M},\frac{1}{f_{2}(w)-M})$. Substituting (6.85) into (6.84), the formulas of $\lambda$ and $\mu$ can be obtained accordingly. *

We remark that (6.1) can be applied to any $w\in\Omega$. For spin-1 case, the closed-form solution of (6.75) is computable.

Lemma 6.3**.**

When $F=1$, given any nonzero $w\in\Omega$, the optimal solution $z=(z_{1};z_{0};z_{-1})$ of (6.75) is

(1) If $M=0$, then $z_{0}=w_{0}/t$ and

[TABLE]

with $t=\sqrt{\|w_{0}\|_{2}^{2}+(\|w_{1}\|_{2}+\|w_{-1}\|_{2})^{2}/2}$.

(2) If $M>0$ and $w_{0}=0$, then $z_{0}=0$,

[TABLE]

(3) If $M>0$ and $w_{-1}=0$, then $z_{-1}=0$,

[TABLE]

(4) If $M>0$ and $\|w_{0}\|,\|w_{-1}\|>0$, then $z_{l}=w_{l}/(1-\mu-l\lambda)$ ($l=\pm 1,0$) with

[TABLE]

where

[TABLE]

with $\beta$ depending on $\|w_{l}\|$ as

[TABLE]

and

[TABLE]

[TABLE]

Proof 6.4**.**

The first three cases are straightforward to verify. Here we only present the proof for case (4). If $M>0$ and $\|w_{0}\|_{2},\|w_{-1}\|_{2}>0$, then $z_{l}=w_{l}/(1-\mu-l\lambda)$ ($l=0,\pm 1$). Let $\|z_{0}\|_{2}=s$, we have from (6.80) that $1-\mu=\|w_{0}\|_{2}/s$ and

[TABLE]

which implies

[TABLE]

*There exists a unique solution $s\in(0,\sqrt{1-M})$ and the Lagrange multipliers can be identified. *

6.2 Orthogonal Retraction

Inspired by the projective retraction, we can consider $\psi$ of the form

[TABLE]

with undetermined positive coefficients $\sigma_{F},...,\sigma_{-F}$. Besides the constraints

[TABLE]

we have to introduce additional $2F-1$ conditions to uniquely determine the $2F+1$ coefficients.

In [9], Bao and Lim proposed the condition $\sigma_{1}\sigma_{-1}=\sigma_{0}^{2}$ for spin-1 BEC. It can be generalized to

[TABLE]

for spin-$F$ cases. The mapping $R$ characterized by (6.70) and (6.99) - (6.101) is called the orthogonal retraction in this paper.

Lemma 6.5**.**

For an arbitrary point $w\in\Omega$, $\psi(w)$ defined by (6.99) - (6.101) reads

[TABLE]

where $r$ is the unique positive zero of the polynomial

[TABLE]

Proof 6.6**.**

Let $r=\sigma_{1}^{2}/\sigma_{0}^{2}>0$, then

[TABLE]

Substituting (6.104) into (6.100) yields

[TABLE]

Introducing $h_{3}(t)=h_{2}(t)t^{-f_{1}(w)-F}$, we have

[TABLE]

and

[TABLE]

which implies the function $h_{3}(\cdot)$ has at least one positive zero.

Let $\displaystyle\sum_{1}$ and $\displaystyle\sum_{2}$ denote $\displaystyle\sum_{1\leq l<M-f_{1}(w)}$ and $\displaystyle\sum_{M-f_{1}(w)<l\leq f_{2}(w)-f_{1}(w)}$ respectively. At a zero $r_{0}$ of $h_{3}$, we have

[TABLE]

which leads to

[TABLE]

*From (6.108) we can see that $h_{3}$ has exactly one positive zero, and $h_{2}$ has exactly one positive zero too. Substituting (6.105) into (6.104) leads to the formulas of the coefficients. *

Noticing that in spin-1 cases, $h_{2}$ degenerates to a quadratic polynomial, and the orthogonal retraction has a closed form solution.

The well-definedness of the orthogonal retraction is guaranteed by following theorem [2]:

Theorem 6.7**.**

Let $\mathcal{M}$ be an embedded manifold of a vector space $\mathcal{E}$ and let $\mathcal{N}$ be an abstract manifold such that $\dim(\mathcal{M})+\dim(\mathcal{N})=\dim(\mathcal{E})$. Assume that there is a diffeomorphism

[TABLE]

where $\mathcal{E}_{\ast}$ is an open subset of $\mathcal{E}$, with a neural element $e$ satisfying

[TABLE]

Then the mapping

[TABLE]

*where $\pi_{1}:\mathcal{M}\times\mathcal{N}\rightarrow\mathcal{M}$ is the projection onto the first component, defines a retraction on $\mathcal{M}$. *

Lemma 6.8**.**

*The orthogonal retraction is a well-defined retraction on $\mathcal{M}$. *

Proof 6.9**.**

Take $\mathcal{N}=\mathcal{R}_{+}^{2}$, and $\mathcal{N}$ is a manifold satisfying $\dim(\mathcal{M})+\dim(\mathcal{N})=\dim(\mathcal{R}^{N})$. Define the mapping $\varphi:\mathcal{M}\times\mathcal{N}\rightarrow\Omega$ as

[TABLE]

Lemma 6.5 shows for any $w\in\Omega$ there exists a unique $u=\psi(w),v=(\sigma_{0},r)^{T}$ such that $\varphi(u,v)=w$, thus $\varphi$ is a bijection. It is obvious to see that $\varphi$ is smooth on $\mathcal{M}\times\mathcal{N}$, and $\varphi(u,\mathbf{1})=u,~{}\forall u\in\mathcal{M}$.

From Lemma 6.5, we have

[TABLE]

where $\sigma_{0}=\sqrt{\frac{1}{\sum_{l=-F}^{F}\|w_{l}\|_{2}^{2}r^{l}}}$ and $r$ is characterized by the equation

[TABLE]

Since $h_{2}^{\prime}(r)=(h_{3}(t)\cdot t^{f_{1}(w)+F})^{\prime}|_{t=r}=h_{3}^{\prime}(r)\cdot r^{f_{1}(w)+F}>0$, it follows from the implicit function theorem that $r$, when considered as a function of $w$, is smooth. Then $\varphi^{-1}$ is also a smooth function at every $w\in\Omega$, which makes $\varphi$ a diffeomorphism. Thus the orthogonal retraction, given by

[TABLE]

*is a retraction on $\mathcal{M}$. *

6.3 Closed-form Retraction

In the projective retraction, the coefficients take the form

[TABLE]

When $w\rightarrow\mathcal{M}$, we have $\sigma_{l}\rightarrow 1$ and $\mu,\lambda\rightarrow 0$, and

[TABLE]

Thus we can approximate the projective retraction by taking

[TABLE]

As shown below, (6.115) has a closed-form formula, and the mapping $R$ characterized by (6.70), (6.99), (6.100) and (6.115) is called the closed-form retraction in this paper. Firstly, we introduce

[TABLE]

Apparently, $S$ is an open set and $\mathcal{M}\subset S\subset\Omega$. We next discuss the computation of $\psi(w)$ for $w\in S$.

Lemma 6.10**.**

For an arbitrary point $w\in S$, $\psi(w)$ defined by (6.99), (6.100) and (6.115) reads

[TABLE]

Proof 6.11**.**

Substituting $\eqref{eq:retr3}$ into $\eqref{eq:oretr2}$ yields

[TABLE]

and the solution is given by

[TABLE]

*The condition $w\in S$ ensures $1+2\mu+2l\lambda>0$ for $l=F,...,-F$. In view of the retraction (6.115), we obtain the formula (6.117). *

Lemma 6.12**.**

*The closed-form retraction is a well-defined retraction. *

Proof 6.13**.**

Denote

[TABLE]

$\mathcal{N}$* is an open subset of $\mathcal{R}^{2}$ and therefore a 2-dimensional manifold. Define the mapping $\varphi:\mathcal{M}\times\mathcal{N}\rightarrow\mathcal{R}^{N}$ as*

[TABLE]

For an arbitrary point $(u,v)\in\mathcal{M}\times\mathcal{N}$, let $w=\varphi(u,v)$.

• •

If $v_{2}=0$, then $w=\frac{1}{\sqrt{1+2v_{1}}}u$ ($u\in\mathcal{M}$), and

[TABLE]

• •

If $v_{2}\neq 0$, then

[TABLE]

Noticing that

[TABLE]

we have

[TABLE]

On one hand, above analysis shows $\varphi(\mathcal{M}\times\mathcal{N})\subset S$; on the other hand, Lemma 6.10 indicates that for any $w\in S$, there exists a unique $u=\psi(w),v=(\mu,\lambda)^{T}$ such that $\varphi(u,v)=w$. Hence $\varphi$ is a bijection from $\mathcal{M}\times\mathcal{N}$ to $S$. It is straightforward to see that $\varphi$ and $\varphi^{-1}$ are both smooth functions, and $\varphi(u,\mathbf{0})=u,\forall u\in\mathcal{M}$. Thus from Theorem 6.7 we know that the closed-form retraction, given by

[TABLE]

*is a retraction on $\mathcal{M}$. *

7 Numerical Experiments

In this section, we first compare the performance of Algorithm 1 with RGBB and the Riemannian trust region method (RTR) [1] by testing some BEC examples. RGBB and RTR are also expedited with the mesh refinement technique. We present numerical results of these algorithms under the three different retractions defined in Section 6. Then we apply Algorithm 1 to compute the ground states of spin-2 and spin-3 BEC with different parameters. All codes are written in MATLAB. All experiments were performed on a workstation with Intel Xenon E5-2680 v3 processors at 2.50GHz($\times 12$) and 128GB memory running CentOS 6.8 and MATLAB R2018b.

In the spin-1 BEC, the initial data is chosen as $\Phi_{0}({\bf x})=U\phi_{0}({\bf x})$, where

[TABLE]

$U=\left(\frac{\sqrt{1+3M}}{2},\sqrt{\frac{1-M}{2}},\frac{\sqrt{1-M}}{2}\right)^{T}$ for the ferromagnetic interaction ($\beta_{1}\leq 0$); $U=\left(\sqrt{\frac{1+M}{2}},0,\sqrt{\frac{1-M}{2}}\right)^{T}$ for the antiferromagnetic interaction ($\beta_{1}>0$) [9].

In the spin-2 BEC, the initial data is chosen as $\Phi_{0}({\bf x})=U\phi_{0}({\bf x})$, where

[TABLE]

with $m_{1}=\sqrt{2+M}$ and $m_{2}=\sqrt{2-M}$ for the ferromagnetic interaction ($\beta_{1}<0$ and $\beta_{2}>20\beta_{1}$). And $U=\left(\frac{\sqrt{2+M}}{2},0,0,0,\frac{\sqrt{2-M}}{2}\right)^{T}$ for the nematic interaction ($\beta_{2}<0$ and $\beta_{2}<20\beta_{1}$), $U=\left(\sqrt{\frac{M+1}{3}},0,0,\sqrt{\frac{2-M}{3}},0\right)^{T}$ for the cyclic interaction ($\beta_{1}>0$ and $\beta_{2}>0$) [10]. In the spin-3 BEC, the initial data is chosen as $\Phi_{0}({\bf x})=U\phi_{0}({\bf x})$, where $U\in\mathcal{R}^{7}$ is taken as the random vector. In all the examples, we take $p=q=0$.

7.1 Performance of algorithms

In RGBB we used all of the default parameters. As for RTR, we added a rule $\|r_{j+1}\|_{2}\leq\min\{0.1,0.1\|r_{0}\|_{2}\}$ into the stopping criterion of the truncated CG method. All other default settings of RTR were used. For ARNT, we set $\eta_{1}=0.01,\eta_{2}=0.9,\gamma_{0}=0.2,\gamma_{1}=1,\gamma_{2}=10$, and $\sigma_{k}=\hat{\sigma}_{k}\|\mathrm{grad}\tilde{E}(u_{k})\|_{2}$, where $\hat{\sigma}_{k}$ is updated by the procedure in Algorithm 1 with $\hat{\sigma}_{0}=1$. Furthermore, when an estimation of the absolute value of the negative curvature, denoted by $\sigma_{est}$, is available at the $k$-th subproblem, we can calculate

[TABLE]

with some small $\tilde{\gamma}\geq 0$. Then the parameter $\sigma_{k+1}$ is reset to $\sigma_{k+1}^{new}$.

On the finest mesh, all algorithms terminate when either $\|\mathrm{grad}\tilde{E}(u_{k})\|_{2}\leq 10^{-6}$ or the number of iterations reaches 10000, while on the coarse meshes they all terminate when $\|\mathrm{grad}E(x_{k})\|_{2}\leq 10^{-5}$. In the implementation of ARNT, RGBB stops when either $\|\mathrm{grad}\tilde{E}(u_{k})\|_{2}\leq 10^{-2}$ or the number of iterations reaches 2000. The maximum number of inner iterations in ARNT is chosen adaptively depending on $\|\mathrm{grad}\tilde{E}(u_{k})\|_{2}$.

In the subsequent tables, the columns “f”, “nrmG” and “time” display the final objective function value, the final norm of the Riemannian gradient and the total CPU time each algorithm spent to reach the stopping criterion. The column “iter” reports the number of iterations (the average numbers of inner iterations) on the finest mesh. The choice of retractions is shown in the column “retr”, where R1, R2 and R3 denote the projective retraction, the orthogonal retraction and the closed-form retraction, respectively.

We present results of following cases for spin-1, spin-2 and spin-3 BEC:

• •

Spin-1 BEC [9]

• 2D. Antiferromagnetic case. $V(x,y)=\frac{1}{2}(x^{2}+y^{2})+10[\sin^{2}(\frac{\pi x}{4})+\sin^{2}(\frac{\pi y}{4})]$, $\beta_{0}=300$, $\beta_{1}=100$, $U=[-16,16]\times[-16,16],n=2^{9}.$
• 3D. Ferromagnetic case. $V(x,y,z)=\frac{1}{2}\sum_{\alpha=x,y,z}\left(\alpha^{2}+200\sin^{2}(\frac{\pi\alpha}{2})\right)$, $\beta_{0}=880$, $\beta_{1}=-4.1$, $U=[-16,16]\times[-16,16]\times[-16,16],n=2^{8}.$
• •

Spin-2 BEC [18]

• 2D. Antiferromagnetic case. $V(x,y)=\frac{1}{2}(x^{2}+y^{2})+10\left[\sin^{2}(\frac{\pi x}{2})+\sin^{2}(\frac{\pi y}{2})\right]$, $\beta_{0}=243$, $\beta_{1}=12.1$, $\beta_{2}=-13$, $U=[-8,8]\times[-8,8],n=2^{8}.$
• 3D. Cyclic case. $V(x,y,z)=\frac{1}{2}(x^{2}+y^{2}+z^{2})+100\left[\sin^{2}(\frac{\pi x}{2})+\sin^{2}(\frac{\pi y}{2})+\sin^{2}(\frac{\pi z}{2})\right]$, $\beta_{0}=183.9$, $\beta_{1}=26.8$, $\beta_{2}=134.7$, $U=[-16,16]\times[-16,16]\times[-16,16],n=2^{8}.$
• •

Spin-3 BEC

• 2D. $V(x,y)=\frac{1}{2}(x^{2}+y^{2})+10\left[\sin^{2}(\frac{\pi x}{2})+\sin^{2}(\frac{\pi y}{2})\right]$, $\beta_{0}=100$, $\beta_{1}=1,\beta_{2}=10,\beta_{3}=1$, $U=[-8,8]\times[-8,8],n=2^{8}.$
• 3D. $V(x,y,z)=\frac{1}{2}(x^{2}+y^{2}+z^{2})+100\left[\sin^{2}(\frac{\pi x}{2})+\sin^{2}(\frac{\pi y}{2})+\sin^{2}(\frac{\pi z}{2})\right]$, $\beta_{0}=100$, $\beta_{1}=1,\beta_{2}=10,\beta_{3}=1$, $U=[-8,8]\times[-8,8]\times[-8,8],n=2^{7}.$

The detailed numerical results are reported in Tables 1-6. In most cases, all three algorithms converge to points with the same function values. For spin-1 and spin-2 cases, the choice of different retractions has small impact on the numerical performance, and the second-order algorithms ARNT and RTR exhibit higher stability than the first-order algorithm RGBB in response to the change of retractions. In the 3D case of spin-3 BEC (Table 6), RTR converges to a larger function value than ARNT when $M=1.5$ using retractions R2 and R3. Overall, taking both numerical stability and time cost into consideration, ARNT shows the best performance.

7.2 Application in spin-2 BEC

In this section, we apply the ARNT method with the projective retraction to compute the ground state of a spin-2 BEC in 1-3 dimensions and under different interactions. Specifically, the following cases are considered [18]:

• •

1D, $V(x)=\frac{1}{2}x^{2}+25\sin^{2}(\frac{\pi x}{4}).$

• Case I (ferromagnetic).   $\beta_{0}=130.6$, $\beta_{1}=-25.4,\beta_{2}=-125.3$, $U=[-8,8],n=2^{8}.$
• Case II (antiferromagnetic).  $\beta_{0}=243$, $\beta_{1}=12.1,\beta_{2}=-13$, $U=[-16,16],n=2^{9}.$
• Case III (cyclic). $\beta_{0}=183.9$, $\beta_{1}=26.8,\beta_{2}=134.7$, $U=[-16,16],n=2^{9}.$
• •

2D, $V(x,y)=\frac{1}{2}(x^{2}+y^{2})+10\left[\sin^{2}(\frac{\pi x}{2})+\sin^{2}(\frac{\pi y}{2})\right]$, $U=[-8,8]\times[-8,8],n=2^{8}.$

• Case I (ferromagnetic).   $\beta_{0}=130.6$, $\beta_{1}=-25.4,\beta_{2}=-125.3$.
• Case II (antiferromagnetic).  $\beta_{0}=243$, $\beta_{1}=12.1,\beta_{2}=-13$.
• Case III (cyclic). $\beta_{0}=183.9$, $\beta_{1}=26.8,\beta_{2}=134.7$.
• •

3D, $V(x,y,z)=\frac{1}{2}(x^{2}+y^{2}+z^{2})+100\left[\sin^{2}(\frac{\pi x}{2})+\sin^{2}(\frac{\pi y}{2})+\sin^{2}(\frac{\pi z}{2})\right]$

• Case I (ferromagnetic).   $\beta_{0}=130.6$, $\beta_{1}=-25.4,\beta_{2}=-125.3$, $U=[-8,8]\times[-16,16]\times[-16,16],n=2^{7}.$
• Case II (antiferromagnetic).  $\beta_{0}=243$, $\beta_{1}=12.1,\beta_{2}=-13$, $U=[-16,16]\times[-16,16]\times[-16,16],n=2^{8}.$
• Case III (cyclic). $\beta_{0}=183.9$, $\beta_{1}=26.8,\beta_{2}=134.7$, $U=[-16,16]\times[-16,16]\times[-16,16],n=2^{8}.$

The ground state energies in above cases are listed in Table 7. Under ferromagnetic interaction (Case I) the energy remains constant when $M$ changes; under antiferromagnetic interaction (Case II), the energy gets higher as $M$ increases; under cyclic interaction (Case III), the energy decreases first and then goes up as $M$ increases.

The wave functions of the ground states computed by ARNT are given in Figure 1- 3. The peak function value under ferromagnetic interaction is lower than that under the other two types of interactions. By comparing Figures 1, 2 and 3, we can find a common property: when $M>0$, in the ground states under nematic interaction, the components $\phi_{1},\phi_{0},\phi_{-1}$ are always zero-valued functions; and in the ground states under cyclic interaction, the components $\phi_{1},\phi_{0},\phi_{-2}$ are always zero-valued functions.

7.3 Application in spin-3 BEC

In this section, we apply the ARNT method with the projective retraction to compute the ground state of a spin-3 BEC in 1-3 dimensions under different interactions. In detail, the following cases are considered:

• •

1D, $V(x)=\frac{1}{2}x^{2}+25\sin^{2}(\frac{\pi x}{4})$, $U=[-8,8],n=2^{8}.$

• Case I.   $\beta_{0}=100$, $\beta_{1}=1,\beta_{2}=-10,\beta_{3}=-1$.
• Case II.   $\beta_{0}=100$, $\beta_{1}=1,\beta_{2}=10,\beta_{3}=1$.
• •

2D, $V(x,y)=\frac{1}{2}(x^{2}+y^{2})+10\left[\sin^{2}(\frac{\pi x}{2})+\sin^{2}(\frac{\pi y}{2})\right]$, $U=[-8,8]\times[-8,8],n=2^{8}.$

• Case I.   $\beta_{0}=100$, $\beta_{1}=1,\beta_{2}=-10,\beta_{3}=-1$.
• Case II.   $\beta_{0}=100$, $\beta_{1}=1,\beta_{2}=10,\beta_{3}=1$.
• •

3D, $V(x,y,z)=\frac{1}{2}(x^{2}+y^{2}+z^{2})+100\left[\sin^{2}(\frac{\pi x}{2})+\sin^{2}(\frac{\pi y}{2})+\sin^{2}(\frac{\pi z}{2})\right]$, $U=[-8,8]\times[-8,8]\times[-8,8],n=2^{7}.$

• Case I.   $\beta_{0}=100$, $\beta_{1}=1,\beta_{2}=-10,\beta_{3}=-1$.
• Case II.   $\beta_{0}=100$, $\beta_{1}=1,\beta_{2}=10,\beta_{3}=1$.

The ground state energies in above cases are listed in Table 8. In each case, the energy increases as $M$ increases. The wave functions of the ground states computed by ARNT are given in Figures 4-6. By comparing the figures, we can see that in Case I, when $M>0$, the components $\phi_{2},\phi_{1},\phi_{0},\phi_{-1},\phi_{-2}$ are always close to zero; in Case II, the components $\phi_{2},\phi_{1},\phi_{-1},\phi_{-2}$ are always close to zero ($\infty$-norm less than $10^{-6}$).

8 Conclusions

The Fourier pseudospectral method was adopted to discretize the energy functional and constraints for computing the ground states of spin-$F$ Bose-Einstein condensate (BEC). The original variational problem was reduced to a finite dimensional Riemannian optimization problem. An adaptive regularized Newton method, combined with a Riemannian gradient method and a cascadic multigrid technique, was designed to solve the discretized problem. Three different retractions were proposed to implement the optimization algorithms on the manifold. Comparison with the Riemannian gradient method and trust-region method for different retractions and parameters showed that our approach is more efficient and stable. Extensive numerical examples of spin-2 and spin-3 BEC in 1D, 2D and 3D with optical lattice potential and various interaction demonstrated the robustness of our approach. The energy and wave functions of ground states are reported to reveal some interesting physical phenomena. Our method is the first one to explore spin-3 BEC computationally. Although the spin-3 cases discussed in this paper are relatively simple, our algorithm is also applicable for cases with more diverse parameters.