1. Introduction
Realization of spinor Bose–Einstein condensates (BECs) opens a new research area in ultracold atomic systems. The idea of spinor condensations was proposed by Ho and Ohmi.[1–4] In 1998, BECs were realized in a system of spin-1 23Na confined in an optical dipole trap in experiment.[5] The directions of the atomic spins may change due to the interparticle interaction, so the order parameters of a spin-f BEC have 2f + 1 components that can vary on space and time, which present rich variety of spin textures. In addition, there are some new phenomena that are not present in single-component and two-component BECs, such as ground-state phase diagram, spin dynamics, and spin-exchanged collision.[6,7] Such novel properties make the spinor BEC as an attractive area of research for quantum gases. At present, many studies have been carried out.[8–14] However, there is few work on vortexes of spinor BECs with time–space modulated interaction presenting.
Quantum vortexes are topological singular points of the quantum fluid wave function. The phase of the wave function will have a periodic change around these singular points. It can be observed in many fields of physics, such as super-fluids and optical fields. Since BECs in a dilute atomic gas were observed in 1995, experimenters have been searching a method to create vortexes in this system. At present, there are many methods to create vortexes, for example, using a focused laser beams to stir a BEC of 87Rb confined in a magnetic trap can produce a vortex.[15] The dynamics of quantized vortexes is essential for understanding different superfluid phenomena such as novel quantum phases[16] and quantum turbulence.[17] What’s more, vortexes are also topological defects that play a vital role in coherent and dissipative properties of superfluid systems.[18]
Inspired by importance of vortices and our understanding of dynamical behaviors of nonlinear excitation in physical systems, we make the theoretical analysis and numerical studies of the spinor BECs with time–space modulated interactions. Recent experiments show that the effective scattering length can be tuned by Feschbach resonance.[19] This brought about a good proposal for manipulation of nonlinear excitations and matter waves by controlling the time-dependent or space-dependent scattering strength.[20–26] As for spinor BECs, one can use an optically induced Feshbach resonance[27] or a confinement induced resonance[28] to tune nonlinearities. In this paper, two families of exact vortex solutions of spinor BECs with time–space modulated interactions are constructed explicitly by similarity transformation. Physical interpretations and numerical simulations of the obtained results are illustrated. Stable and unstable vortices for different topological charges and quantum numbers are shown by tuning the frequency of external potential and spatiotemporally interaction.
2. The model and analytic vortex solutions
Here we focus on BECs of alkali atoms 87Rb in the F = 1 hyperfine state, restricted to two-dimensional space by purely optical means. Under the mean-field approximation, the dynamics of the spinor condensates can be described by the three-component Gross–Pitaevskii equation (GPE) with an angular momentum rotational term
$$ \begin{eqnarray}\begin{array}{lll}{\rm{i}}\displaystyle \frac{\partial {\psi }_{\pm 1}}{\partial t} & = & \left(-\displaystyle \frac{{\nabla }^{2}}{2}+({g}_{{\rm{n}}}+{g}_{{\rm{s}}})(|{\psi }_{\pm 1}{|}^{2}+|{\psi }_{0}{|}^{2})\right.\\ & & +\left.({g}_{{\rm{n}}}-{g}_{{\rm{s}}})|{\psi }_{\mp 1}{|}^{2}+{V}_{{\rm{ext}}}+{E}_{\pm 1}-\varOmega {L}_{z}\right){\psi }_{\pm 1}\\ & & +{g}_{{\rm{s}}}{\psi }_{\mp 1}^{* }{\psi }_{0}^{2},\\ {\rm{i}}\displaystyle \frac{\partial {\psi }_{0}}{\partial t} & = & \left(-\displaystyle \frac{{\nabla }^{2}}{2}+({g}_{{\rm{n}}}+{g}_{{\rm{s}}})(|{\psi }_{1}{|}^{2}+|{\psi }_{-1}{|}^{2})+{g}_{n}|{\psi }_{0}{|}^{2}\right.\\ & & \left.+{V}_{{\rm{ext}}}+{E}_{0}-\varOmega {L}_{z}\right){\psi }_{0}+2{g}_{{\rm{s}}}{\psi }_{-1}{\psi }_{0}^{* }{\psi }_{1}.\end{array}\end{eqnarray}$$ (1)
Here,
is the dimensionless wave function of the rotating spin-1 BECs and Ej ∈ R is the dimensionless Zeeman energy of spin component mF = –1,0,1. Vext is the external trapping potential and it can be assumed as with m being the mass of 87Rb atoms, ωr and ωz the confining frequencies in the transverse and axial directions. We assume that the harmonic trapping frequencies satisfy ωz ≫ ω⊥, then the condensates are pressed into a pancake. The harmonic trapping potential is reduced to its 2D form . Lz is the z-component of the dimensionless angular momentum rotation defined as Lz = –i (x∂y – y∂x). The strength of the interaction is given by , gs = ± 1; a0 and a2 are the s-wave scattering lengths of scattering channels with total hyperfine spin-0 and spin-2, respectively.[15] The interaction will be ferromagnetic if (such as 87Rb), and antiferromagnetic if (such as 23Na). Specially, gn = 216, gs = –1 for 87Rb and gn = 32, gs = 1 for 23Na. In our case, the scattering lengths depend on time and space, that is to say, gn = gn(x,y,t), gs = gs(x,y,t), which can be realized by controlling the induced Feschbach resonances optically or confinement induced resonances in the real BECs experiments.[27–29] The length and time are measured in units of and , respectively.Introducing the polar coordinate transformation x = r cos θ, y = r sinθ, where θ is the azimuthal angel. Equation (1) can be rewritten as
$$ \begin{eqnarray}\begin{array}{ll}{\rm{i}}\displaystyle \frac{\partial {\psi }_{\pm 1}}{\partial t}= & -\displaystyle \frac{1}{2}{\psi }_{\pm 1rr}-\displaystyle \frac{1}{2r}{\psi }_{\pm 1r}-\displaystyle \frac{1}{{r}^{2}}{\psi }_{\pm 1\theta \theta }-{\rm{i}}\varOmega {\psi }_{\pm 1\theta }\\ & +[({g}_{{\rm{n}}}+{g}_{{\rm{s}}})(|{\psi }_{\pm 1}{|}^{2}+|{\psi }_{0}{|}^{2})+({g}_{{\rm{n}}}-{g}_{{\rm{s}}})|{\psi }_{\mp 1}{|}^{2}\\ & +{V}_{{\rm{ext}}}+{E}_{1}-\varOmega {L}_{z}]{\psi }_{\pm 1}+{g}_{{\rm{s}}}{\psi }_{\mp 1}^{* }{\psi }_{0}^{2},\\ {\rm{i}}\displaystyle \frac{\partial {\psi }_{0}}{\partial t}= & -\displaystyle \frac{1}{2}{\psi }_{0rr}-\displaystyle \frac{1}{2r}{\psi }_{0r}-\displaystyle \frac{1}{{r}^{2}}{\psi }_{1\theta \theta }-{\rm{i}}\varOmega {\psi }_{0\theta }+[({g}_{{\rm{n}}}+{g}_{{\rm{s}}})\\ & \times (|{\psi }_{1}{|}^{2}+|{\psi }_{-1}{|}^{2})+{g}_{n}|{\psi }_{0}{|}^{2}+{V}_{{\rm{ext}}}+{E}_{0}-\varOmega {L}_{z}]{\psi }_{0}\\ & +2{g}_{{\rm{s}}}{\psi }_{-1}{\psi }_{0}^{* }{\psi }_{1}.\end{array}\end{eqnarray}$$ (2)
In order to obtain the vortex solutions to Eq. (2) for limr → ∞ψ±1,0 = 0, we assume the solution in the form
$$ \begin{eqnarray}{\varPsi }_{\pm 1,0}={\beta }_{\pm 1,0}(r,t){U}_{\pm 1,0}(X(r,t)){{\rm{e}}}^{{\rm{i}}\alpha \theta +i{\phi }_{\pm 1,0}(r,t)},\end{eqnarray}$$ (3)
where α±1,0 denotes the topological charge relevant to the angular momentum and , β±1,0(r,t) is the amplitude of each wave, and X(r,t) denotes the intermediate variable reflecting the charges of the wave function U±1,0. Substituting Eq. (3) into Eq. (2) yields the ordinary differential equations (ODEs)
$$ \begin{eqnarray}\begin{array}{lll} & & {U}_{1XX}+{b}_{11}{U}_{1}^{3}+{b}_{12}{U}_{1}{U}_{-1}^{2}+{b}_{13}{U}_{0}^{2}{U}_{1}+{b}_{14}{U}_{-1}{U}_{0}^{2}=0,\\ & & {U}_{-1,XX}+{b}_{21}{U}_{1}^{2}{U}_{-1}+{b}_{22}{U}_{-1}^{3}+{b}_{23}{U}_{0}^{2}{U}_{-1}+{b}_{24}{U}_{1}{U}_{0}^{2}=0,\\ & & {U}_{0XX}+{b}_{31}{U}_{-1}^{2}{U}_{0}+{b}_{32}{U}_{1}^{2}{U}_{0}+{b}_{33}{U}_{0}^{3}+{b}_{34}{U}_{1}{U}_{-1}{U}_{0}=0,\end{array}\end{eqnarray}$$ (4)
and the partial differential equations (PDEs) related to β±1,0(r,t), ϕ±1(r,t) and X(r,t)
$$ \begin{eqnarray}\begin{array}{\l}{r}^{2}{\beta }_{irr}+r{\beta }_{ir}-(2{\omega }^{2}{r}^{4}+2{E}_{i}{r}^{2}+2{r}^{2}{\phi }_{it}\\ \,+{r}^{2}{\phi }_{ir}^{2}+2{\alpha }_{i}^{2}+2\varOmega {r}^{2}{\alpha }_{i}){\beta }_{i}=0,\\ r{\beta }_{i}{X}_{r,r}+(2r{\beta }_{ir}+{\beta }_{i}){X}_{r}=0,\\ 2r{\beta }_{it}+2r{\phi }_{ir}{\beta }_{ir}+(r{\phi }_{irr}+{\phi }_{ir}){\beta }_{i}=0,\\ {X}_{t}+{\phi }_{ir}{X}_{r}=0,\,(i=1,-1)\\ 4{r}^{2}{\beta }_{0rr}+4r{\beta }_{0r}-[2{\alpha }_{1}^{2}+2{\alpha }_{2}^{2}+{r}^{2}({\phi }_{1r}^{2}+{\phi }_{2r}^{2})+2{r}^{2}{\phi }_{1r}{\phi }_{2r}\\ \,+4\varOmega {r}^{2}({\alpha }_{1}+{\alpha }_{2})+4{r}^{2}({\phi }_{1t}+{\phi }_{2t})\\ \,+8{\omega }^{2}{r}^{4}+8{E}_{0}{r}^{2}+4{\alpha }_{1}{\alpha }_{2}]{\beta }_{0}=0,\\ r{\beta }_{0}{X}_{r,r}+(2r{\beta }_{0r}+{\beta }_{0}){X}_{r}=0,\\ 2r({\phi }_{1r}+{\phi }_{2r}){\beta }_{0r}+4r{\beta }_{0t}+(r{\phi }_{1rr}+r{\phi }_{2rr}+{\phi }_{1r}+{\phi }_{2r}){\beta }_{0}=0,\\ ({\phi }_{1r}+{\phi }_{2r}){X}_{r}+2{X}_{t}=0,\end{array}\end{eqnarray}$$ (5)
where bi,j (i = 1,2,3,j = 1,2,3,4) are constants and satisfy b14 = 2(b11 – b13) and 2b33 = b11 + b13.To obtain explicit solutions, we choose ω2 as the form
$$ \begin{eqnarray}{\omega }^{2}={\omega }_{r}/{\omega }_{z}={\omega }_{0}^{2}+\epsilon \cos (\bar{\omega }t),\end{eqnarray}$$ (6)
where ϵ ∈ (–1,1) and ω0, . Under the condition E1 + E–1 = E0, solving ODEs (4) gives
$$ \begin{eqnarray}\begin{array}{lll} & & {U}_{\pm 1,0}^{(1)}(X)={C}_{\pm 1,0}{v}_{0}CN({v}_{0}X+K(\sqrt{2}/2),\sqrt{2}/2),\\ & & \end{array}\end{eqnarray}$$ (7)
$$ \begin{eqnarray}\begin{array}{lll} & & \\ & & {U}_{\pm 1,0}^{(2)}(X)={D}_{\pm 1,0}{v}_{1}SD({v}_{1}X,\sqrt{2}/2),\end{array}\end{eqnarray}$$ (8)
where C± 1,0 and D± 1,0 are all constants related to the coefficients bij; CN and SD = SN/DN are Jacobi elliptic functions.In order to solve the PDEs, we take ϕ1(r,t) = f11r2 + f12, ϕ–1(r,t) = f21r2 + f22. Here, fij (i = 1,2, j = 1,2) are all functions of r and t; f11 and f21 represent the frequency chirps; f12 and f22 are phases, respectively. Now, solving PDEs (5) gives
Here bij,ci,μ1,μ2,ω,Ω and α are all constants. Thus, we obtain two families of exact vortex solutions to Eq. (2):
$$ \begin{eqnarray}\begin{array}{lll}{\varPsi }_{\pm 1,0} & = & {C}_{\pm 1,0}{v}_{0}\beta (r,t){{\rm{e}}}^{{\rm{i}}\alpha \theta +{\rm{i}}{\phi }_{\pm 1,0}(r,t)}CN({v}_{0}X\\ & & +K(\sqrt{2}/2),\sqrt{2}/2),\end{array}\end{eqnarray}$$ (9)
$$ \begin{eqnarray}\begin{array}{lll}{\varPsi }_{\pm 1,0} & = & {D}_{\pm 1,0}{v}_{1}\beta (r,t){{\rm{e}}}^{{\rm{i}}\alpha \theta +{\rm{i}}{\phi }_{\pm 1,0}(r,t)}SD({v}_{1}X,\sqrt{2}/2).\,\end{array}\end{eqnarray}$$ (10)
When imposing the bounded condition lim|r|→∞Ψ±1,0(r,θ,t) = 0, we have
where the natural number n is the quantum number, and
3. Structures of the vortex solutions
The structures of the vortex states for the exact vortex solutions (9) and (10) can be controlled by tuning the frequency of the trapped potential, which depend on the radial quantum number n and the topological charge α. The structures of the vortex states for the exact vortex solutions (9) and (10) are similar. We only take the vortex solutions (9) as an example to illustrate the structures of the vortex states in the following sections.
In real experiment, we take spin-1 BECs of alkali 87Rb atoms with N = 1.8 × 105 atoms. The radial frequency and the axial frequency are taken as ω⊥ = 2π × 18 Hz and ωz = 2π × 628 Hz.[30] For ε = 0, it can be found that the radial frequency is time-independent from Eq. (6). In Fig. 1, we show the density distributions and phase diagrams at t = 0 of the vortex solution (9) for different radial quantum numbers n and fixed topological charge α = 1. The first row in Fig. 1 corresponding to n = 1 is the lowest energy states. The second and third rows corresponding to n = 2 and n = 3, respectively, are two excited states. The last row in Fig. 1 shows the corresponding phase diagrams of the |Ψ± 1,0|2 for n = 1. It can be seen that the numbers of the ring structure of the vortex solutions increase one by one with different radial quantum numbers n. In Fig. 2, we show the properties of the vortex solution (9) by taking different values of topological charge α and the fixed radial quantum number n = 1. Figure 2 shows that the density profiles grow more and more localized with increasing topological charge α and the vortex expands outward with the increasing topological charge α.

Figure 1.The density distributions |Ψ± 1,0|2 and phase diagrams at t = 0 of the vortex solution (9) of the rotating spin-1 BECs for topological charge α = 1. (a1)–(a4) Evolution of the density distributions |Ψ1(x,y,0)|2 for different radial quantum numbers n and the phase diagram. (b1)–(b4) Evolution of the density distributions |Ψ–1(x,y,0)|2 for different radial quantum numbers n and the phase diagram. (c1)–(c4) Evolution of the density distributions |Ψ0(x,y,0)|2 for different radial quantum numbers n and the phase diagram. The parameters are taken as Ω = 0.5, μ1 = 8, μ2 = 4, ϵ = 0 and Ω0 = 1.

Figure 2.The density profiles |Ψ± 1,0|2 and phase diagrams at t = 0 for the vortex solution (9) of the rotating spin-1 BECs for different topological charge α and fixed quantum number n = 1. The blue lines denote the density profiles |Ψ± 1,0|2 for the topological charge α = 1, respectively. The green lines denote the density profiles |Ψ± 1,0|2 for the topological charge α = 2 and the red lines denote the density profiles |Ψ± 1,0|2 for the topological charge α = 3, respectively.
4. Dynamic stability analysis
It is important to generate the stable states in Bose–Einstein condensates and there are many methods to produce the stable states. For example, Saito and Ueda obtained the stable matter-wave bright soliton by tuning the strength of interactions via making use of Feshbach resonance.[31] Now we investigate the dynamical stability of the vortex solution (9) by performing some numerical simulations. Here we run the numerical simulations using the split-step Fourier transformation. The domain is composed of 400 × 400 grid points and the step sized of time integration is τ = 0.0001. The initial values are taken as Ψ±1,0(x,y,0). First, we consider the rotating spin-1 BECs with harmonic potential, i.e., ε = 0 in Eq. (6). Based on Eq. (2), the evolution of the density and phase diagram of the vortex solution (9) for fixed radial quantum numbers n and increasing topological charge α are illustrated in Fig. 3. Figure 4 shows the evolution of the density and phase diagram of the vortex solution (9) for the fixed topological charge α and different radial quantum numbers n. Figures 3 and 4 show that the vortex solution is stable only for n = α = 1 and the ring structures are destroyed when α ≥ 2.
![The evolution of density distributions |ψ± 1,0|2 and phase diagrams for the vortex solution (9) of the rotating spin-1 BECs for different values of topological charge α and fixed quantum number n = 1. The first column is stable vortex for the topological charge α = 1. The second and third columns are unstable vortex for the topological charge α = 2 and α = 3, respectively. The parameters are Ω = 0.5, μ1 = 8, μ2 = 4, ϵ = 0 and ω0 = 1. The domain is (x,y) ∈ [–5,5] × [–5,5] for all the cases.](/Images/icon/loading.gif)
Figure 3.The evolution of density distributions |ψ± 1,0|2 and phase diagrams for the vortex solution (9) of the rotating spin-1 BECs for different values of topological charge α and fixed quantum number n = 1. The first column is stable vortex for the topological charge α = 1. The second and third columns are unstable vortex for the topological charge α = 2 and α = 3, respectively. The parameters are Ω = 0.5, μ1 = 8, μ2 = 4, ϵ = 0 and ω0 = 1. The domain is (x,y) ∈ [–5,5] × [–5,5] for all the cases.

Figure 4.The evolution of density distributions |ψ± 1,0|2 and phase diagrams for the vortex solution (9) of the rotating spin-1 BECs for the fixed topological charge α and different quantum number n = 1. (a1)–(c3) Evolution of density distributions |ψ± 1,0|2 and phase diagrams for the quantum number n = 1. (d1)–(f3) Evolution of density distributions |ψ± 1,0|2 and phase diagrams for the quantum number n = 2. The parameters are Ω = 0.5, μ1 = 8, μ2 = 4, ϵ = 0 and ω0 = 1.
When the harmonic trap depends on the time, that is to say, ϵ ≠ 0 in Eq. (6), the time evolution of the density and phase diagram of the vortex solution (9) with two different values of topological charge α and radial quantum numbers n is concerned in Fig. 5. It is observed that the time-dependent frequency of the external potential affects the density distribution of the vortex states and the time of vortex decay from Figs. 3–5. The vortex states correspond to many high-energy collective excitation modes because of the dependence of interaction parameter and the frequency of the external potential on time and space. For the lowest mode with the topological charge α = 1 and quantum number n = 1, the vortex state is stable against a Gaussian noise. When the topological charge α = 2 or the quantum number n = 2, the vortex states have a tendency to decay due to the high-energy collective excitation. These characteristics can be observed from Figs. 3–5.

Figure 5.Time evolution of density distributions |ψ± 1,0|2 and phase diagrams for the vortex solution (9) of the rotating spin-1 BECs for different values of topological charge α and quantum number n. The first three lines show the evolution of density distributions |ψ± 1,0|2 and phase diagrams for n = α = 1. The three lines in the middle show the evolution of density distributions |ψ± 1,0|2 and phase diagrams for n = 1 and α = 2. The last three lines demonstrate the evolution of density distributions |ψ± 1,0|2 and phase diagrams for n = 2 and α = 1. The parameters are Ω = 0.5, μ1 = 8, μ2 = 4, ϵ = 0.1 and ω0 = 0.2.
5. Conclusions
In this paper, we focus on the quantized vortices in spinor BECs with time–space modulated interactions and take the alkali atoms 87Rb in the F = 1 hyperfine state as an example to show how to produce vortexes. Under the mean-field approximation, the dynamics of spinor condensates can be described by three-component GPE with an angular momentum rotational term. Two types of exact vortex solutions expressed by the Jacobi elliptic function are constructed by similarity transformation. The dynamical behaviors of the vortex solutions are analyzed by numerical stimulation. The results show that stable vortex states can be obtained by adjusting the frequency of the external potential and the spatiotemporally modulated interaction. We hope that the vortex state of spinor BECs with time–space modulated interactions can be realized in future experiments and help us to understand this behavior further.
Acknowledgment
Acknowledgment. The authors thank Professor W. M. Liu for his guidance.