Ground states, solitons and spin textures in spin-1 Bose–Einstein condensates

Abstract

We present an overview of our recent theoretical studies on the quantum phenomena of the spin-1 Bose–Einstein condensates, including the phase diagram, soliton solutions and the formation of the topological spin textures. A brief exploration of the effects of spin–orbit coupling on the ground-state properties is given. We put forward proposals by using the transmission spectra of an optical cavity to probe the quantum ground states: the ferromagnetic and polar phases. Quasi-one-dimension solitons and ring dark solitons are studied. It is predicted that characteristics of the magnetic solitons in optical lattice can be tuned by controlling the long-range light-induced and static magnetic dipole-dipole interactions; solutions of single-component magnetic and single-, two-, three-components polar solitons are found; ring dark solitons in spin-1 condensates are predicted to live longer lifetimes than that in their scalar counterparts. In the formation of spin textures, we have considered the theoretical model of a rapidly quenched and fast rotating trapped spin-1 Bose–Einstein condensate, whose dynamics can be studied by solving the stochastic projected Gross–Pitaevskii equations. Spontaneous generation of nontrivial topological defects, such as the hexagonal lattice skyrmions and square lattice of half-quantized vortices was predicted. In particular, crystallization of merons (half skyrmions) can be generated in the presence of spin–orbit coupling.

Graphical abstract

Keywords

Bose–Einstein condensate / spinor / vortex lattice / soliton / spin–orbit coupling

1 Introduction

Since the experimental realization of Bose–Einstein condensates (BECs), intensive interest in this field has been growing and the condensate has been an ideal platform for various effects of quantum many-body interaction. This opens a new field to explore various types of excitations. BECs in optical traps liberate the spin degrees of freedom of the atoms and distinguish themselves from their scalar counterparts by their much richer internal structures and dynamics  [1–4]. Thus, the ultra-cold atomic gas is composed of atoms in components with different magnetic quantum numbers and such a condensate is referred to as a spinor condensate. In addition, BECs have a high degree of experimental manipulability. Thus spinor BEC is identified as unique playground for exploring quantum phenomena related to the symmetry breaking, topological excitations, and the interplay between superfluidity and magnetism. Based on the single mode approximation, the quantum eigenstates of a spinor BEC have been extensively studied with sophisticated algebraic techniques  [5–7]. A rich variety of phenomena, such as spin domains  [8] and textures  [2] and spin-mixing  [9] arise in the condensate held by an optical potential.

At zero external magnetic field, the ground state of the spin-1 condensate is determined by the spin-dependent interaction term in the Hamiltonian [4, 10–12]. Vortices excitations have been extensively explored, such as fractional vortices [13–15], coreless vortices [3, 16–19], and Monopoles [20, 21]. Soliton dynamics and Bright-dark soliton complexes are investigated in Refs. [22–24], and the existence of multiple bright solitons was shown by proving that the system possesses a completely integrable point.

In this article, we review our recent theoretical advances in the stationary and dynamical characteristics of the spin-1 BECs. In Section 2, we derive the mean field Hamiltonian and equations from the second quantized form of the Hamiltonian to describe the stationary and dynamic characteristics of the spin-1 condensate. Section 3 is devoted to investigating the ground states of the spin-1 condensates with external parabolic potential and proposals to distinguish experimentally the ferromagnetic and polar phases of the condensate. We show the soliton solutions and dynamics of one or two dimensional spin-1 BECs in Section 4. Section 5 focuses on the nontrivial topological structure: skyrmion lattices. Finally, we summarize our results in Section 6.

2 Derivation of the mean-field model

In cold atomic vapours, particle separations are usually an order of magnitude larger than the length scales associated with the atom–atom interaction. Thus, we can approximately take account of the binary collisions between atoms. In addition, because the ultra-cold atoms have low velocities, only the s-wave scattering length a_s is necessary to describe the interaction characteristics of the colliding atoms.

The single-particle Hamiltonian of the spin-1 atoms in the second quantized form is  [25]

\begin{array}{}\text{(1)}& \begin{array}{cc}{\widehat{H}}_0=\hfill & {\displaystyle \int }\mathrm{d}r{\displaystyle \sum _{k,l=-1}^1}{\widehat{\psi }}_k^{†}(-{\displaystyle \frac{{\mathrm{\hslash }}^2{\nabla }^2}{2M}}\hfill \\ \hfill & +U(r)-p{(F_z)}_{kl}+q{(F_z^2)}_{kl})\widehat{\psi }{}_l\hfill \end{array}\end{array}

where {\widehat{\psi }}_k is the field operator with k=0,\pm 1 denoting the magnetic quantum number and M is the mass of the atom. F_z is the z-component of the spin matrix and U(r) is the one-body potential. The coefficient in the linear Zeeman term is p=-g{\mu }_{B}B, where g, {\mu }_B and B are the landé hyperfine g factor, Bohr magneton and magnetic field that is assumed to be applied in the z-direction, respectively. In the quadratic Zeeman term, q={(g{\mu }_{B}B)}^2/(E_k-E_i), where E_i and E_k are the initial and the intermediate energies, respectively. In the following, we focus on the derivation of the interaction part of the Hamiltonian in the second quantized form.

Exchange symmetry states that the many-body wave function does not change upon the exchange of two identical particles for bosons. Since only the s-wave scattering is considered, the orbital part of the wave function does not change under this exchange of two identical particles, indicating that confinements need to be exerted on the spin part of the many-body wave function, i.e., the total spin of two interacting particles must be even. To understand this issue more clearly, we can recall the spin part of the binary wave function for spin-1 atoms.

Considering two colliding atoms whose spin operators are denoted by {\widehat{F}}_1 and {\widehat{F}}_2, the total spin operator is

\begin{array}{}\text{(2)}& \widehat{F}={\widehat{F}}_1+{\widehat{F}}_2\end{array}

We denote the simultaneous eigenstates of {\widehat{F}}_1^2, {\widehat{F}}_2^2, {\widehat{F}}_{1z} and {\widehat{F}}_{2z} by |f_1,m_1;f_2,m_2⟩ (m_1,m_2=1,0,-1), which satisfy

\begin{array}{}\text{(3)}& \begin{array}{c}{\widehat{F}}_1^2|f_1,m_1;f_2,m_2⟩=2{\mathrm{\hslash }}^2|f_1,m_1;f_2,m_2⟩\hfill \\ {\widehat{F}}_2^2|f_1,m_1;f_2,m_2⟩=2{\mathrm{\hslash }}^2|f_1,m_1;f_2,m_2⟩\hfill \\ {\widehat{F}}_{1z}|f_1,m_1;f_2,m_2⟩=m_1\mathrm{\hslash }|f_1,m_1;f_2,m_2⟩\hfill \\ {\widehat{F}}_{2z}|f_1,m_1;f_2,m_2⟩=m_2\mathrm{\hslash }|f_1,m_1;f_2,m_2⟩\hfill \end{array}\end{array}

The simultaneous eigenstates of {\widehat{F}}^2 and {\widehat{F}}_1^2, {\widehat{F}}_2^2 and F_z are |f_1,f_2;f,m⟩ (f=|f_1-f_2|,\mathrm{\cdots },f_1+f_2-1,f_1+f_2 and m=-f,-f+1,\mathrm{\cdots },f), which satisfy

\begin{array}{}\text{(4)}& \begin{array}{c}{\widehat{F}}_1^2|f_1,f_2;f,m⟩=2{\mathrm{\hslash }}^2|f_1,f_2;f,m⟩\hfill \\ {\widehat{F}}_2^2|f_1,f_2;f,m⟩=2{\mathrm{\hslash }}^2|f_1,f_2;f,m⟩\hfill \\ {\widehat{F}}^2|f_1,f_2;f,m⟩=f(f+1){\mathrm{\hslash }}^2|f_1,f_2;f,m⟩\hfill \\ {\widehat{F}}_z|f_1,f_2;f,m⟩=m\mathrm{\hslash }|f_1,f_2;f,m⟩\hfill \end{array}\end{array}

The Clebsch–Gordan coefficients connect these two set of eigenstates. Here for spin-1 atoms, f=0,1,2. However, the f=1 states change sign upon the exchange of two identical atoms. Actually, for f=1, the three states with magnetic number m=1,0,-1 are

\begin{array}{}\text{(5)}& \begin{array}{c}|f_1,f_2;1,1⟩={\displaystyle \frac{1}{\sqrt{2}}}\left[|f_1,1;f_2,0⟩-|f_1,0;f_2,1⟩\right]\hfill \\ |f_1,f_2;1,0⟩={\displaystyle \frac{1}{\sqrt{2}}}\left[|f_1,1;f_2,-1⟩-|f_1,-1;f_2,1⟩\right]\hfill \\ |f_1,f_2;1,-1⟩={\displaystyle \frac{1}{\sqrt{2}}}\left[|f_1,0;f_2,-1⟩-|f_1,-1;f_2,0⟩\right]\hfill \end{array}\end{array}

Thus, the total spin of two interacting particles must be even, i.e., f=0,2. Consequently, there are two scattering lengths a₀ and a₂ to describe the scattering characteristics corresponding to these two channels of total spin f=0,2. Following Ref. [3], the form of the interaction can be written as

\begin{array}{}\text{(6)}& {\widehat{V}}_{int}=4\pi {\mathrm{\hslash }}^2/M\delta (r_1-r_2)\left(a_0{\mathcal{P}}_0+a_2{\mathcal{P}}_2\right)\end{array}

where, {\mathcal{P}}_0 and {\mathcal{P}}_2are the projection operators which project the binary interaction into two channels of total spin f=0,2, respectively. Next, we need to find ways to express the projection operators. Fortunately, we have the composition law as

\begin{array}{}\text{(7)}& F^2={(F_1+F_2)}^2\end{array}

Thus, we have

\begin{array}{}\text{(8)}& F_1\cdot F_2={\{}_{\begin{array}{cc}\mathrm{1 }& f=2\end{array}}^{\begin{array}{cc}-2& f=0\end{array}}\end{array}

Then we can express F_1\cdot F_2 in terms of the projection operators like

\begin{array}{}\text{(9)}& F_1\cdot F_2=-2{\mathcal{P}}_0+{\mathcal{P}}_2\end{array}

On the other hand, by referring to the completeness of the projection operators {\mathcal{P}}_0+{\mathcal{P}}_2=1, we finally obtain

\begin{array}{}\text{(10)}& \begin{array}{c}{\mathcal{P}}_0=\left(1-F_1\cdot F_2\right)/3\hfill \\ {\mathcal{P}}_2=\left(2+F_1\cdot F_2\right)/3\hfill \end{array}\end{array}

After the substitution of Eq. (10) into Eq. (6), we have

\begin{array}{}\text{(11)}& {\widehat{V}}_{int}=\delta (r_1-r_2)\left(g_n+g_s{F}_1\cdot F_2\right)\end{array}

where g_n=4\pi {\mathrm{\hslash }}^2(a_0+2a_2)/(3M), g_s=4\pi {\mathrm{\hslash }}^2(a_2-a_0)/(3M). Then we obtain the interaction part of the total Hamiltonian in the second quantized form

\begin{array}{}\text{(12)}& H_{int}={\displaystyle \int dr\frac{1}{2}(g_n{\widehat{\psi }}_k^{†}{\widehat{\psi }}_{k^{\prime }}^{†}{\widehat{\psi }}_{k^{\prime }}{\widehat{\psi }}_k+g_s{\widehat{\psi }}_k^{†}{\widehat{\psi }}_{k^{\prime }}^{†}{F}_{kl}\cdot F_{k^{\prime }{l}^{\prime }}{\widehat{\psi }}_{l^{\prime }}{\widehat{\psi }}_l)}\end{array}

where F=(F_x,F_y,F_z) (F_{x,y,z}are the spin-1 matrices). In the interaction part of the Hamiltonian, the Einstein summation convention is used (k,k^{\prime }=1,0,-1).

In the mean-field approximation, the Bose–Einstein condensate wavefunction {\psi }_k=⟨{\widehat{\psi }}_k⟩. Thus, the mean-field Hamiltonian of the spin-1 Bose–Einstein condensate is

\begin{array}{}\text{(13)}& \begin{array}{cc}H=\hfill & {\displaystyle \int }\mathrm{d}r[{\psi }_k^{\ast }(-{\displaystyle \frac{{\mathrm{\hslash }}^2{\nabla }^2}{2M}}+U(r)-p{(F_z)}_{kl}+q{(F_z^2)}_{kl}){\psi }_l\hfill \\ \hfill & +{\displaystyle \frac{1}{2}}(g_n{\psi }_k^{\ast }{\psi }_{k^{\prime }}^{\ast }{\psi }_{k^{\prime }}{\psi }_k+g_s{\psi }_k^{\ast }{\psi }_{k^{\prime }}^{\ast }{F}_{kl}\cdot F_{k^{\prime }{l}^{\prime }}{\psi }_{l^{\prime }}{\psi }_l)]\hfill \end{array}\end{array}

In the following sections, the time-splitting spectral method [26] is used to solve the coupled nonlinear Gross–Pitaevskii equations at zero temperature. For dynamics of the condensate at finite temperature, we solve the stochastic projected Gross–Pitaevskii equations [27–33].

3 Ground states

3.1 Ground states of the trapped condensate

The density distribution of one dimension spin-1 bosonic condensate in its ground state can be evaluated within a modified Gross–Pitaevskii theory and obvious Fermi-like distribution emerges for each component [34]. In this section we focus on pancake BECs composed of ²³Na and ⁸⁷Rb atoms. As pointed out in Refs. [3, 4], the ground state can be a polar or ferromagnetic state according to the spin-dependent interaction term g_s{\psi }_k^{\ast }{\psi }_{k^{\prime }}^{\ast }{F}_{kl}\cdot F_{k^{\prime }{l}^{\prime }}{\psi }_{l^{\prime }}{\psi }_l/2. The polar state will be realized with g_s>0, and the ferromagnetic state with .

Following the interaction parameters in Ref. [3], we choose a_0=46a_B, a_2=52a_B (a_B is the Bohr radius) for ²³Na condensate and a_0=110a_B, a_2=107a_B for ⁸⁷Rb condensate. Thus, the ground state is polar for the ²³Na condensate and ferromagnetic for the ⁸⁷Rb condensate. In Fig. 1, we show the quasi-two-dimension density distribution with a total atom number N=1\times {10}^4.

In fact, all spinors are related to each other by gauge transformation {\mathrm{e}}^{\mathrm{i}\theta } and spin rotations R(\alpha ,\beta ,\gamma ) (\alpha ,\beta ,\gamma are Euler angles), and all of them are degenerate. Thus the density configuration shown in Fig. 1 is not unique. Corresponding to the density distribution of the ⁸⁷Rb condensate, Fig. 2 shows the ferromagnetic spin structure, showing the spin magnetic moments align along x direction.

As pointed out in Refs. [35–37], we should keep in mind that nonlinear modulation instability of soliton structures in the spinor BECs would occur. Spinor BEC with ferromagnetic interactions can spontaneously generates a topological spin texture [38]. Even for the ground state, the ferromagnetic phase of the condensate may experience modulational instability, leading to a fragmentation of the spin domains. The effect of the modulational instability is shown in Fig. 3 for a pancake ferromagnetic condensate.

3.2 Condensate with spin–orbit coupling

The spin–orbit coupling is an essential mechanism for spintronics devices, and many fundamental phenomena in condensed matter physics and atomic physics have a firm relationship with it. For instance, spin–orbit coupling causes the appearance of the quantum spin Hall effect in electronic condensed matter systems [39]. Artificial external Abelian and non-Abelian gauge potentials coupled to neutral atoms have been generated recently by controlling the atom–light interaction [40–45], and this provides the possibility of spin control in neutral atom systems through laser fields. In the presence of the spin–orbit coupling, the neutral atoms in the effective abelian and non-abelian gauge fields behave like electrons in an electro–magnetic field. Thus, spin–orbit coupled ultra-cold atom systems have attracted intensive attention since last few years [46–63].

In the presence of the Rashba spin–orbit coupling, the Hamiltonian can be written as

\begin{array}{}\text{(14)}& \begin{array}{cc}H=\hfill & {\displaystyle \int }\mathrm{d}r[{\psi }_k^{\ast }(-{\displaystyle \frac{{\mathrm{\hslash }}^2{\nabla }^2}{2M}}+U(r)+{\displaystyle \frac{\mathrm{\hslash }\kappa }{M}}({\partial }_x{(F_x)}_{kl}\hfill \\ \hfill & +{\partial }_y{(F_y)}_{kl}))\psi {}_l+{\displaystyle \frac{1}{2}}(g_n{\psi }_k^{\ast }{\psi }_{k^{\prime }}^{\ast }{\psi }_{k^{\prime }}{\psi }_k\hfill \\ \hfill & +g_s{\psi }_k^{\ast }{\psi }_{k^{\prime }}^{\ast }{F}_{kl}\cdot F_{k^{\prime }{l}^{\prime }}{\psi }_{l^{\prime }}{\psi }_l)]\hfill \end{array}\end{array}

where κ characterizes the strength of the spin–orbit coupling strength, which is experimentally related to the wavelength of the laser beams and exact experimental setups. As proposed in Ref. [64], \kappa =2\pi \mathrm{\hslash }\sin (\theta /2)/\lambda, (\theta and \lambda are the angle between two Raman beams and the wavelength of the laser, respectively.). In Ref. [65], the tetrapod setup is used to generate the Rashba spin–orbit coupled BEC and it requires \kappa =\sqrt{2}\pi \mathrm{\hslash }/\lambda.

Experimentally, the quasi-two-dimension condensate can be realized by imposing a strong harmonic confinement, i.e., {\omega }_z\gg {\omega }_{\perp }={\omega }_x={\omega }_y. After the integral along z direction, the effective two dimensional Hamiltonian within the mean-field theory we can be expressed as

\begin{array}{}\text{(15)}& \begin{array}{cc}H=\hfill & {\displaystyle \int }\mathrm{d}r\{{\displaystyle \sum _{m=0,\pm 1}}[-{\displaystyle \frac{{\mathrm{\hslash }}^2}{2M}}{\psi }_m^{\ast }{\nabla }^2{\psi }_m\hfill \\ \hfill & +{\displaystyle \frac{1}{2}}M{\omega }_{\perp }^2(x^2+y^2)|{\psi }_m{|}^2]+{\displaystyle \frac{\mathrm{\hslash }\kappa }{M}}[{\psi }_1^{\ast }(-\mathrm{i}{\partial }_x-{\partial }_y){\psi }_0\hfill \\ \hfill & +{\psi }_0^{\ast }(-\mathrm{i}{\partial }_x-{\partial }_y){\psi }_{-1}+{\psi }_0^{\ast }(-\mathrm{i}{\partial }_x+{\partial }_y){\psi }_1\hfill \\ \hfill & +{\psi }_{-1}^{\ast }(-\mathrm{i}{\partial }_x+{\partial }_y){\psi }_0]+{\displaystyle \frac{g_n^{\prime }}{2}}(|{\psi }_1{|}^2+|{\psi }_0{|}^2+|{\psi }_{-1}{|}^2){}^2\hfill \\ \hfill & +{\displaystyle \frac{g_s^{\prime }}{2}}[{(|{\psi }_1{|}^2-|{\psi }_{-1}{|}^2)}^2+2|{\psi }_1^{\ast }{\psi }_0+{\psi }_0^{\ast }{\psi }_{-1}{|}^2]\}\hfill \end{array}\end{array}

where g_n^{\prime }=\sqrt{\frac{M{\omega }_z}{2\pi \mathrm{\hslash }}}{g}_n, g_s^{\prime }=\sqrt{\frac{M{\omega }_z}{2\pi \mathrm{\hslash }}}{g}_s.

The spin–orbit coupled spinor condensate can develop a plane wave phase [see Fig. 4(a)] or stripe phase [see Fig. 4(b)] due to the interaction energy [66–68]. In this condition, the characteristic length of the harmonic external potential is much greater than the wavelength of the stripe, meaning a weak external trap compared to the spin–orbit coupling strength. Trapped spin-1 BECs with neither too strong nor too weak spin–orbit coupling strength, however, are less studied. In this case, the atom–atom interaction should be large enough to produce trapped plane wave or stripe state. In the trapped plane wave or stripe state, we find there exist nontrivial detailed structures, as we will see in the following paragraphs.

The spin is accompanied by a magnetic moment, resulting in magnetism. The Zeeman field can be used to explore quantum characters of the spinor BECs, for instance, applied external Zeeman field can induce the coherent dynamics in the spin-1 BECs [69]. Zeeman field can produce a rich ground state phase diagram via tuning their magnetization [4, 70]. As pointed out in Ref. [4], a broken-axisymmetry state rises as an intermediate phase between the ferromagnetic and polar phases, showing the magnetic system responds to the external magnetic field in a very unique manner. Thus, an external Zeeman field on the spin–orbit coupled BECs.

For the quasi-two-dimensional Rashba type spin–orbit coupled spin-1 BEC within the Zeeman field, the energy functional in the mean-field approximation can be written as [71]

\begin{array}{}\text{(16)}& \begin{array}{cc}H=\hfill & {\displaystyle \int }\mathrm{d}r\{{\displaystyle \sum _m}{\psi }_m^{\ast }(-{\displaystyle \frac{{\mathrm{\hslash }}^2{\nabla }^2}{2M}}-pm+q{m}^2){\psi }_m\hfill \\ \hfill & +{\displaystyle \frac{\mathrm{\hslash }\kappa }{M}}[{\psi }_1^{\ast }(-\mathrm{i}{\partial }_x-{\partial }_y){\psi }_0+{\psi }_0^{\ast }(-\mathrm{i}{\partial }_x-{\partial }_y){\psi }_{-1}+h.c.]\hfill \\ \hfill & +{\displaystyle \frac{g_n^{\prime }}{2}}{n}^2+{\displaystyle \frac{g_s^{\prime }}{2}}[{(n_1-n_{-1})}^2+2|{\psi }_1^{\ast }{\psi }_0+{\psi }_0^{\ast }{\psi }_{-1}{|}^2]\}\hfill \end{array}\end{array}

where m=1,0,-1 and the density distribution of the component m reads n_m={|{\psi }_m|}^2 with the condensation wave function {\psi }_m, thus the total atomic density is n={\sum }_m{n}_m.

In the presence of the linear Zeeman field only, the single-particle energy spectrum is

\begin{array}{}\text{(17)}& \begin{array}{c}{E}_0=E_k\hfill \\ E_{\pm }=E_k\pm \sqrt{p^2+{\displaystyle \frac{4{\kappa }^2}{M}}{E}_k}\hfill \end{array}\end{array}

And in the presence of the quadratic Zeeman field only, we have

\begin{array}{}\text{(18)}& \begin{array}{c}{E}_0=E_k+q\hfill \\ E_{\pm }=E_k+{\displaystyle \frac{1}{2}}q\pm {\displaystyle \frac{1}{2}}\sqrt{q^2+{\displaystyle \frac{16{\kappa }^2}{M}}{E}_k}\hfill \end{array}\end{array}

where E_k={\mathrm{\hslash }}^2{k}_{\perp }^2/(2M), and the subscript \pm denotes the helicity branches of spin parallel or antiparallel to wave vector. The single-particle ground states are of the negative helicity with the momentums k_{\perp }=\sqrt{\frac{2{\kappa }^2}{{\mathrm{\hslash }}^2}-\frac{M^2{p}^2}{2{\mathrm{\hslash }}^2{\kappa }^2}} in linear Zeeman field and k_{\perp }=\sqrt{32{\kappa }^4-2q^2{M}^2}/(4\mathrm{\hslash }\kappa ) in quadratic Zeeman field, respectively. Corresponding to the negatice helicity branch, the single-particle eigenstates are given by

\begin{array}{}\text{(19)}& {\Psi }_k=\left(\begin{array}{c}{\alpha }_1\\ -{\alpha }_0{\mathrm{e}}^{\mathrm{i}\theta }\\ {\alpha }_{-1}{\mathrm{e}}^{\mathrm{i2}\theta }\end{array}\right){\mathrm{e}}^{\mathrm{i}k\cdot r}\end{array}

where \theta =\arctan (k_y/k_x), and {\alpha }_m ({\sum }_m{\alpha }_m^2=1) are

\begin{array}{}\text{(20)}& \begin{array}{cc}L:\hfill & {\alpha }_1={\displaystyle \frac{2{\kappa }^2+Mp}{4{\kappa }^2}},{\alpha }_{-1}={\displaystyle \frac{2{\kappa }^2-Mp}{4{\kappa }^2}}\hfill \\ \hfill & {\alpha }_0=\sqrt{2{\alpha }_1{\alpha }_{-1}}\hfill \\ Q:\hfill & {\alpha }_1={\alpha }_{-1}={\displaystyle \frac{\sqrt{4{\kappa }^2-Mq}}{4\kappa }}\hfill \\ \hfill & {\alpha }_0={\displaystyle \frac{\sqrt{2}\sqrt{4{\kappa }^2+Mq}}{4\kappa }}\hfill \end{array}\end{array}

where L and Q represent the energy spectrums in linear and quadratic Zeeman fields, respectively. From the single-particle energy spectrum above, we can see the Zeeman effects shift the momentum of ground state compared to the case in absence of external magnetic field.

In two dimensional homogeneous system, the ground states of the spin–orbit coupled condensate are circularly degenerate. The atomic interactions (the interaction preserves the symmetries of single particle Hamiltonian) intend to couple different single particle states and the ground state which preserves the {\mathcal{C}}_{nz} symmetry can be approximately described by the linear superposition of single-particle eigenstates on the degenerate momentum ring, \Psi ={\sum }_n{\mathcal{A}}_n{\Psi }_k with {\mathcal{A}}_n satisfying {\sum }_n{|{\mathcal{A}}_n|}^2=1 (n is an non-negative integer [67]). By substituting the ansatz into Eq. (16), it is shown that these high-symmetry states with n⩾3 are energetically unfavored. Thus, we can assume

\begin{array}{}\text{(21)}& \Psi =\mathcal{A}{\Psi }_k+\mathcal{B}{\Psi }_{-k}\end{array}

where, the coefficients \mathcal{A} and \mathcal{B} are real constants satisfying the normalization condition {\mathcal{A}}^2+{\mathcal{B}}^2=1. The approximate ground state can be obtained by minimizing the energy function in Eq. (16) after the substitution of Eq. (21).

In the presence of the Zeeman field, the boundary between plane wave and stripe phases are determined by \xi =(2-M^2{p}^2/{\kappa }^4)g_s^{\prime }+\frac{M^4{p}^4}{16{\kappa }^8}(g_s^{\prime }-g_n^{\prime }) (in the presence of the linear Zeeman field only) or \xi =2g_s^{\prime }-\frac{M^2{q}^2}{16{\kappa }^4}(2g_s^{\prime }+g_n^{\prime }) (in the presence of the quadratic Zeeman field only), where \xi =[2{\alpha }_0^2({\alpha }_1^2+6{\alpha }_1{\alpha }_{-1}+{\alpha }_{-1}^2)-{({\alpha }_1^2-{\alpha }_{-1}^2)}^2]g_s^{\prime }-{({\alpha }_1^2-{\alpha }_0^2+{\alpha }_{-1}^2)}^2{g}_n^{\prime }. In Fig. 5 (the numerical results obtained by imaginary time evolution are also shown), we show the phase diagram with respect to the parameters (g_s^{\prime },p) and (g_s^{\prime },q). We can observe that the stripe phase only exists in the case of g_s^{\prime }>0. In addition, the phase transition from the stripe to plane wave phase would occur when the strength of magnetic field surpasses the critical value p_0=2{\kappa }^2\sqrt{(\sqrt{2g_s^{\prime 2}+2g_s^{\prime }{g}_n^{\prime }}-2g_s^{\prime })/[M^2(g_n^{\prime }-g_s^{\prime })]} for the linear Zeeman field case, and q_0=4{\kappa }^2\sqrt{2g_s^{\prime }/[(g_n^{\prime }+2g_s^{\prime })M^2]} for the quadratic case.

In particular, the system can be polarized by the magnetic field if the strength of linear (or quadratic) Zeeman effect attain the critical values p_m=2{\kappa }^2/M (or q_m=4{\kappa }^2/M), and all the atoms occupy the state {\psi }_1 (or {\psi }_0) with zero momentum.

3.3 Testing the ground state

In Ref. [72], the superfluidity–Mott-insulator transition of cold atoms in optical lattices was probed by using cavity transmission spectra. As shown above in the spin-1 condensate in the absence of the spin–orbit coupling, there are two possibilities for the ground state, i.e. the ferromagnetic state and the polar state. These two phases are different in the aspect of rotational properties. The polar state is a SU(2) singlet and thus isotropic. However, the ferromagnetic state is a SU(2) multiplet heavily degenerate, and thus can be well directed. From the point view of the quantum phase transition, it is relevant to identify these two phases experimentally [73]. In this section, we focus on the method to distinguish the ferromagnetic and polar quantum ground states by using cavity quantum electrodynamics techniques.

The spin-1 BECs of certain type with total atom number N can be loaded inside an optical cavity [74, 75], and then probed with laser beam which is {\sigma }_{+} polarized. We choose the bare frequency of the cavity mode {\omega }_c to resolve the alkali D line (fine structure) and lie near the D₁ transition. In addition, we need that the detuning between the D₁ transition and the cavity mode frequency is relatively large, so that the hyperfine structures of the D₁ line would not be resolved. The effective Hamiltonian to describe the atom-cavity photon interaction can be expressed as (see the appendix in Ref. [76])

\begin{array}{}\text{(22)}& {\widehat{H}}_{int}={\displaystyle \sum _m}\left(1+{\displaystyle \frac{m}{2}}\right)U_0{\displaystyle \int dx{u}^2(x){\widehat{\psi }}_m^{†}(x){\widehat{\psi }}_m(x)a^{†}a}\end{array}

where {\widehat{\psi }}_m(x) (m=0,\pm 1) is the field annihilation operator for an atom in the hyperfine state |F=1,m⟩ at position x (the quantization direction is along the cavity axis and is denoted as z-direction), and a is the annihilation operator of the cavity photon; u(x) is the cavity mode function whose maximum is normalized to unity; U₀ describes the strength of the ac-Stark shifts of the atomic levels per photon in the viewpoint of the atoms or alternatively, the magnitude of the cavity frequency shift per atom in the viewpoint of the photons. Based upon the single mode approximation (SMA), the atomic field operators {\widehat{\psi }}_m(x)=\varphi (x)c_m [\varphi (x) denotes the common spatial mode] and c_m is the annihilation operator corresponding to the spin component with magnetic quantum number m. We can rewrite the Hamiltonian (22) as

\begin{array}{}\text{(23)}& H_{int}=U_d\left(\widehat{N}+{\displaystyle \frac{1}{2}}\widehat{M}\right)a^{†}a\end{array}

where \widehat{N}={\sum }_m{c}_m^{†}{c}_m is the total atom number and \widehat{M}={\sum }_{m}m{c}_m^{†}{c}_m the total magnetization; U_d\equiv U_0\int dx{u}^2(x){\varphi }^2(x), which can be approximated by U_d\simeq U_0/2 [76].

We assume that the frequency of the probe laser is {\omega }_p, and the Heisenberg equation of the cavity mode annihilation operator in the frame rotating at {\omega }_p is

\begin{array}{}\text{(24)}& \dot{a}=-\mathrm{i}\left[{\omega }_c-{\omega }_p+U_d\left(\widehat{N}+{\displaystyle \frac{1}{2}}\widehat{M}\right)\right]a-{\kappa }_{0}a+\eta \end{array}

where {\kappa }_0 and η are the damping rate of the cavity mode and the driving amplitude respectively. There exists stationary solution of the photon number operator and the corresponding expectation value is

\begin{array}{}\text{(25)}& n_{st}=⟨{\widehat{n}}_{st}⟩={\displaystyle \sum _{M=-N}^N}{f}_M{\displaystyle \frac{{\eta }^2}{{\kappa }_0^2+{(\Delta +U_{d}M/2)}^2}}\end{array}

where \Delta ={\omega }_c+U_{d}N-{\omega }_pis the shifted cavity-probe detuning and f_M the probability for \widehat{M} taking value M in the prescribed state satisfying {\sum }_M{f}_M=1. From Eq. (25), we can find that the cavity transmission spectra which is proportional to n_{st} contains information of the atomic statistics embodied in f_M.

Now Eq. (25) can be used to identify the two phases of the spin-1 BECs: polar and ferromagnetic phases. For the polar phase, the total spin is f=0 and we have the intracavity photon number like

\begin{array}{}\text{(26)}& n_{st}(\Delta )={\displaystyle \frac{{\eta }^2}{{\kappa }_0^2+{\Delta }^2}}\end{array}

We should notice that it is of the same shape as that of a bare cavity. For the ferromagnetic phase, we start our discussion by considering the most experimentally accessible states [77] – angular momentum coherent states, which can be obtained by rotating the z-aligned state |N,N⟩. The coherent state can be expressed as

\begin{array}{}\text{(27)}& {\mathrm{e}}^{-\mathrm{i}\theta \zeta \cdot F}|N,N⟩={\displaystyle \frac{{(x_{+1}{c}_{+1}^{†}+x_0{c}_0^{†}+x_{-1}{c}_{-1}^{†})}^N}{\sqrt{N!}}}|0⟩\end{array}

where \zeta =(-\sin \varphi ,\cos \varphi ,0)is the axis of rotation and \theta is the angle of rotation. x_i’s are determined by {\mathrm{e}}^{-\mathrm{i}\theta \zeta \cdot F}{c}_{+1}^{†}{\mathrm{e}}^{\mathrm{i}\theta \zeta \cdot F}=x_{+1}{c}_{+1}^{†}+x_0{c}_0^{†}+x_{-1}{c}_{-1}^{†}. The intracavity photon number n_{st} can be obtained as

\begin{array}{}\text{(28)}& n_{st}(\Delta )={\displaystyle \frac{{\eta }^2}{2{\kappa }_0}}{\displaystyle {\int }_{-\mathrm{\infty }}^{+\mathrm{\infty }}}dt{\mathrm{e}}^{\mathrm{i}t(\Delta -{\Delta }_c)}\times {\mathrm{e}}^{-{\kappa }_0|t|-t^2{U}_d^2{\sigma }^2/8}\end{array}

where {\Delta }_c=-M_c{U}_d/2(M_c=N\cos \theta) and \sigma =\sqrt{N}{\sigma }_0 ({\sigma }_0=\sqrt{\frac{1}{2}{\sin }^2\theta }).

From Eq. (28), we obtain the following conclusions: (1) Due to the average deviation of M from zero, the shape of n_{st} as a function of \Delta is shifted by {\Delta }_c comparing with that in Eq. (26). (2) Due to the fluctuations of M around the mean, the shape of n_{st} is broadened.

In Fig. 6, the intracavity photon number as a function of the detuning \Delta is plotted. For the polar state, the curve (solid line) is simple Lorentzian distribution centered at {\Delta }_c=0 and with width 2{\kappa }_0. For the ferromagnetic ground states, however, the center and width of the curve vary widely depending on the inclination angle θ (we have artificially shifted the curves to a common center; see the labels of the horizontal axes). Consequently, the polar and ferromagnetic ground states can be well distinguished by referring to the different shapes of the corresponding spectra.

4 Solitons

The competition between the nonlinearity and dispersion leads to the appearance of the localized wave structures, such as bright solitons (arise in the BECs with attractive interaction between atoms), dark soliton (arise in the BECs with repulsive interaction between atoms) and vortices. In this section, we will focus on this topic and have a review of our recent theoretical advances in this nonlinear matter wave state, including magnetic solitons [76], one-, two-, three-component solitons of the polar and ferromagnetic types, and ring dark solitons.

4.1 One dimensional solitons

Explorations on the spinor BECs trapped in optical lattices show that if the optical lattices are deep enough for the individual sites to be independent, the system can undergo a ferromagnetic like phase transition that leads to a “macroscopic” magnetization [79, 80]. Thus, the spinor BEC at each lattice site would behave like a spin magnet. As a consequence, the light-induced dipole–dipole interaction and the static magnetic dipole–dipole interaction cannot be ignored. These dipolar interactions can induce the ferromagnetic phase transition and spin-wave excitations [81]. In this section, we will show that the magnetic solitons would arise in the spinor BECs confined in a one-dimensional optical lattice because of the interactions between the spin waves. The characteristics of these magnetic soliton can be tuned by the light-induced and the magnetic dipole–dipole interactions.

There are three types of two-body interactions for spinor BECs trapped in optical lattices: spin-change collision, magnetic dipole–dipole interaction and light-induced dipole–dipole interaction. The Hamiltonian can be written as

\begin{array}{}\text{(29)}& \begin{array}{cc}H=\hfill & {\displaystyle \sum _n}{\displaystyle \int dr{\widehat{\psi }}_n^{†}(r)[-\frac{{\mathrm{\hslash }}^2{\nabla }^2}{2m}+V_L(r)]{\widehat{\psi }}_n(r)}\hfill \\ \hfill & +{\displaystyle \sum _{n,m,n^{\prime },m^{\prime }}}{\displaystyle \int }\mathrm{d}r\mathrm{d}{r}^{\prime }{\widehat{\psi }}_n^{†}(r){\widehat{\psi }}_m^{†}(r^{\prime })[V_{n{n}^{\prime }m{m}^{\prime }}^{coll}(r,r^{\prime })\hfill \\ \hfill & +V_{n{n}^{\prime }m{m}^{\prime }}^{d-d}(r,r^{\prime })]\widehat{\psi }{}_{m^{\prime }}(r^{\prime })\widehat{\psi }{}_{n^{\prime }}(r)+H{}_B\hfill \end{array}\end{array}

where V_L(r) is the lattice potential, and the index n,m,n^{\prime },m^{\prime }=-1,0,1 denotes the Zeeman sublevels of the ground state of the atoms. V_{n{n}^{\prime }m{m}^{\prime }}^{coll}(r,r^{\prime }) describes the two-body ground-state collisions, V_{n{n}^{\prime }m{m}^{\prime }}^{d-d}(r,r^{\prime }) include magnetic dipole–dipole interaction and light induced dipole–dipole interaction, and H_B represents the external magnetic interaction.

When the depth of the optical lattice is deep enough, the spinor atomic field operator can be expressed as {\widehat{\psi }}_m(r)={\sum }_i{\varphi }_i(r)C_m(i), where {\varphi }_i(r) is the condensate wave function for the ith microtrap and the operators {\widehat{C}}_m(i) satisfy the bosonic commutation relations [{\widehat{C}}_m(i),{\widehat{C}}_n^{†}(j)]={\delta }_{mn}{\delta }_{ij}. By using the tight-binding approximation and the Holstein–Primakoff transformation method [82], we can obtain the Hamiltonian in the first-order approximation:

\begin{array}{}\text{(30)}& \begin{array}{cc}H=\hfill & -{\gamma }_{B}NS{B}_z-J^{z}ZN{S}^2+{\displaystyle \sum _i}{\gamma }_B{B}_z{a}_i^{+}{a}_i\hfill \\ \hfill & +2Z{J}^{z}S{\displaystyle \sum _i}{a}_i^{+}{a}_i-2SJ{\displaystyle \sum _{i\delta }}(a_i^{+}{a}_{i+\delta }+a_i{a}_{i+\delta }^{+})\hfill \end{array}\end{array}

where {\gamma }_B=-{\mu }_B{g}_F, with {\mu }_B being the Bohr magneton and g_F the Landé g factor; N={\sum }_i{N}_i is the total atomic number in the lattice; S and Z are the local maximum spin and coordination number of the lattice; B_z is the magnitude of the magnetic field which is along z direction; the creation and annihilation boson operators a,a^{+} are introduced to represent the spin operators. We have assumed that J=J_{ij}=J_{ij}^{ld}-\frac{1}{4}{J}_{ij}^{md},J_{ij}^{md} is a coefficient which represent magnetic dipole-dipole interaction and J_{ij}^{ld} is a coefficient which represents light induced dipole–dipole interaction.

Then we can obtain the dispersion relation of the spin wave:

\begin{array}{}\text{(31)}& \mathrm{\hslash }{\omega }_k=-{\gamma }_B{B}_z+{\displaystyle \sum _{\delta >0}}[4J^{md}S-8SJ\cos (\delta k)]\end{array}

We can find that the dispersion relation of the spin wave can be easily tuned by the dipole–dipole interaction in the optical lattice.

We can get further on our way to include higher order terms. Thus the hamiltonian can be expressed as (omitting the fourth order term):

\begin{array}{}\text{(32)}& \begin{array}{cc}H=\hfill & {\lambda }_a^{{}^{\prime }}NS(S+1)-{\gamma }_B{B}_{z}SN+{\gamma }_B{B}_z{\displaystyle \sum _i}{a}_i^{+}{a}_i\hfill \\ \hfill & -{\displaystyle \sum _i}{\displaystyle \sum _{j\ne i}}(J_{ij}^{md}{S}^2-J_{ij}^{md}S{a}_j^{+}{a}_j-J_{ij}^{md}S{a}_i^{+}{a}_i\hfill \\ \hfill & +J_{ij}^{md}{a}_i^{+}{a}_i{a}_j^{+}{a}_j)-{\displaystyle \sum _i}{\displaystyle \sum _{j\ne i}}(2J_{ij}S{a}_i^{+}{a}_j\hfill \\ \hfill & +2J_{ij}S{a}_i{a}_j^{+})+{\displaystyle \frac{1}{2}}{\displaystyle \sum _i}{\displaystyle \sum _{j\ne i}}J{}_{ij}(a_i^{+}{a}_i^{+}{a}_i{a}_j\hfill \\ \hfill & +a_i^{+}{a}_j^{+}{a}_j{a}_j+a_j{a}_i^{+}{a}_i^{+}{a}_i+a_j^{+}{a}_j{a}_j{a}_i^{+})+\mathrm{\cdots }\hfill \end{array}\end{array}

In the following part, we will find that the hight order terms will bring into appearance the soliton type nonlinear excitations.

In the ideal case, we can assume that the spinor condensate in the optical lattice should initially be prepared in a spin-coherent state |\psi ⟩=|\{{\psi }_l\}⟩={\mathrm{\Pi }}_l|{\psi }_l⟩, with |{\psi }_l⟩=\exp (-{|{\psi }_l|}^2/2)\exp (-{\psi }_l{a}^{+})|0⟩. The vacuum state |0⟩ is the ground state of the BEC in the optical lattice, i.e., |0⟩ = |GS⟩=|N,-N⟩ [82]. The dynamics of the probability amplitude describing the nonlinear dynamics of coherent spin excitations on the lattice l can be derived as

\begin{array}{}\text{(33)}& \begin{array}{cc}\mathrm{i}\mathrm{\hslash }{\displaystyle \frac{\partial {\psi }_l}{\partial t}}=\hfill & ({\gamma }_B{B}_z+4J_{MD}S-8JS){\psi }_l+4JS(2{\psi }_l\hfill \\ \hfill & -{\psi }_{l+1}-{\psi }_{l-1})-2J{}_{MD}(|{\psi }_{l-1}{|}^2+|{\psi }_{l+1}{|}^2)\psi {}_l\hfill \\ \hfill & +J(|{\psi }_{l+1}{|}^2{\psi }_{l+1}+|{\psi }_{l-1}{|}^2{\psi }_{l-1}+2|{\psi }_l{|}^2{\psi }_{l+1}\hfill \\ \hfill & +2|{\psi }_l{|}^2{\psi }_{l-1})+J({\psi }_{l+1}^{+}+{\psi }_{l-1}^{+}){\psi }_l^2\hfill \end{array}\end{array}

The nearest-neighbor interactions are considered only, which is reasonable for the BECs in one dimension optical lattice [83]).

In the long-wavelength limit and by considering the continuum limit approximation we consider {\psi }_l,{\psi }_{l+1},{\psi }_{l-1}\to \psi (z,t). We then obtain the following soliton solutions which are characterized by parameters {\xi }_1=4JS=4S(J_{LD}-J_{MD}/4) and {\xi }_2=(4J_{MD}-8J)=6J_{MD}-8J_{LD}, where J_{MD} and J_{LD} are the average values of J_{ij}^{md} and J_{ij}^{ld}, respectively.

(i) When , , and {\xi }_2>0, there is single dark soliton solution:

\begin{array}{}\text{(34)}& \begin{array}{cc}\psi (z,t)=\hfill & \sqrt{-{\displaystyle \frac{{\xi }_2}{{\xi }_1}}}\tanh \left[-{\displaystyle \frac{{\xi }_2}{\sqrt{2}{\xi }_1}}(z-vt)\right]\hfill \\ \hfill & \cdot \exp [\mathrm{i}({\gamma }_{0}z-{\omega }_{1}t)]\hfill \end{array}\end{array}

where \mathrm{\hslash }{\omega }_1={\gamma }_B{B}_z+4J^{z}S-8JS-{\xi }_1{v}^2/4-2{\xi }_2^2/|{\xi }_1| is the single dark soliton energy, and {\gamma }_0=v/2, where v is the velocity of the soliton.

(ii) When \frac{3}{4}{J}_{MD}>J_{LD}>\frac{1}{4}{J}_{MD}, {\xi }_1>0, and {\xi }_2>0, there is bright magnetic soltion solution like

\begin{array}{}\text{(35)}& \psi (z,t)=\sqrt{{\displaystyle \frac{{\xi }_2}{{\xi }_1}}}sech\left[{\displaystyle \frac{{\xi }_2}{{\xi }_1}}(z-vt)\right]\exp [\mathrm{i}({\gamma }_{0}z-{\omega }_{2}t)]\end{array}

where \mathrm{\hslash }{\omega }_2={\gamma }_B{B}_z+4J^{z}S-8JS+{\xi }_1{v}^2/4-2{\xi }_2^2/{\xi }_1,{\gamma }_0=v/2,v is the velocity of the soliton.

For both solutions, the energy of the soliton is smaller than the spin wave. Most importantly, the width and amplitude of the above nonlinear matter waves can be directly manipulated by tuning the light induced or magnetic dipole–dipole interactions.

Next, we are going to review the soliton solutions in Ref. [35] obtained by solving the one dimension nonlinearly coupled Gross–Pitaevskii equations, which can be derived from the three dimension nonlinearly coupled Gross–Pitaevskii equations. The one dimension nonlinearly coupled Gross–Pitaevskii equations are [22, 84]

\begin{array}{}\text{(36)}& \begin{array}{cc}\mathrm{i}\mathrm{\hslash }{\partial }_t{\psi }_{\pm 1}=\hfill & -{\displaystyle \frac{{\mathrm{\hslash }}^2}{2M}}{\partial }_x^2{\psi }_{\pm 1}+(c_0+c_2)({|{\psi }_{\pm 1}|}^2+{|{\psi }_0|}^2){\psi }_{\pm 1}\hfill \\ \hfill & +(c_0-c_2){|{\psi }_{\mp 1}|}^2{\psi }_{\pm 1}+c_2{\psi }_{\mp 1}^{\ast }{\psi }_0^2\hfill \\ \mathrm{i}\mathrm{\hslash }{\partial }_t{\psi }_0=\hfill & -{\displaystyle \frac{{\mathrm{\hslash }}^2}{2M}}{\partial }_x^2{\psi }_0+(c_0+c_2)({|{\psi }_1|}^2+{|{\psi }_{-1}|}^2){\psi }_0\hfill \\ \hfill & +c_0{|{\psi }_0|}^2{\psi }_0+2c_2{\psi }_0^{\ast }{\psi }_1{\psi }_{-1}\hfill \end{array}\end{array}

where c_0=\left\{4{\mathrm{\hslash }}^2{a}_0/[M{a}_{\perp }^2(1-c{a}_0/a_{\perp })]+24{\mathrm{\hslash }}^2{a}_2/[M{a}_{\perp }^2(1-c{a}_2/a_{\perp })]\right\}/3 and c_2=\left\{4{\mathrm{\hslash }}^2{a}_2/[M{a}_{\perp }^2(1-c{a}_2/a_{\perp })]-4{\mathrm{\hslash }}^2{a}_0/[M{a}_{\perp }^2(1-c{a}_0/a_{\perp })]\right\}/3 denote the mean-field interaction and the spin-exchange interaction, respectively [85]. a_{\perp } is the characteristic length of the condensate in the transverse direction and c=-\zeta (1/2)=1.46. Introducing the transformation {({\psi }_1,{\psi }_0,{\psi }_{-1})}^{\mathrm{T}}\to {({\varphi }_1,\sqrt{2}{\varphi }_0,{\varphi }_{-1})}^{\mathrm{T}} and the time and length are measured in units of \mathrm{\hslash }/|c_0| and \sqrt{{\mathrm{\hslash }}^2/(2M)|c_0|}, respectively, then Eqs. (36) become

\begin{array}{}\text{(37)}& \begin{array}{cc}\mathrm{i}{\partial }_t{\varphi }_{\pm 1}=\hfill & -{\partial }_x^2{\varphi }_{\pm 1}-(\nu +a)({|{\varphi }_{\pm 1}|}^2+2{|{\varphi }_0|}^2){\varphi }_{\pm 1}\hfill \\ \hfill & -(\nu -a){|{\varphi }_{\mp 1}|}^2{\varphi }_{\pm 1}-2a{\varphi }_{\mp 1}^{\ast }{\varphi }_0^2\hfill \\ \mathrm{i}{\partial }_t{\varphi }_0=\hfill & -{\partial }_x^2{\varphi }_0-(\nu +a)({|{\varphi }_1|}^2+{|{\varphi }_{-1}|}^2){\varphi }_0\hfill \\ \hfill & -2\nu {|{\varphi }_0|}^2{\varphi }_0-2a{\varphi }_0^{\ast }{\varphi }_1{\varphi }_{-1}\hfill \end{array}\end{array}

where \nu =-c_0/|c_0| and a=-c_2/|c_0|.

For the particular integrable case \nu =a=1, solitons in the model of spinor BEC have been studied in Ref. [22]. For \nu \ne a, there are exact soliton solutions of polar or ferromagnetic type.

(i) Single-component ferromagnetic soliton (a+\nu >0)

\begin{array}{}\text{(38)}& {\varphi }_{-1}={\varphi }_0=0,{\varphi }_{+1}=\sqrt{{\displaystyle \frac{-2\mu }{\nu +a}}}{\mathrm{e}}^{-\mathrm{i}\mu t}\mathrm{sech}(\sqrt{-\mu }x)\end{array}

Negative chemical potential μ is the intrinsic parameter of this soliton family. From solution (38), moving soliton solutions can be obtained by means of Galilean transformation. By using linearized equations for small perturbations, this solution is demonstrated to be stable against infinitesimal perturbations.

(ii) Single-component polar soliton (\nu =+1)

\begin{array}{}\text{(39)}& {\varphi }_{\pm 1}=0,{\varphi }_0=\sqrt{-\mu }{\mathrm{e}}^{-\mathrm{i}\mu t}\mathrm{sech}(\sqrt{-\mu }x)\end{array}

By means of direct simulations with a small random perturbation being distributed uniformly between 0 and 0.03, we find that this soliton solution is always unstable in the region of .

(iii) Two-component polar soliton (\nu =+1)

\begin{array}{}\text{(40)}& {\varphi }_0=0,{\varphi }_{\pm 1}=\sqrt{-\mu }{\mathrm{e}}^{-\mathrm{i}\mu t}\mathrm{sech}(\sqrt{-\mu }x)\end{array}

Direct simulations show that the soliton is indeed unstable in the region of |a|⩾1, and stable if .

(iv) Three-component polar solitons

\begin{array}{}\text{(41)}& \begin{array}{c}{\varphi }_0={({\mu }_{+1}{\mu }_{-1})}^{1/4}\mathrm{sech}(\sqrt{-\mu }x){\mathrm{e}}^{-\mathrm{i}\mu t}\hfill \\ {\varphi }_{\pm 1}=\sqrt{-{\mu }_{\pm 1}}{\mathrm{e}}^{-\mathrm{i}\mu t}\mathrm{sech}(\sqrt{-\mu }x)\hfill \end{array}\end{array}

This three-component polar soliton is completely stable and this point has been checked in detail by seeding small random perturbations in all the components.

By comparing respective values of the Hamiltonian for these solutions when the atom number is fixed, we find that the single-component ferromagnetic soliton and the stable three-component polar soliton simultaneously provide for the minimum of energy in the case of \nu =+1 and a>0. For \nu =+1 and (attractive spin-independent and repulsive spin-exchange interactions), the two-component polar soliton plays the role of the ground state. Finally, for \nu =+1 and , there are no stable solitons (among the ones considered above).

4.2 Two dimensional solitons

Solitons in quasi-one dimension can live long lifetimes [86, 87]. But it is not the case in high dimensional cases. In two dimensional case, most types of solitons only live a short lifetime due to the dynamical instability arising from the higher dimensionality. As pointed out in Ref. [88], dark lumplike and dark stripe solitons will decay into vortex pairs under small transverse perturbations. However, there is candidate for observing long-lived dark soliton: ring dark soliton. The concept of ring dark soliton was first introduced and experimentally realized in nonlinear optics [89–91]. In this section, we will go through the dynamical generation and evolution of the ring dark solitons in F=1 ²³Na and ⁸⁷Rb condensates by modeling the phase engineering technique in the experiments [92].

Matter-wave dark solitons can be created experimentally by using the phase engineering and density engineering techniques [93–97]. As a great progress in this technique, Wright etal. [98] have realized spatially varying control of both the amplitude and the phase of the spinor order parameter.

In the experiment, a far-off resonant laser pulse is passed through an absorption plate (to modulate the spatial distribution of the intensity of the laser beam), and then space-dependent conservative Stark-shift potential can be created, leading to a space-dependent phase shift in the condensate order parameter. In Fig. 7 we show the schematic for the generation of ring dark solitons by using the phase engineering technique. The laser beam illuminated onto the condensate generate an effective potential, resulting in phase distributions in the condensate order parameter. The duration of the laser beam pulse is denoted as \delta t, which must be much smaller than the correlation time of the condensate. The phase step generated on the condensate order parameter can be adjusted by varying \delta t. The following model can be used to describe the phase distribution after the illustration of the laser beam pulse:

\begin{array}{}\text{(42)}& \varphi (x,y)={\displaystyle \frac{\mathrm{\Delta }\varphi }{2}}\times \left[1+\tanh \left({\displaystyle \frac{R_0-\sqrt{x^2+y^2}}{\mathrm{\Delta }x}}\right)\right]\end{array}

where \mathrm{\Delta }\varphi and R₀ characterize the total phase difference and the radius of the ring dark soliton, respectively. \mathrm{\Delta }x reveals the width of the potential edge (induced by the laser beam).

We choose that {\omega }_{\perp }={\omega }_x={\omega }_y=(2\pi )\times 12 Hz and {\omega }_z=(2\pi )\times 372 Hz, respectively, which is feasible experimentally. And the quasi-two dimension multicomponent Gross–Pitaevskii equations are

\begin{array}{}\text{(43)}& \begin{array}{cc}\mathrm{i}\mathrm{\hslash }{\displaystyle \frac{\partial {\psi }_1}{\partial t}}=\hfill & \left(-{\displaystyle \frac{{\mathrm{\hslash }}^2}{2M}}{\nabla }^2+U\right){\psi }_1+(g_n^{\prime }+g_s^{\prime })({|{\psi }_1|}^2+{|{\psi }_0|}^2){\psi }_1\hfill \\ \hfill & +(g_n^{\prime }-g_s^{\prime }){|{\psi }_{-1}|}^2{\psi }_1+g_s^{\prime }{\psi }_0^2{\psi }_{-1}^{\ast }\hfill \\ \mathrm{i}\mathrm{\hslash }{\displaystyle \frac{\partial {\psi }_0}{\partial t}}=\hfill & (-{\displaystyle \frac{{\mathrm{\hslash }}^2}{2M}}{\nabla }^2+U){\psi }_0+(g_n^{\prime }+g_s^{\prime })(|{\psi }_1{|}^2\hfill \\ \hfill & +|{\psi }_{-1}{|}^2)\psi {}_0+g{}_n{}^{\prime }|\psi {}_0|{}^2\psi {}_0+2g{}_s{}^{\prime }\psi {}_1\psi {}_0{}^{\ast }\psi {}_{-1}\hfill \\ \mathrm{i}\mathrm{\hslash }{\displaystyle \frac{\partial {\psi }_{-1}}{\partial t}}=\hfill & (-{\displaystyle \frac{{\mathrm{\hslash }}^2}{2M}}{\nabla }^2+U){\psi }_{-1}+(g_n^{\prime }+g_s^{\prime })(|{\psi }_{-1}{|}^2\hfill \\ \hfill & +|{\psi }_0{|}^2)\psi {}_{-1}+(g_n^{\prime }-g_s^{\prime })|\psi {}_1|{}^2\psi {}_{-1}+g{}_n{}^{\prime }\psi {}_0{}^2\psi {}_1{}^{\ast }\hfill \end{array}\end{array}

where U=\frac{1}{2}M{\Omega }^2(x^2+y^2) and \Omega ={\omega }_{\perp }/{\omega }_z. Based on the above two dimensional Gross–Pitaevskii equations (43), we are going to observe the dynamical properties of ring dark solitons in spin-1 BECs (the ground states are obtained by imaginary time evolution method).

The numerical calculations show that the lifetime of the ring dark solitons in spin-1 condensate is largely prolonged. If each component of the condensate puts forth the same parameterized ring dark solitons, then the total density distribution of the condensate exhibits the same profile as that in each component. According to Refs. [99–101] and our numerical calculations, we find that in the ⁸⁷Rb condensate, the lifetime of the nearly black ring soliton in spinor BEC without Feshbach resonance management is prolonged evidently compared to that in scalar BEC, and this lifetime is comparable with the lifetime of a shallower ring dark soliton (the total phase difference is \mathrm{\Delta }\varphi =0.707\pi) in Ref. [99]. By utilizing the Feshbach technique or varying the parameters of the external potentials and the spin-dependent, -independent interactions, it would be beneficial to prolong the lifetime of the ring dark solitons.

As for the decaying dynamics of the ring dark solitons, differences come out in ²³Na and ⁸⁷Rb condensates. We show this point in Fig. 8. At t=34.22 ms, nearly black ring dark solitons have formed in both condensates. However, ring dark solitons in the ²³Na, live a shorter lifetime compared to that in the ⁸⁷Rb condensate. Then the ring dark solitons break into x-like or +-like vortex pairs. The characteristic oscillation properties of the vortex pairs are the same as found by Ref. [100] in scalar BEC. The vortex and antivortex pairs oscillate between the x- and +-like configurations and arrange themselves on a ring that periodically shrinks and expands.

When the ring dark soliton is generated in certain components of the condensate, ring dark solitons in other components will be induced. The mutual filling of the condensates from different components supplies as effective potentials, resulting in a coexistence structure which lives hundreds of milliseconds. Ring dark solitons in both the ²³Na and ⁸⁷Rb condensates show a darkness-transferring phenomenon between different components. Figure 9 shows the dynamical process of the transferring for the ²³Na condensate.

The oscillations of the coexistence structure exhibit periodicity roughly and the oscillation frequency is about 0.58\Omega /\sqrt{2}. We should notice that this oscillation frequency is different from that found in scalar BEC [100, 102]. The periodical property can be found from Fig. 10, which shows the radius of the ring dark soliton varying with respect to the time. Line a corresponds to the case of generating the same parameterized ring dark solitons in all three components of the condensate and line b the case of generating only one obvious ring dark soliton in one component of the condensate. The oscillation frequency of line a is {\omega }_o\approx 2.14\Omega /\sqrt{2}, which is roughly consistent with that in Ref. [100] for the scalar BEC). The oscillation frequency of line b is approximately {\omega }_o\approx 0.58\Omega /\sqrt{2} for line b, indicating a prolonged oscillation period. Compared to the 1D condensate, the ring dark soliton oscillates between the maximum and minimum radius with approximately a two-times-larger frequency in the 2D condensate. For the coexistence structure composed of interdependent ring dark solitons from different components of the condensate, the oscillation shows periodic character in a rough way because there are interactions and radiation of the ring dark solitons. In fact, the embedded particles of another component intend to hinder the oscillation of the ring dark solitons, resulting in a low oscillation frequency of the ring dark solitons. Numerical calculations show that ring dark solitons suffers an earlier collapse if an obviously off-centered initial ring dark soliton is used. In addition, ring dark solitons on the coexistence structure of the ferromagnetic condensate incline to decay into vortex and antivortex pairs much earlier than that of the polar condensate.

5 Skyrmion lattices

As a common feature of the fluid systems, vortices have received extensive explorations in diverse fields of science. In fluids without viscosity, quantized vortices come up, which has been well known in liquid helium [103].

For the interactions in BECs can be accurately controlled by means of Feshbach technique [104, 105], explorations of such topological structures in BECs have certain advantages. In Ref. [3], coreless (or Skyrmion) vortices in spin-1 BECs was proposed and such topological structures have attracted much attention [16–19]. In this section we will focus on the topological structures in spin-1 BECs by means of solving the stochastic projected Gross–Pitaevskii equations, which can simulate the dynamics of the condensate at finite temperatures [33, 106, 107].

The nonequilibrium dynamics of a quenched rotating spinor BEC is described by the coupled SPGPEs,

\begin{array}{}\text{(44)}& \mathrm{d}{\psi }_j=\mathcal{P}\left[-{\displaystyle \frac{\mathrm{i}}{\mathrm{\hslash }}}{H}_j^{GP}{\psi }_j+{\displaystyle \frac{{\gamma }_j}{k_{B}T}}(\mu -H_j^{GP}){\psi }_j\right]\mathrm{d}t+\mathrm{d}{W}_j\end{array}

where

\begin{array}{}\text{(45)}& \begin{array}{cc}{H}_j^{GP}{\psi }_j=\hfill & \left[-{\displaystyle \frac{{\mathrm{\hslash }}^2{\nabla }^2}{2M}}+V(r)+g_n{\displaystyle \sum _{m=0,\pm 1}}{|{\psi }_m|}^2\right]{\psi }_j\hfill \\ \hfill & +g_s{\displaystyle \sum _{\zeta =x,y,z}}{\displaystyle \sum _{n,k,l=0\pm 1}}{(F_{\zeta })}_{jn}{(F_{\zeta })}_{kl}{\psi }_n{\psi }_k^{\ast }{\psi }_l\hfill \\ \hfill & -{\Omega }_r{L}_z{\psi }_j+{\displaystyle \sum _{\zeta =x,y}}{\displaystyle \sum _{n=0\pm 1}}{\kappa }_{\zeta }{(F_{\zeta })}_{jn}{\widehat{p}}_{\zeta }{\psi }_n\hfill \end{array}\end{array}

{\Omega }_r, L_z=-\mathrm{i}\mathrm{\hslash }(x{\partial }_y-y{\partial }_x) denote the rotation frequency and the z component of the orbital angular momentum, respectively. \widehat{p} and {\kappa }_{\zeta } are the momentum operator and the spin–orbit coupling coefficient. In Eq. (44), T, μ and {\gamma }_j denote the final temperature, the chemical potential and the growth rate of the j-th component; the complex-valued \mathrm{d}{W}_j stands for the white noise associated with the condensate growth, which is characterized by the fluctuation-dissipation relation ⟨\mathrm{d}{W}_i^{\ast }(x,t)\mathrm{d}{W}_j(x^{\prime },t)⟩=2{\gamma }_j{\delta }_{ij}{\delta }_C(x-x^{\prime })\mathrm{d}t, where {\delta }_C is the Dirac delta function for the condensate band field [27]. The projection operator \mathcal{P} projects the dynamics of the spinor BEC into the energy region below the cutoff energy E_R.

5.1 Spin textures

Topological defects may appear as a consequence of spontaneously broken symmetries. In spinor BECs, it is possible to create topological defects. The interplay between the superfluidity [108] and the magnetism leads to a variety of nontrivial states, for instance, skyrmions, Mermin–Ho vortices, Anderson–Thouless vortices and so on. Topological defects can be created in the phase transitions at finite temperatures according to the Kibble–Zurek scenario and this is because of the spontaneous symmetry breaking and thermal fluctuations around the critical point. By thermally quenching a rotating spin-1 BEC, we find topological defect crystallization of skyrmions appears. For convenience, U(r)=M{\omega }^2(x^2+y^2)/2 is considered and the oscillator units are chosen in the numerical computations, where the length, time, and energy are respectively scaled in units of \sqrt{\mathrm{\hslash }/(M\omega )}, 1/\omega and \mathrm{\hslash }\omega.

For the spinor BEC of ⁸⁷Rb, we show the time evolutions of the density profiles in Fig. 10. We can find that the condensates grow up in the evolution of the evaporative cooling. At the same time, the number of the vortices gets larger and vortices start closely binding up forming into vortex trimers in each component of the condensate. In the interwoven lattice structures, vortex core in one component of the condensate are filled with atoms coming from other components. In other words, quantized vortices avoid overlapping with each other. To characterize these structures, we calculate the incompressible kinetic energy per particle for each component [109] and find the incompressible kinetic energys of different components are different from each other.

In Fig. 11, we plot the spin textures and the topological charge distributions. From Fig. 11(a), we can find that a hexagonal lattice is formed, and Fig. 11(b) shows the spatial variation of the local spin moments parallel to the axis of rotation. Two unit cell of the spin orientations are presented in Fig. 11(c), and This vortex-like arrangement of magnetic moments is in fact the configuration of a skyrmion [110–112]. The topological charge [111] density in Fig. 11(d) shows a hexagonal lattice structure. Numerical integration over the primitive unit cell shows that each skyrmion carries a topological charge Q=1.

For the spinor BEC of ²³Na, the vortices lattice is quite different from that of the ⁸⁷Rb. In the spinor BEC of ²³Na, a square lattice of tightly bound vortex dimers emerges in each component. Figure 12 shows the corresponding spin textures and topological charge density distributions. For the spin texture, we find a mosaic of magnetic domains with staggered magnetization is formed, and almost all spins in the domains lie in the x–y plane. It should be noted that the domain walls are narrow and the spins reverse their magnetization through a Bloch wall transition. The staggered magnetic domains usually serve as smoking gun of half quantized vortices [113]. We apply the transformation {({\psi }_1,{\psi }_0,{\psi }_{-1})}^{\mathrm{T}}\to \widehat{G}(\theta )\widehat{R}(n,\chi ){({\psi }_1,{\psi }_0,{\psi }_{-1})}^{\mathrm{T}}, where \widehat{G}(\theta ) and \widehat{R}(n,\chi ) are the gauge transformation and the spin rotation, respectively. Without loss of generality, we can choose n=\cos \varphi (r)e_x+\sin \varphi (r)e_y. Numerically, we find that the spin texture will not be changed with \widehat{G}(\pi )\widehat{R}(n,\pi ). Requiring \widehat{G}(\pi )\widehat{R}(n,\pi ){({\psi }_1,{\psi }_0,{\psi }_{-1})}^{\mathrm{T}}={\mathrm{e}}^{\mathrm{i2}\pi }{({\psi }_1,{\psi }_0,{\psi }_{-1})}^{\mathrm{T}}, we can obtain {\mathrm{e}}^{-2\mathrm{i}\varphi (r)}{\psi }_{-1}={\psi }_1. In Fig. 12(d), we plot the |{\psi }_1-{\mathrm{e}}^{-2\mathrm{i}\varphi (r)}{\psi }_{-1}|, where the dark shaded region indicates the locations of the polar states while the brightest region means the positions disobeying the transformation, that is, the core of the half quantized vortices.

So we can conclude that crystallization of skyrmions and half quantized vortices will arise in the spinor BECs.

5.2 Spin–orbit coupled condensate

Spin–orbit coupled spin-1 BECs can lead to half-skyrmion-like excitations (merons), which are peculiar topological defects [114–116]. It is expected that the isolated merons can only be observed when particular boundary conditions are imposed for a meron carrying half a unit of topological charge. Merons can be generated in superfluid ³He-A contained in a rotating cylinder [117, 118] or spinor BECs with a constraining magnetic field [119]. we show in the following paragraphs that stable collective excitations (the crystalline order of merons or other isolated variants) can be created in a rapidly quenched spinor BEC with spin–orbit coupling without rotation ({\kappa }_{\zeta }\ne 0 and {\Omega }_r=0).

We begin with the condition of {\kappa }_x\ne {\kappa }_y. In Fig. 13, density and phase profiles show stripe structures in the process of the condensation. However, the stripe structure in the density profiles inclined to fade away while the system approaches equilibrium. The stripe structure in the phase profiles remains, showing that the equilibrium state is in the plane wave phase [67]. In this case, spin domains appear shortly after the quench starts and arrays of spin vortices serve as the domain walls. As the system is approaching to the equilibrium state, the domain walls are expelled outward and finally disappear. The plane wave state is stable despite of finite fluctuations in the density and spin textures. The equilibrium state favors the plane wave phase if {\kappa }_x\ne {\kappa }_y or . If {\kappa }_x={\kappa }_y=\kappa >0.8, periodic structures composed of grids of dark solitons in {\psi }_{\pm 1} components and vortices lattice in the {\psi }_0 component would appear. Corresponding to this density configuration, the spin texture consists of two interlacing square lattices of spin vortex with opposite vorticities, which we denote as the spin vortex lattice. For \kappa >0.8, the equilibrium state will end up in the plane wave phase or the spin vortex lattice phase. However, when \kappa >4 is satisfied, the equilibrium chooses the plane wave state.

The energy of the spin vortex lattice is numerically found greater than that of the plane wave phase. All the spin vortices in the spin vortex lattice exhibit similar structure: The central spin points to the z axis, and the others increasingly tilt outward and lie on the xy plane at the boundary. This is actually the Mermin–Ho vortex [16, 120, 121], or meron. A meron and its neighboring counterpart will pair up to form vortex-dipole to lower the total energy because of the ferromagnetic nature of the condensate. In addition, the density distribution of the spin vortex lattice state in the momentum space is characterized by four sharp peaks, showing the state is composed of four plane waves.

Next, we consider the spin–orbit coupled spinor BEC of ²³Na with the same parameter setting as that in the above explorations of the spinor BEC of ⁸⁷Rb. The numerical calculations show that there is no spin texture if the polar spin-1 gas is unpolarized initially. However, if the spin-1 gas is polarized, spin textures appears. Figure 14 shows the equilibrium spin texture for this spin-polarized state with a polarization ratio M/N=0.87 (M is the total magnetic momentum).

Compared to that of a meron, the spin texture in Fig. 14(a) is arranged in a quite different manner. The spins in the boundary always point to z direction. Meanwhile, the inner ones tilt gradually and eventually lie on the xy plane. The topological charge of such a configuration is Q=-1/2. In Fig. 14(c) and (d), there are three inverted merons with distorted cores. However, The spin configuration with only one inverted [shown in Fig. 14(a) and (b)] has an lower energy than that with three inverted meroms. Actually, we find the isolated inverted meron is robust against the thermal fluctuations. Thus, we conclude that crystallization of merons are expected to be emerge in the ferromagnetic spinor BEC. On the other hand, isolated inverted merons can be generated in the initially polarized polar spinor BECs, such as the spin-1 BEC of ²³Na.

In fact, half-skyrmion excitation depends on the combination of spin–orbit coupling and rotation. As found in Ref. [107], the half-skyrmions in this case originates from a dipole structure of spin which is always embedded in three vortices constructed by each condensate component, respectively. In the rotating spin–orbit coupled spin-1 BECs, the half-skyrmion (meron) can appear without pairing. For the half-skyrmion lattice in rotating spin-1 BEC, the numerical results is different from the periodic vortex lattice in Ref. [27] or the the vortices trimers in Ref. [33]. Instead, the vortices encircle the center with several circles.

In Fig. 15, we show the spin texture and the topological charge density in the spinor BEC of ⁸⁷Rb. It is shown that there are many half-skyrmions rounding the center. The topological charge density varies between -1 and 0, indicating that our results in Fig. 15(a) are half-skyrmions but not meron-antimeron pairs. We conclude that the half-skyrmions (meron) can occur but do not in the form of the meron-antimeron pairs.

6 Conclusions

In summary, we have investigated both the static and dynamic matter wave structures in spin-1 BECs. The ground state of the spin-orbit coupled spinor condensate is investigated in the presence of the external Zeeman field. By using the cavity quantum electrodynamics techniques, the ferromagnetic and polar quantum ground states are possible to be distinguished. For the polar state, the intracavity photon number (as a function of detuning) is simple Lorentzian distribution; for the ferromagnetic ground states, however, the center and width of the curve vary widely depending on the inclination angle. Many kinds of soliton solutions of the one dimensional spin-1 condensate are found, and their stability is analyzed. Two dimensional ring dark solitons are studied, and they are found live longer lifetimes than that in the scalar BECs. In addition, the ring dark soliton in one component of the condensate inclines to induce ring dark solitons in other components, forming into long lived coexistence structures, which show quasi-periodic oscillations. Skyrmion lattices formed in rotating spin-1 BECs or spin–orbit coupled spinor BECs are investigated. It is predicted that crystallization of skyrmions and merons will arise in the spinor BECs.

References

[1] D. M. Stamper-Kurn, M. R. Andrews, A. P. Chikkatur, S. Inouye, H. J. Miesner, J. Stenger, and W. Ketterle, Phys. Rev. Lett., 1998, 80(10): 2027

[2] T. Ohmi and K.  Machida, J. Phys. Soc. Jpn., 1998, 67(6): 1822

[3] T. L. Ho, Phys. Rev. Lett., 1998, 81(4): 742

[4] J. Stenger, S. Inouye, D. M. Stamper-Kurn, H. J. Miesner, A. P. Chikkatur, and W. Ketterle, Nature, 1998, 396(6709): 345

[5] C. K. Law, H. Pu, and N. P. Bigelow, Phys. Rev. Lett., 1998, 81(24): 5257

[6] M. Koashi and M.  Ueda, Phys. Rev. Lett., 2000, 84(6): 1066

[7] T. L. Ho and L.  Yin, Phys. Rev. Lett., 2000, 84(11): 2302

[8] H. J. Miesner, D. M. Stamper-Kurn, J. Stenger, S. Inouye, A. P. Chikkatur, and W. Ketterle, Phys. Rev. Lett., 1999, 82(11): 2228

[9] H. Pu, C. K. Law, S. Raghavan, J. H. Eberly, and N. P. Bigelow, Phys. Rev. A, 1999, 60(2): 1463

[10] M. S. Chang, C. D. Hamley, M. D. Barrett, J. A. Sauer, K. M. Fortier, W. Zhang, L. You, and M. S. Chapman, Phys. Rev. Lett., 2004, 92(14): 140403

[11] A. Widera, F. Gerbier, S. Fölling, T. Gericke, O. Mandel, and I. Bloch, New J. Phys., 2006, 8(8): 152

[12] N. Klausen, J. Bohn, and C. Greene, Phys. Rev. A, 2001, 64(5): 053602

[13] F. Zhou, Phys. Rev. Lett., 2001, 87(8): 080401

[14] G. W. Semenoff and F.  Zhou, Phys. Rev. Lett., 2007, 98(10): 100401

[15] T. Isoshima, K. Machida, and T. Ohmi, J. Phys. Soc. Jpn., 2001, 70(6): 1604

[16] T. Mizushima, K. Machida, and T. Kita, Phys. Rev. Lett., 2002, 89(3): 030401

[17] U. Leonhardt and G. E.  Volovik, JETP Lett., 2000, 72(2): 46

[18] T. Isoshima and K.  Machida, Phys. Rev. A, 2002, 66(2): 023602

[19] T. Isoshima, K. Machida, and T. Ohmi, J. Phys. Soc. Jpn., 2001, 70(6): 1604

[20] H. T. C. Stoof, E. Vliegen, and U. Al Khawaja, Phys. Rev. Lett., 2001, 87(12): 120407

[21] J. Ruostekoski and J. R.  Anglin, Phys. Rev. Lett., 2003, 91(19): 190402

[22] J. Ieda, T. Miyakawa, and M. Wadati, Phys. Rev. Lett., 2004, 93(19): 194102

[23] E. V. Doktorov, J. Wang, and J. Yang, Phys. Rev. A, 2008, 77(4): 043617

[24] H. E. Nistazakis, D. J. Frantzeskakis, P. G. Kevrekidis, B. A. Malomed, and R. Carretero-González, Phys. Rev. A, 2008, 77(3): 033612

[25] M. Ueda and Y.  Kawaguchi, Phys. Rep., 2012, 520: 253

[26] W. Z. Bao and Y. Z.  Zhang, Methods Appl. Anal., 2010, 17: 049

[27] A. S. Bradley, C. Gardiner, and M. Davis, Phys. Rev. A, 2008, 77(3): 033616

[28] S. J. Rooney, A. S. Bradley, and P. B. Blakie, Phys. Rev. A, 2010, 81(2): 023630

[29] S. J. Rooney, P. B. Blakie, B. P. Anderson, and A. S. Bradley, Phys. Rev. A, 2011, 84(2): 023637

[30] A. Sinatra, C. Lobo, and Y. Castin, J. Phys. B, 2002, 35(17): 3599

[31] C. W. Gardiner and P.  Zoller, Phys. Rev. A, 1998, 58(1): 1050

[32] C. W. Gardiner and P.  Zoller, Phys. Rev. A, 2000, 61(3): 033601

[33] S. W. Su, C. H. Hsueh, I. K. Liu, T. L. Horng, Y. C. Tsai, S. C. Gou, and W. M. Liu, Phys. Rev. A, 2011, 84(2): 023601

[34] Y. Hao, Y. Zhang, J. Q. Liang, and S. Chen, Phys. Rev. A, 2006, 73(5): 053605

[35] L. Li, Z. Li, B. A. Malomed, D. Mihalache, and W. M. Liu, Phys. Rev. A, 2005, 72(3): 033611

[36] N. P. Robins, W. Zhang, E. A. Ostrovskaya, and Y. S. Kivshar, Phys. Rev. A, 2001, 64(2): 021601

[37] W. Zhang, D. L. Zhou, M. S. Chang, M. S. Chapman, and L. You, Phys. Rev. Lett., 2005, 95(18): 180403

[38] H. Saito, Y. Kawaguchi, and M. Ueda, Phys. Rev. Lett., 2006, 96(6): 065302

[39] S. Murakami, N. Nagaosa, and S. C. Zhang, Phys. Rev. B, 2004, 69(23): 235206

[40] Y. J. Lin, R. L. Compton, A. R. Perry, W. D. Phillips, J. V. Porto, and I. B. Spielman, Phys. Rev. Lett., 2009, 102(13): 130401

[41] Y. J. Lin, R. L. Compton, K. Jiménez-García, W. D. Phillips, J. V. Porto, and I. B. Spielman, Nat. Phys., 2011, 7(7): 531

[42] Y. J. Lin, K. Jiménez-García, and I. B. Spielman, Nature, 2011, 471(7336): 83

[43] M. Aidelsburger, M. Atala, S. Nascimbene, S. Trotzky, Y. A. Chen, and I. Bloch, Phys. Rev. Lett., 2011, 107(25): 255301

[44] Z. Fu, P. Wang, S. Chai, L. Huang, and J. Zhang, Phys. Rev. A, 2011, 84(4): 043609

[45] D. L. Campbell, G. Juzeliūnas, and I. B. Spielman, Phys. Rev. A, 2011, 84(2): 025602

[46] T. D. Stanescu, B. Anderson, and V. Galitski, Phys. Rev. A, 2008, 78(2): 023616

[47] S. Gopalakrishnan, A. Lamacraft, and P. M. Goldbart, Phys. Rev. A, 2011, 84(6): 061604

[48] S. Sinha, R. Nath, and L. Santos, Phys. Rev. Lett., 2011, 107(27): 270401

[49] H. Hu, B. Ramachandhran, H. Pu, and X. J. Liu, Phys. Rev. Lett., 2012, 108(1): 010402

[50] J. Radić, T. A. Sedrakyan, I. B. Spielman, and V. Galitski, Phys. Rev. A, 2011, 84(6): 063604

[51] H. Hu and X. J.  Liu, Phys. Rev. A, 2012, 85(1): 013619

[52] D. W. Zhang, Z. Y. Xue, H. Yan, Z. D. Wang, and S. L. Zhu, Phys. Rev. A, 2012, 85(1): 013628

[53] B. Ramachandhran, B. Opanchuk, X. J. Liu, H. Pu, P. D. Drummond, and H. Hu, Phys. Rev. A, 2012, 85(2): 023606

[54] R. Barnett, S. Powell, T. Grass, M. Lewenstein, and S. Das-Sarma, Phys. Rev. A, 2012, 85(2): 023615

[55] W. Zheng and Z. B.  Li, Phys. Rev. A, 2012, 85(5): 053607

[56] C. M. Jian and H.  Zhai, Phys. Rev. B, 2011, 84(6): 060508

[57] Y. X. Du, H. Yan, D. W. Zhang, C. J. Shan, and S. L. Zhu, Phys. Rev. A, 2012, 85(4): 043619

[58] M. Merkl, A. Jacob, F. E. Zimmer, P. Öhberg, and L. Santos, Phys. Rev. Lett., 2010, 104(7): 073603

[59] Y. P. Zhang, L. Mao, and C. W. Zhang, Phys. Rev. Lett., 2012, 108(3): 035302

[60] Y. Deng, J. Cheng, H. Jing, C. P. Sun, and S. Yi, Phys. Rev. Lett., 2012, 108(12): 125301

[61] X. Q. Xu and J. H.  Han, Phys. Rev. Lett., 2011, 107(20): 200401

[62] X. F. Zhou, J. Zhou, and C. J. Wu, Phys. Rev. A, 2011, 84(6): 063624

[63] Z. F. Xu, R. Lü, and L. You, Phys. Rev. A, 2011, 83(5): 053602

[64] P. Wang, Z. Q. Yu, Z. Fu, J. Miao, L. Huang, S. Chai, H. Zhai, and J. Zhang, Phys. Rev. Lett., 2012, 109(9): 095301

[65] G. Juzeliūnas, J. Ruseckas, and J. Dalibard, Phys. Rev. A, 2010, 81(5): 053403

[66] T. L. Ho and S. Z.  Zhang, Phys. Rev. Lett., 2011, 107(15): 150403

[67] C. J. Wang, C. Gao, C. M. Jian, and H. Zhai, Phys. Rev. Lett., 2010, 105(16): 160403

[68] S. K. Yip, Phys. Rev. A, 2011, 83(4): 043616

[69] M. S. Chang, Q. S. Qin, W. X. Zhang, L. You, and M. S. Chapman, Nat. Phys., 2005, 1(2): 111

[70] K. Murata, H. Saito, and M. Ueda, Phys. Rev. A, 2007, 75(1): 013607

[71] L. Wen, Q. Sun, H. Q. Wang, A. C. Ji, and W. M. Liu, Phys. Rev. A, 2012, 86(4): 043602

[72] I. B. Mekhov, C. Maschler, and H. Ritsch, Nat. Phys., 2007, 3(5): 319

[73] S. Ashhab and A. J.  Leggett, Phys. Rev. A, 2002, 65(2): 023604

[74] F. Brennecke, T. Donner, S. Ritter, T. Bourdel, M. Köhl, and T. Esslinger, Nature, 2007, 450(7167): 268

[75] Y. Colombe, T. Steinmetz, G. Dubois, F. Linke, D. Hunger, and J. Reichel, Nature, 2007, 450(7167): 272

[76] J. M. Zhang, S. Cui, H. Jing, D. L. Zhou, and W. M. Liu, Phys. Rev. A, 2009, 80(4): 043623

[77] J. Kronjäger, C. Becker, P. Navez, K. Bongs, and K. Sengstock, Phys. Rev. Lett., 2006, 97(11): 110404

[78] Z. W. Xie, W. Zhang, S. T. Chui, and W. M. Liu, Phys. Rev. A, 2004, 69(5): 053609

[79] H. Pu, W. P. Zhang, and P. Meystre, Phys. Rev. Lett., 2001, 87(14): 140405

[80] K. Cross, C. P. Search, H. Pu, W. P. Zhang, and P. Meystre, Phys. Rev. Lett., 2002, 88: 033603

[81] W. Zhang, H. Pu, C. P. Search, and P. Meystre, Phys. Rev. Lett., 2002, 88(6): 060401

[82] D. I. Pushkarov and Kh. I.  Pushkarov, Phys. Lett. A, 1977, 61(5): 339

[83] V. V. Konotop, M. Salerno, and S. Takeno, Phys. Rev. E, 1997, 56(6): 7240

[84] T. Bergeman, M. G. Moore, and M. Olshanii, Phys. Rev. Lett., 2003, 91(16): 163201

[85] M. Olshanii, Phys. Rev. Lett., 1998, 81(5): 938

[86] C. Becker, S. Stellmer, P. Soltan-Panahi, S. Dörscher, M. Baumert, E. Richter, J. Kronjäger, K. Bongs, and K. Sengstock, Nat. Phys., 2008, 4(6): 496

[87] A. Weller, J. P. Ronzheimer, C. Gross, J. Esteve, M. K. Oberthaler, G. Theocharis, and P. G. Kevrekidis, Phys. Rev. Lett., 2008, 101(13): 130401

[88] G. X. Huang, V. A. Makarov, and M. G. Velarde, Phys. Rev. A, 2003, 67(2): 023604

[89] Y. S. Kivshar and X.  Yang, Phys. Rev. E, 1994, 50(1): R40

[90] H. E. Nistazakis, D. J. Frantzeskakis, B. A. Malomed, and P. G. Kevrekidis, Phys. Lett. A, 2001, 285(3–4): 157

[91] A. Dreischuh, D. Neshev, G. G. Paulus, F. Grasbon, and H. Walther, Phys. Rev. E, 2002, 66(6): 066611

[92] S. W. Song, D. S. Wang, H. Wang, and W. M. Liu, Phys. Rev. A, 2012, 85(6): 063617

[93] L. Dobrek, M. Gajda, M. Lewenstein, K. Sengstock, G. Birkl, and W. Ertmer, Phys. Rev. A, 1999, 60(5): R3381

[94] J. Denschlag, J. E. Simsarian, D. L. Feder, C. W. Clark, L. A. Collins, J. Cubizolles, L. Deng, E. W. Hagley, K. Helmerson, W. P. Reinhardt, S. L. Rolston, B. I. Schneider, and W. D. Phillips, Science, 2000, 287(5450): 97

[95] Z. Dutton, M. Budde, C. Slowe, and L. V. Hau, Science, 2001, 293(5530): 663

[96] N. S. Ginsberg, J. Brand, and L. V. Hau, Phys. Rev. Lett., 2005, 94(4): 040403

[97] I. Shomroni, E. Lahoud, S. Levy, and J. Steinhauer, Nat. Phys., 2009, 5(3): 193

[98] K. C. Wright, L. S. Leslie, A. Hansen, and N. P. Bigelow, Phys. Rev. Lett., 2009, 102(3): 030405

[99] X. H. Hu, X. F. Zhang, D. Zhao, H. G. Luo, and W. M. Liu, Phys. Rev. A, 2009, 79(2): 023619

[100] G. Theocharis, D. J. Frantzeskakis, P. G. Kevrekidis, B. A. Malomed, and Y. S. Kivshar, Phys. Rev. Lett., 2003, 90(12): 120403

[101] X. F. Zhang, X. H. Hu, D. S. Wang, X. X. Liu, and W. M. Liu, Front. Phys., 2011, 6(1): 46

[102] V. V. Konotop and L.  Pitaevskii, Phys. Rev. Lett., 2004, 93(24): 240403

[103] L. P. Pitaevskii and E. M.  Lifshitz, Statistical Physics (part 2), 3rd Ed., Oxford: Pergamon, 1980, section 29

[104] S. Inouye, M. Andrews, J. Stenger, H. Miesner, D. Stamper-Kurn, and W. Ketterle, Nature, 1998, 392: 151

[105] C. Chin, R. Grimm, P. Julienne, and E. Tiesinga, Rev. Mod. Phys., 1998, 70: 1003

[106] S. W. Su, I. K. Liu, Y. C. Tsai, W. M. Liu, and S. C. Gou, Phys. Rev. A, 2012, 86(2): 023601

[107] C. F. Liu and W. M.  Liu, Phys. Rev. A, 2012, 86(3): 033602

[108] A. Y. Cherny, J. S. Caux, and J. Brand, Front. Phys., 2012, 7(1): 54

[109] T. L. Horng, C. H. Hsueh, S. W. Su, Y. M. Kao, and S. C. Gou, Phys. Rev. A, 2009, 80(2): 023618

[110] U. K. Rößler, A. N. Bogdanov, and C. Pfleiderer, Nature, 2006, 442(7104): 797

[111] S. Mühlbauer, B. Binz, F. Jonietz, C. Pfleiderer, A. Rosch, A. Neubauer, R. Georgii, and P. Boni, Science, 2009, 323(5916): 915

[112] X. Z. Yu, Y. Onose, N. Kanazawa, J. H. Park, J. H. Han, Y. Matsui, N. Nagaosa, and Y. Tokura, Nature, 2010, 465(7300): 901

[113] D. A. Ivanov, Phys. Rev. Lett., 2001, 86(2): 268

[114] A. Actor, Rev. Mod. Phys., 1979, 51(3): 461

[115] L. Brey, H. A. Fertig, R. Cote, and A. H. MacDonald, Phys. Rev. B, 1996, 54(23): 16888

[116] G. E. Brown and M.  Rho, The Multifaced Skyrmion, Singapore: World Scientific, 2010

[117] V. M. H. Ruutu, J. Kopu, M. Krusius, U. Parts, B. Placais, E. V. Thuneberg, and W. Xu, Phys. Rev. Lett., 1997, 79(25): 5058

[118] R. Ishiguro, O. Ishikawa, M. Yamashita, Y. Sasaki, K. Fukuda, M. Kubota, H. Ishimoto, R. Packard, T. Takagi, T. Ohmi, and T. Mizusaki, Phys. Rev. Lett., 2004, 93(12): 125301

[119] A. E. Leanhardt, Y. Shin, D. Kielpinski, D. E. Pritchard, and W. Ketterle, Phys. Rev. Lett., 2003, 90(14): 140403

[120] T. Mizushima, K. Machida, and T. Kita, Phys. Rev. A, 2002, 66(5): 053610

[121] T. Kita, T. Mizushima, and K. Machida, Phys. Rev. A, 2002, 66: 061601(R)

Rights & permissions

© Higher Education Press and Springer-Verlag Berlin Heidelberg 2013