Vortex pairs in a spin-orbit coupled Bose-Einstein condensate

TL;DR

This paper explores how Rashba spin-orbit coupling influences vortex behavior in a two-component Bose-Einstein condensate, revealing deformed mass currents and unique vortex pair dynamics not seen in scalar condensates.

Contribution

It provides the first detailed analysis of vortex pairs in spin-orbit coupled BECs, highlighting altered vortex dynamics and stationary states due to spin-orbit effects.

Findings

Deformed mass current around vortex cores in the plane-wave phase.

Reduced vortex-antivortex pair velocity with spin-orbit coupling.

Existence of stationary vortex pair states.

Abstract

Static and dynamic properties of vortices in a two-component Bose-Einstein condensate with Rashba spin-orbit coupling are investigated. The mass current around a vortex core in the plane-wave phase is found to be deformed by the spin-orbit coupling, and this makes the dynamics of the vortex pairs quite different from those in a scalar Bose-Einstein condensate. The velocity of a vortex-antivortex pair is much smaller than that without spin-orbit coupling, and there exist stationary states. Two vortices with the same circulation move away from each other or unite to form a stationary state.

Equations

$$E_0\left[\mathbf{\Psi }\right]=\int d\mathbf{r}{\mathbf{\Psi }}^{†}\left(\frac{{\mathbf{p}}^2}{2m}-\frac{\hslash k_0}{m}\mathbf{p}\cdot {\mathbf{\sigma }}_{\perp }\right)\mathbf{\Psi },$$

$$E_{\mathrm{int}}\left[\mathbf{\Psi }\right]=\int d\mathbf{r}\left(\frac{{\text{g}}_0}{2}\sum _{j=1}^2|{\psi }_j{|}^4+{\text{g}}_{12}|{\psi }_1{|}^2|{\psi }_2{|}^2\right),$$

$$E[\mathbf{\Psi }]=E_0[\mathbf{\Psi }]+E_{\mathrm{int}}[\mathbf{\Psi }].$$

$$i\frac{\partial {\psi }_1}{\partial t}=-\frac{1}{2}{\nabla }^2{\psi }_1+i\kappa {\partial }_{-}{\psi }_2+\left(|{\psi }_1{|}^2+\gamma |{\psi }_2{|}^2\right){\psi }_1,$$

$$i\frac{\partial {\psi }_2}{\partial t}=-\frac{1}{2}{\nabla }^2{\psi }_2+i\kappa {\partial }_{+}{\psi }_1+\left(\gamma |{\psi }_1{|}^2+|{\psi }_2{|}^2\right){\psi }_2,$$

$$\bm{\Psi}(\bm{r})=\frac{1}{\sqrt{2}}\left(\begin{array}[]{c}e^{i\kappa x}\\ e^{i\kappa x}\end{array}\right),$$

$$v_{\xi }(\mathbf{r})$$

$$S_{\xi }(\mathbf{r})=\frac{1}{\rho (\mathbf{r})}{\mathbf{\Psi }(\mathbf{r})}^{†}{\sigma }_{\xi }\mathbf{\Psi }(\mathbf{r})(\xi =x,y,z)$$

$$\bm{\Psi}(\bm{r})=\frac{1}{\sqrt{2}}\left(\begin{array}[]{c}e^{i[\Phi(\bm{r})+\kappa x]}\\ e^{i[\Phi(\bm{r})+\kappa x]}\end{array}\right),$$

$${\psi }_1(x,y)\to {\psi }_2(x,-y),{\psi }_2(x,y)\to {\psi }_1(x,-y).$$

$$\bm{\Psi}(\bm{r})=\frac{1}{\sqrt{2}}\left(\begin{array}[]{c}e^{i[\Phi_{1}(\bm{r})+\kappa x]}\\ e^{i[\Phi_{2}(\bm{r})+\kappa x]}\end{array}\right).$$

$$E_0$$

$${\Phi }_1(\mathbf{r})={\tan }^{-1}\left(\lambda \frac{y-{\delta }_y/2}{x-{\delta }_x/2}\right),$$

$${\Phi }_2(\mathbf{r})={\tan }^{-1}\left(\lambda \frac{y+{\delta }_y/2}{x+{\delta }_x/2}\right),$$

$$E_0$$

$$\Phi (\mathbf{r})=\sum _{j=1}^2{n}_j{\tan }^{-1}\frac{y-y_j}{x-x_j},$$

$$\bm{\Psi}(\bm{r})=\left(\begin{array}[]{c}\sqrt{\rho_{1}(\bm{r})}e^{i[\Phi_{1}(\bm{r})+\kappa x]}\\ \sqrt{\rho_{2}(\bm{r})}e^{i[\Phi_{2}(\bm{r})+\kappa x]}\end{array}\right),$$

$${\Phi }_1(\mathbf{r})={\tan }^{-1}\left(\frac{y-\delta /2}{x}\right)-{\tan }^{-1}\left(\frac{y-\delta /2-d_y}{x}\right),$$

$${\Phi }_2(\mathbf{r})={\tan }^{-1}\left(\frac{y+\delta /2}{x}\right)-{\tan }^{-1}\left(\frac{y+\delta /2-d_y}{x}\right),$$

$${\rho }_1(\mathbf{r})$$

$${\rho }_2(\mathbf{r})$$

Full text

I Introduction

Topological excitations in superfluids originate from the intertwining between internal and external degrees of freedom in the order parameters. The simplest example is a quantized vortex in a scalar superfluid, in which the complex order parameter with the U(1) manifold winds around the vortex core, producing azimuthal superflow Onsager ; Feynman . For the order parameters with spin degrees of freedom, a rich variety of topological excitations are possible; these include skyrmions Leslie ; Choi , monopoles Ray , half-quantum vortices Seo , and knots Hall . Because of the close relationship between the spin and motional degrees of freedom in the topological excitations, we expect that their static and dynamic properties are significantly altered if there exists coupling between them, that is, if there exists spin-orbit coupling (SOC).

Recently, Bose-Einstein condensates (BECs) of ultracold atomic gases with SOC have been realized experimentally Lin ; Zhang ; Campbell ; Wu ; JLi ; in these experiments, the atomic spin or quasispin was coupled with the atomic momentum using Raman laser beams. Numerous theoretical studies have been performed to evaluate the static properties of topological excitations in spin-orbit (SO) coupled BECs, e.g., vortex arrays Wang , vortices in rotating systems XQXu ; Radic ; Zhou ; Liu , half-quantum vortices Sinha ; Lobanov , skyrmions Hu ; Kawakami ; CFLiu ; Chen ; YLi , topological spin textures Kawakami2 ; Xu ; Ruo ; Han ; Han2 , dipole-induced topological structures Deng ; Wilson ; Kato , and solitons with vortices Sakaguchi ; YCZhang . However, there have been only a few studies on their dynamics in SO-coupled BECs. The dynamics of a single quantized vortex in a harmonic trap was considered in Refs. Fetter ; Kasamatsu .

In this paper, we investigate the dynamics of a quantized vortex pair in a quasispin-$1/2$ BEC with Rashba SOC. When a singly quantized vortex is created in a uniform plane-wave state, the phase distribution around the vortex core is significantly altered by the SOC; this indicates that the mass current around the vortex is quite different from that without SOC and affects the dynamics of a vortex pair. As a result, a vortex-antivortex pair will be stationary or will travel much more slowly than one without SOC. The dynamics of a vortex-vortex pair with the same circulation are also quite different from those without SOC; the vortices move away from each other, or they approach each other and unite to form a stationary state.

This paper is organized as follows. The problem is formulated in Sec. II. The static properties of a single vortex are studied in Sec. III. The dynamics of a vortex-antivortex pair and those of a vortex-vortex pair with the same circulation are investigated in Secs. IV.1 and  IV.2, respectively. Conclusions are presented in Sec. V.

II Formulation of the problem

We consider a two-dimensional (2D) quasispin-1/2 BEC in a uniform system with Rashba SOC. Within the framework of mean-field theory, the system can be described by the order parameter $\bm{\Psi}(\bm{r})=[\psi_{1}(\bm{r}),\psi_{2}(\bm{r})]^{T}$, where $T$ denotes the transpose. The kinetic and SOC energies are given by

[TABLE]

where $m$ is the atomic mass, $k_{0}$ is the strength of the SOC, and $\bm{\sigma}_{\perp}=(\sigma_{x},\sigma_{y})$ are the $2\times 2$ Pauli matrices. The $s$-wave contact interaction energy is written as

[TABLE]

where $\textmd{g}_{0}$ and $\textmd{g}_{12}$ are the intra- and inter-component interaction coefficients, respectively. The total energy is given by

[TABLE]

In this paper, we consider an infinite system in which the atomic density $\bm{\Psi}^{\dagger}\bm{\Psi}$ far from vortices is a constant, $n_{0}$. In the following, we normalize the length, velocity, time, and energy by the healing length $\hbar/\sqrt{m\textmd{g}_{0}n_{0}}$, the sound velocity $\sqrt{\textmd{g}_{0}n_{0}/m}$, the characteristic time scale $\hbar/(\textmd{g}_{0}n_{0})$, and the chemical potential $\textmd{g}_{0}n_{0}$. The dimensionless coupled Gross-Pitaevskii (GP) equations, $i\partial\bm{\Psi}/\partial t=\delta E[\bm{\Psi}]/\delta\bm{\Psi}$, have the form

[TABLE]

where $\partial_{\pm}=\partial/\partial x\pm i\partial/\partial y$, $\kappa=\hbar k_{0}/\sqrt{m\textmd{g}_{0}n_{0}}$, and the ratio between the inter- and intra-component interactions is $\gamma=\textmd{g}_{12}/\textmd{g}_{0}$. The ground state is the plane-wave state for $\gamma<1$ and the stripe state for $\gamma>1$ Wang , which breaks the rotational symmetry of the system. In the following discussion, we will focus on the miscible case, $\gamma<1$, and the ground state is given by the plane-wave state,

[TABLE]

where the wave vector is chosen to be in the $x$ direction.

The velocity field is useful for understanding the dynamics of vortices. From the equation of continuity $\partial\rho/\partial t+\bm{\nabla}\cdot(\rho\bm{v})=0$ with atomic density $\rho=|\psi_{1}|^{2}+|\psi_{2}|^{2}$, we obtain the velocity field as

[TABLE]

where

[TABLE]

is the pseudospin density. The first term in Eq. (6) corresponds to the canonical part related to the superfluid velocity, and the second term corresponds to the gauge part induced by the SOC. The velocity field vanishes for the vortex-free ground state in Eq. (5), since the first and second terms in Eq. (6) cancel each other.

We numerically solve Eq. (4) by the pseudospectral method with the fourth-order Runge-Kutta scheme. In the imaginary-time propagation, on the left-hand side of Eq. (4), $i$ is replaced with $-1$. The numerical space is taken to be $400\times 400$, which is sufficiently large, and the effect of the periodic boundary condition can be neglected.

III Single vortex

We begin with a single vortex state, in which each component contains a singly quantized vortex. The initial state of the imaginary-time propagation is

[TABLE]

where $\Phi(\bm{r})=\tan^{-1}(y/x)$. After sufficiently long imaginary-time propagation, we obtain the stable stationary state, as shown in Fig. 1. Figures 1(a) and 1(b) show the density and phase distributions of the stationary state. We note that the phase defect in component 1 (2) is shifted in the $+y$ ($-y$) direction. We define the distance between the phase defects as $\delta$. The vortex core in each component is occupied by the other component. This structure can therefore be regarded as a pair of half-quantum vortices; nevertheless, we will refer to it as a “single vortex” in this paper. In the absence of SOC, such a pair of half-quantum vortices repel each other and cannot form a stationary state Eto . A similar structure is also found in a one-dimensional SOC system Kasamatsu . Figure 1(c) shows the spin distribution, and we can see a spin vortex near the origin. The dependence of the vortex shift $\delta$ on the SOC strength $\kappa$ is shown in Fig. 1(d), which implies $\delta\simeq 1/\kappa$.

Figure 2 shows the velocity field $\bm{v}(\bm{r})$ of the single-vortex state. The velocity field is greatly deformed by the SOC, compared with the rotationally symmetric velocity field without SOC. We note that the deformation of the velocity field extends over a wide range, and the upper region ($y\gtrsim 10$) exhibits a uniformly leftward velocity field, while the lower region ($y\lesssim 10$) is rightward. In these regions, $|\bm{v}|\lesssim 0.01$, which is much smaller than that without SOC, $|\bm{v}|=1/r$. This effect of SOC is also seen in Fig. 1(b), where the phase in the upper and lower regions is almost $\propto e^{i\kappa x}$, i.e., the $2\pi$ phase rotation around the vortex core is strongly compressed around the $x$-axis. The velocity field near the vortex core exhibits complicated structures containing multiple circulations, as shown in the right-hand panels in Fig. 2.

Due to the symmetry of the GP equation in Eq. (4), the single-vortex state with clockwise circulation can be obtained from that with counterclockwise circulation by the following transformation:

[TABLE]

By this transformation, the winding number of the vortex is inverted without changing the direction of the plane wave $e^{i\kappa x}$. Applying the transformation to the state shown in Fig. 1, we find that the vortex core in component 1 (2) shifts in the $+y$ ($-y$) direction also for the clockwise vortex. The velocity field and the pseudospin density are transformed as $v_{x}(x,y)\rightarrow v_{x}(x,-y)$, $v_{y}(x,y)\rightarrow-v_{y}(x,-y)$, $S_{x}(x,y)\rightarrow S_{x}(x,-y)$, and $S_{y}(x,y)\rightarrow-S_{y}(x,-y)$.

For a better understanding of the numerical results, we perform variational analysis. The variational wave function is

[TABLE]

Substitution of this wave function into Eq. (1) yields

[TABLE]

where $\chi=(\Phi_{1}+\Phi_{2})/2$, $\phi=\Phi_{1}-\Phi_{2}$, and the constant term is neglected. The first and second lines in Eq. (11) correspond to the kinetic and SOC energies, respectively. From the numerical results that the cores are shifted by $\delta$ and that the $2\pi$ phase rotation around the vortex core is compressed in the $y$ direction, the phases in Eq. (10) are assumed to be

[TABLE]

where $\lambda$ and $\bm{\delta}$ are variational parameters. We substitute these phases into Eq. (11) and integrate with respect to $\theta$. Because of the complicated structure near the vortex cores, we consider the region in which $r\gg 1$. The energy is

[TABLE]

which is minimized by $\delta_{x}=0$ and $\delta_{y}=1/\kappa$. Thus, the energy is lowered by the displacement of the vortex cores in the $y$ direction, and the displacement $\delta_{y}$ is estimated to be $1/\kappa$; this is in good agreement with the numerical results shown in Fig. 1(d). The energy in Eq. (13) decreases as $\lambda$ increases, and this accounts for the compressed $2\pi$ phase rotation. A better variational wave function will allow us to determine the value of $\lambda$. We note that the term $\propto\delta_{y}$ in Eq. (13) originates from the last term in the integrand of Eq. (11), which thus plays an important role in the vortex deformation due to the SOC.

IV Vortex pair

First, for clarity, we define the positions of the vortices and the distances between them, as shown in Fig. 3. The position of the phase defect of the $j$th vortex in component $i$ is denoted by $(x_{ij},y_{ij})$. As shown in Fig. 1(b), in each vortex, the cores in the two components are shifted by $\delta$ in the $y$ direction, and then $x_{1j}=x_{2j}$ and $y_{1j}-y_{2j}=\delta$. The position of the single vortex is defined by $(x_{j},y_{j})=((x_{1j}+x_{2j})/2,(y_{1j}+y_{2j})/2)$. For a vortex pair, the index $j$ is taken in such a way that $y_{1}>y_{2}$. The distance between the vortices is defined by $(d_{x},d_{y})=(x_{1}-x_{2},y_{1}-y_{2})$ and $d=(d_{x}^{2}+d_{y}^{2})^{1/2}$. The center of the vortex pair is defined by $(x_{c},y_{c})=((x_{1}+x_{2})/2,(y_{1}+y_{2})/2)$. The winding numbers of the first and second vortices are denoted by $\langle n_{1},n_{2}\rangle$. In the following subsections, we will consider the vortex pairs $\langle\pm 1,\mp 1\rangle$ and $\langle\pm 1,\pm 1\rangle$, which we call vortex-antivortex pairs and vortex-vortex pairs, respectively.

IV.1 Vortex-antivortex pair

In the absence of SOC, a vortex-antivortex pair travels at a constant velocity or is annihilated Aioi ; C. Rorai . A vortex-antivortex pair is stationary only in a trap potential Crasovan , and there is no stationary state in a uniform system.

In the presence of the SOC, our numerical results show that stable stationary vortex-antivortex pairs can be formed with a proper choice of the distance between vortices $d$; an example is shown in Fig. 4. We prepare the initial state in Eq. (8) with

[TABLE]

where, for this example, $n_{1}=\pm 1$, $n_{2}=\mp 1$, $x_{1}=x_{2}=0$, and $y_{1}=-y_{2}=d_{i}/2$ with $d_{i}$ being the initial distance between vortices. From this initial state, the imaginary-time propagation is performed sufficiently. The stationary state is always reached if the initial distance is $d_{i}\gtrsim 20$. It can be seen in Fig. 4 that the distance between vortices in the stationary state is $d_{y}\simeq 4.1$ for $\langle 1,-1\rangle$ and $d_{y}\simeq 6.5$ for $\langle-1,1\rangle$.

To understand the stabilization mechanism of the stationary vortex-antivortex pairs, we calculate the total energy using a model function given by

[TABLE]

with phases

[TABLE]

and densities

[TABLE]

where $\mathcal{N}(\bm{r})$ is the normalization factor to ensure $\rho_{1}(\bm{r})+\rho_{2}(\bm{r})=1$ and $\nu(x,y)=(x^{2}+y^{2})/(x^{2}+y^{2}+w^{2})$.

We set $\delta=1/\kappa$, and from the numerical results, the radius of the vortex $w$ is estimated to be $2w\simeq 1/\kappa$. Figure 5 shows the total energy as a function of $d_{y}$; this is obtained by substituting Eq. (15) into Eq. (3). It can be seen in Fig. 5 that local energy minima appear on either side of the global minimum and form the energy barriers, that stabilize the vortex-antivortex pair. We note that without SOC, there are no such barriers, as can be seen in Fig. 5(a) for $\kappa=0$. We also note that the barriers do not appear for uniform densities $\rho_{1}=\rho_{2}=1/2$, and the inhomogeneous densities of Eq. (17) are necessary for the barriers to form. Hence, we conclude that this is the combined effect of SOC and the nonlinear interaction.

We now turn our attention to the dynamics of the vortex-antivortex pair. Figure 6 shows the trajectories of the vortex cores, where the initial state is prepared as follows. We first prepare the state in Eq. (8) with the phase in Eq. (14), and then we allow the imaginary-time evolution for a short period (typically, $t\simeq 80$). From this state, the real-time evolution begins. Figures 6(a) and 6(b) show the dynamics of the vertically aligned vortex pair; the distance $d\simeq 10$ is larger than that of the stationary states shown in Fig. 4. The vortex-antivortex pair moves in the $-x$ and $+x$ directions at constant velocity with a fixed distance between vortices. These directions for the propagation agree with those for a scalar BEC. However, the velocities $v_{x}\simeq-0.006$ in Fig. 6(a) and $v_{x}\simeq 0.011$ in Fig. 6(b) are much slower than $v_{x}=1/d\simeq 0.1$, which is that seen in a scalar BEC for the same $d$. Figures 6(c) and 6(d) show the cases of oblique and horizontal alignments. The propagation directions of these vortex pairs are different from those in a scalar BEC. This can be understood by inspecting the velocity field shown in Fig. 2. For example, on the negative $x$-axis in the left-hand panel of Fig. 2(b), the velocity field is towards the lower right, which indicates that a vortex located on the left-hand side of the counterclockwise vortex will feel a mass current in this direction. Similarly, a vortex located on the right-hand side of the clockwise vortex will feel a mass current towards the lower right; this results in the dynamics shown in Fig. 6(d).

Figure 7 shows the velocity $v_{x}$ of the vertically aligned vortex pair (i.e., $d_{x}=0$) as a function of the vortex distance $d_{y}$, which is obtained by a method similar to that used to obtain Fig. 6. Such vortex pairs always travel in the $\pm x$ direction. The $d_{y}$ dependence of the velocity is quite different from that in a scalar BEC. For $\kappa=1$ (Fig. 7(a)), the velocity $v_{x}$ of the $\langle 1,-1\rangle$ pair (red circles) changes from negative to positive as $d_{y}$ increases, and $v_{x}=0$ at $d_{y}\simeq 5$, which corresponds to the stationary state seen in Fig. 4(a). The velocity $v_{x}$ of the $\langle-1,1\rangle$ pair (blue circles) also crosses the $v_{x}=0$ axis at $d_{y}\simeq 7$, which corresponds to the stationary state seen in Fig. 4(b). For $6\lesssim d_{y}\lesssim 8$, the velocity changes from negative to positive and from positive to negative as $d_{y}$ increases. For a relatively large distance between vortices ($d\gtrsim 10$), the propagation directions are the same as those of a scalar BEC, but the $d_{y}$ dependence of $v_{x}$ is weak; this can be understood from the fact that the velocity field is almost uniform far from the vortex core, as shown in Fig. 2. The velocity $|v_{x}|$ is always smaller than that in a scalar BEC for both $\langle 1,-1\rangle$ and $\langle-1,1\rangle$ pairs. There is no stable vortex-antivortex pair for small $d_{y}$; the vortices are unstable against pair annihilation.

The $\langle-1,1\rangle$ pair exhibits interesting dynamics when $\kappa$ is small. As shown in Fig. 7(b), there is no stable $\langle-1,1\rangle$ pair in the region $5.5\lesssim d_{y}\lesssim 9.6$. Figure 8 shows the dynamics of the $\langle-1,1\rangle$ pair with the initial distance $d_{y}=9.6$, where the initial state is prepared by the imaginary-time propagation for a short duration from the initial phase in Eq. (14) with $d_{y}>10$. There is no stable state for $d_{y}=9.6$ according to Fig. 7(b). As the vortex pair travels in the $-x$ direction, the distance $d_{y}$ decreases, and eventually the pair settles into a stable state with $d_{y}\simeq 5.7$; the excess energy is released from the vortex pair as density and spin waves.

IV.2 Vortex-vortex pair

In a scalar BEC, two quantized vortices with the same circulation move around each other. In contrast, the dynamics of vortex-vortex pairs with SOC are significantly different from those in a scalar BEC. Figure 9 shows the trajectories of vortices for the $\langle 1,1\rangle$ pair, where the initial state is prepared by the same method as in Fig. 6. When the initial positions are those shown in Fig. 9(a), they move away from each other. In the case of Fig. 9(b), the two vortices pass each other. The dynamics shown in Figs. 9(c) and 9(d) are more interesting. The two vortices approach each other and unite to form a stationary state, and the excess energy is released as waves. The resultant stationary state is stable and remains at rest, and the two vortices lie in a line perpendicular to the plane wave. In all cases, the center of the pair initially moves in the direction of $+x$. The transformation in Eq. (9) gives the dynamics of $\langle-1,-1\rangle$.

Figures 10(a) and 10(b) summarize the directions of the vortex motion when $d\simeq 10$ for the $\langle 1,-1\rangle$ and $\langle 1,1\rangle$ pairs, respectively. In Fig. 10(a), the motion of the center of the vortex-antivortex pair $(x_{c},y_{c})$ is indicated by the red arrows, and the relative position of the two vortices is nearly constant. In Fig. 10(b), the relative motion of the vortex-vortex pair is indicated by the blue arrow, and the center of the pair always shifts in the direction of the plane wave.

V Conclusions

We have investigated the behaviors of quantized vortices in quasispin-1/2 BECs with Rashba SO coupling in a uniform 2D system, where the atomic interactions satisfy the miscible condition and the ground state is the plane-wave state. We found that the static and dynamic properties of vortices are significantly different from those of a scalar BEC.

For a single vortex state, we found that the vortex cores in two components are shifted in the $\pm y$ directions by $\simeq 1/\kappa$ (Fig. 1). We also found that the phase distribution and velocity field around the vortex are greatly deformed compared with those of a scalar BEC (Figs. 1 and 2), which affects the dynamics of the vortex pairs. The vortex-antivortex pairs have stable stationary states at rest (Fig. 4), and this is in marked contrast to the vortex-antivortex pairs in a scalar BEC, which always travel. The stationary states can be explained by variational analysis (Fig. 5). Other than when in a stationary state, the vortex-antivortex pair travels at a velocity much slower than that for a scalar BEC with the same vortex distance. The dependence of the velocity and moving direction on the vortex location is also quite different from that in the case of a scalar BEC (Figs. 6 and 7). The vortex-vortex pair exhibits interesting dynamics: the vortices pass and move away from each other, or approach each other and combine into a stationary state (Fig. 9).

In experiments, the vortex states shown in this paper may be produced by the phase imprinting technique S. Burger ; J. Denschlag and the ensuing relaxation. The dynamics of vortices can be observed by the destructive imaging Neely or the real-time imaging Freilich . We hope that our numerical results presented in this paper can provide insight into a range of topics in the nonlinear dynamics of SO-coupled BECs.

Acknowledgements.

This work was supported by JSPS KAKENHI Grant Numbers JP16K05505, JP26400414, and JP25103007, by the NMFSEID under Grant No. 61127901, and by the Youth Innovation Promotion Association of CAS under Grant No. 2015334.