Properties of a nematic spin vortex in an antiferromagnetic spin-1 Bose-Einstein condensate

A spin-1 condensate is described by the spinor field

with the three components representing the condensate amplitude in the spin levels $m=1,0,-1$, respectively, where $m$ is the quantum number associated with the $z$-component of spin. In weak fields the short-ranged contact interactions between atoms are rotationally invariant with a Hamiltonian density

Here the first term, with coupling constant $c_{0}$, describes the density dependent interactions, where $n\equiv\bm{\Psi}^{\dagger}\bm{\Psi}$ is the total density. The second term describes the spin-dependent interactions, where $c_{1}$ is the spin-dependent coupling constant, $\vec{F}\equiv\bm{\Psi}^{\dagger}\vec{f}\bm{\Psi}$ is the spin density, and $\vec{f}\equiv(\check{f}_{x},\check{f}_{y},\check{f}_{z})$ are the spin-1 matrices. In addition, we consider the presence of a (uniform) quadratic Zeeman shift. Taking the field to be along $z$ this is described by

In practice the coefficient $q$ is readily changed in experiments using microwave dressing (e.g. see Gerbier et al. (2006); Leslie et al. (2009)). We note that the uniform linear Zeeman term can be removed using a gauge rotation, and can be neglected. The spin properties also depend on the (conserved) $z$-magnetization $M_{z}=\int\differential{V}F_{z}$ of the system. Here we consider only $M_{z}=0$.

The case of antiferromagnetic interactions, where $c_{1}>0$, is realized with ²³Na atoms in their lowest hyperfine manifold. Here the condensate prefers to minimize the spin-density to reduce the spin-dependent interaction energy. For a condensate of uniform density $n_{b}$ and spin-density $\vec{F}=\vec{0}$ the spinor is in a polar state

where the real unit vector $\vec{d}=(d_{x},d_{y},d_{z})$ is the nematic director and $\theta$ is the global phase. Noting that $\bm{\Psi}_{\text{P}}$ is invariant under $\theta\to\theta+\pi$ and $\vec{d}\to-\vec{d}$, we see that $\vec{d}$ defines a preferred axis in spin space, but not a preferred direction along that axis. The ground state orientation of $\vec{d}$ is determined by the quadratic Zeeman energy, which is given by $\mathcal{H}_{\mathrm{QZ}}=qn_{b}(1-d_{z}^{2})$. Thus for $q>0$ the system maximizes $d_{z}^{2}$, by being in the EAP phase, i.e. $\bm{\Psi}_{\text{EAP}}=\sqrt{n_{b}}e^{i\theta}(0,1,0)^{T}$. The case of interest in this paper is the EPP phase for $q<0$ where $\vec{d}=(d_{x},d_{y},0)$. In addition to the global phase, the EPP ground state also breaks a U(1) symmetry in spin space. This can be seen from the director, which can be written as $\vec{d}=(\cos\varphi,\sin\varphi,0)$, i.e. $\bm{\Psi}_{\text{EPP}}=\sqrt{\frac{n_{b}}{2}}e^{i\theta}(-e^{-i\varphi},0,e^{i\varphi})^{T}$, where $\varphi$ is the angle the director takes with respect to the $x$-axis.

Note that in the EPP phase we have $\Psi_{0}=0$ and the system is effectively a two-component condensate. Indeed, several studies of the relevant EPP vortices we present in the next subsection have been performed in a two-component (or binary) condensate. For completeness we briefly mention the mapping of the spinor parameters onto an equivalent binary system. To do this we note that interaction Hamiltonian density $\mathcal{H}_{\mathrm{int}}$ may be expressed in the binary form

with intra-species coupling constant $g_{ii}=c_{0}+c_{1}$ (identical for both components) and interspecies coupling constant $g_{1,-1}=c_{0}-c_{1}$. For the antiferromagnetic case $g_{1,-1}<g_{ii}$ and the components are miscible Mineev (1974).

Here we consider a quasi-2D spinor gas where the order parameter manifold of the EPP phase permits vortices as point defects. To give the basic structure of such vortex states we write the wave function on the $xy$-plane with the vortex core taken at the origin. Sufficiently far from the core the general vortex state is of the form

where $\phi$ is the azimuthal angle in the $xy$-plane [i.e. from $\bm{\Psi}_{\text{EPP}}$, taking $\theta$ and $\varphi$ to vary as $\sigma_{\text{M}}\phi$ and $\sigma_{\text{S}}\phi$ in space, respectively]. Here $\sigma_{\text{M}}$ and $\sigma_{\text{S}}$ are the winding numbers associated with the mass and spin current around the vortex, respectively (e.g. see Kudo and Kawaguchi (2015)). Combining the circulations for the $m=\pm 1$ components we see that $\sigma_{\text{M}}\pm\sigma_{\text{S}}$ must both be integers for the field to be singled valued. There are 8 non-trivial cases with $|\sigma_{\text{M}}\pm\sigma_{\text{S}}|\leq 1$ which define the elementary vortices of interest.

The evolution of a spin-1 condensate is described by the Gross-Pitaevskii equation (GPE)

where the $\check{\mathcal{L}}$ is the nonlinear GPE operator

with $\alpha\in\{x,y,z\}$. Here $\mathds{1}$ denotes the identity matrix in spin space. The nonlinear terms $n$ and $\vec{F}$ are determined using $\bm{\Psi}$. Here we focus our attention on a quasi-2D system with spatial coordinates $\vec{\rho}=(x,y)$, where $c_{0}$ and $c_{1}$ are the quasi-2D coupling constants.

Stationary solutions of the form $\Psi_{m}(\vec{\rho},t)=\psi_{m}(\vec{\rho}\,)e^{-i(\mu+m\lambda)t/\hbar}$ satisfy the time-independent Gross-Pitaevskii equation

with nonlinear terms in $\check{\mathcal{L}}$ now evaluated with $\bm{\psi}$. Here $\mu$ and $\lambda$ are the chemical and magnetic potentials introduced as a Lagrange multipliers to conserve norm $N$ and $z$-magnetization $M_{z}$, respectively¹¹1These are defined by the thermodynamic relations $\mu=\left(\frac{\partial E}{\partial N}\right)_{\!S,V}$, and $\lambda=\left(\frac{\partial E}{\partial M_{z}}\right)_{\!S,V}$. For the bulk EPP phase at $T=0$ (with entropy $S=0$) $E=N(q+\frac{1}{2}c_{0}n_{b})+\frac{1}{2}c_{1}F_{z,b}M_{z}$ [from (3) and (2)], with $N=Vn_{b}$ and $M_{z}=VF_{z,b}$ for the appropriately dimensioned volume $V$.. For a uniform EPP phase condensate of bulk density $n_{b}$ and magnetization density $F_{z,b}$ we have

and $\lambda=c_{1}F_{z,b}$ (e.g. see Kawaguchi and Ueda (2012)).

An EPP phase vortex stationary state takes the form [generalizing Eq. (15)]:

where we have used the radial coordinate $\rho=\sqrt{x^{2}+y^{2}}$. Here we take $\bm{\chi}=[\chi_{1},\chi_{0},\chi_{-1}]^{\mathrm{T}}$ as real, and have explicitly imposed the circulation on each component using $\check{C}=\text{diag}\{e^{i(\sigma_{\text{M}}-\sigma_{\text{S}})\phi},e^{i\sigma_{\text{M}}\phi},e^{i(\sigma_{\text{M}}+\sigma_{\text{S}})\phi}\}$. Substituting (22) into Eq. (20), gives a radial equation for the $\bm{\chi}$:

with

Here

are the radial part of the kinetic energy operator, the effective chemical potential, and the effective linear and quadratic Zeeman shifts, respectively. We have also subtracted $q$ off the single particle energy [see Eq. (24)] for convenience, so that the adjusted chemical potential appearing in Eq. (26) is given by

and is independent of $q$ [cf. Eq.(21)]. Here we take $\mu_{b}=c_{0}n_{b}$ as a useful characteristic energy scale of the system, with associated healing length

In this work we consider a NSV with $(\sigma_{\text{M}},\sigma_{\text{S}})=(0,1)$. From our choice of $\mu_{b}$ being fixed [Eq. (29)], far away from the vortex core the system will approach the bulk value for number density , i.e. $n_{b}$. Additionally, we focus on the case $\lambda=0$ (and thus $\tilde{p}=0$), for which the bulk spin density is $F_{z,b}=0$. In this case the NSV is completely unmagnetized²²2The discrete symmetry associated with the operation $\check{f}_{z}\to-\check{f}_{z}$ is explicitly broken for the “biased” case with $\lambda\neq 0$, where a magnetization ($F_{z}\neq 0$) is generally preferred. In this work, we consider the “unbiased” case with $\lambda=0$. Even for $\lambda=0$, a spontaneous magnetization can occur locally, e.g. $F_{z}\neq 0$ in the core of a HQV. In this case, the spinor field is no longer represented by the real vector $\vec{d}$ like Eq. (10). with $\chi_{1}=-\chi_{-1}$, and the boundary conditions on the $\chi_{m}$ are:

where the prime denotes a derivative with respect to $\rho$. We solve for the stationary state solutions numerically using a gradient flow technique (e.g. see Bao and Cai (2013)) with a finite difference implementation of the derivative operators and boundary conditions. We use an equally spaced radial grid of $N_{\rho}$ points $\rho_{j}=(j-\frac{1}{2})\Delta\rho$ with $1\leq j\leq N_{\rho}$. We choose the point spacing $\Delta\rho$ to be much smaller than $\xi_{b}$ and typically use a maximum radius of $\rho_{\mathrm{max}}\approx 410\xi_{b}$, with $N_{\rho}=8192$. We choose to implement the outer boundary conditions in Eqs. (31) and (32) as $\chi^{\prime}_{m}(\rho_{\mathrm{max}})=0$.

Because the NSV stationary solutions are completely unmagnetized they are independent of the strength of the spin-dependent interaction. However, they depend on the quadratic Zeeman energy, and we find that there is a critical value

[or $q_{c}\approx-0.3414\mu$, see Eq. (29)], which separates the normal-core and polar-core forms of the vortex.

To derive the BdG equations we introduce a time-dependent fluctuation $\delta\bm{\psi}(\vec{\rho},t)=\bm{\psi}(\vec{\rho},t)-\bm{\psi}_{\text{V}}(\vec{\rho}\,)$ on the vortex stationary state $\bm{\psi}_{\text{V}}$ [see Eq. (22)]:

where $\beta_{\nu,\eta}$ are arbitrary (small) linearization amplitudes, and

Here we have introduced the radial quasiparticle amplitudes $\{\tilde{\bm{u}}_{\nu,\eta}(\rho),\tilde{\bm{v}}_{\nu,\eta}(\rho)\}$ and eigenvalues $E_{\nu,\eta}$, with $\eta$ being the quantum number associated with the $z$-component of angular momentum (relative to the condensate) and $\nu$ representing the remaining quantum numbers. Linearizing the time-dependent GPE (18) we obtain that the quasiparticle amplitudes satisfy the Bogoliubov-de Gennes (BdG) equation

where nonlinear terms in

are evaluated with $\bm{\psi}_{\text{V}}$, and

Because of the symmetry of the radial BdG equation for a solution (i) $(E,\eta,\tilde{\mathbf{u}},\tilde{\mathbf{v}})$ there are potentially three additional solutions which relate to the first as: (ii) $(-E,-\eta,\tilde{\mathbf{v}},\tilde{\mathbf{u}})$, (iii) $(E^{*},\eta,\tilde{\mathbf{u}}^{*},\tilde{\mathbf{v}}^{*})$, and (iv) $(-E^{*},-\eta,\tilde{\mathbf{v}}^{*},\tilde{\mathbf{u}}^{*})$ . If the eigenvalue is real, then $\tilde{\mathbf{u}}$ and $\tilde{\mathbf{v}}$ can also be taken real (see Ref. Takahashi et al. (2009)) and only (i) and (ii) are unique solutions. Furthermore quasiparticle amplitudes with real non-zero eigenvalues can be normalized to $\pm 1$ as

In the description of equilibrium condensates only positively normalized quasiparticles are considered to be physical. We note that the partner (ii) to a positively normed quasiparticle (i) has negative norm.

Here we are particularly interested in quasiparticles with complex eigenvalues where all the symmetries (i)-(iv) furnish unique solutions. If $\text{Im}(E)>0$ then the solution (i) is dynamically unstable (exponentially growing in time), and so is the partner solution (iv), while (ii) and (iii) are exponentially decaying solutions. Examining the effect of the unstable mode perturbation $\delta\bm{\psi}$ (36) we see modes (i) and (iv) are identical perturbations, so here we can choose to focus on solution (i).

The BdG results can reveal two types of instabilities for the NSV: (i) A dynamic instability revealed by a solution with a complex eigenvalue indicating that the respective eigenmode will exponentially grow with time. (ii) A Landau instability marked by a (positively normed) solution with a negative real eigenvalue, such that the system could reduce its energy if some dissipative mechanism allowed transfer of population into this state.

We have numerically calculated the BdG spectrum of the NSV over a wide parameter regime. We find that for $q<0$ the NSV only exhibits dynamic instabilities which occur in excitations with $\eta=\pm 1$, arising from modes that are localized in the vortex core⁵⁵5The centrifugal term in the $m=-\eta$ ($m=\eta$) component of the operator $\check{\mathcal{K}}_{\eta}^{+}$ ($\check{\mathcal{K}}_{\eta}^{-}$) vanishes for $\eta=\pm 1$, allowing the corresponding component of the $\bm{u}_{\nu,\eta}$ ($\bm{v}_{\nu,\eta}$) amplitude to develop amplitude in the NSV core. [see Fig. 3(d)]. The growth of these unstable modes causes the $m=\pm 1$ cores of the NSV to separate [see Fig. 3(e)], thus initiating the dissociation of the NSV into two HQVs (see Sec. II.2.3).

In Fig. 3(a) we present a stability phase diagram quantifying the dynamic instability of the NSV (i.e. showing the imaginary part of the eigenvalue of the dynamically unstable mode) as $q$ and $c_{1}$ vary. In general the imaginary part is always relatively small ($\lesssim 10^{-1}\mu_{b}$), so we expect the instability to manifest slowly in the system dynamics (see Sec. V). The dynamic instability is seen to depend on both $q$ and $c_{1}$. It is independent of $q$ for the normal-core NSV (i.e. $q<q_{c}$), but reduces with increasing $q$ in the polar-core regime [Fig. 3(b)]. The instability also decreases with increasing $c_{1}$ [Fig. 3(c)].

For $\eta=\pm 1$ we find that in addition to the unstable modes there are also zero-energy modes. These modes often reveal a broken symmetry in the system. These new zero modes are in addition to the usual $\eta=0$ zero energy mode, which can be written in terms of the condensate mode as $\{\bm{u}_{0,0},\bm{v}_{0,0}\}=\{\bm{\chi},-\bm{\chi}\}$ and is associated with the breaking of the gauge symmetry. The new zero energy modes with $\eta=\pm 1$ are associated with the breaking of translational symmetry in our uniform system by the presence of the NSV. These modes are localized in the core, similar to the unstable modes, but have a flat phase profile [see Fig. 4 and cf. Fig. 3(d)]. Whereas the unstable mode causes the $m=\pm 1$ component vortices to separate, the zero energy mode causes both to translate together. For a scalar condensate a similar zero mode emerges and in the case of a vortex line (i.e. finite $z$-extent), it is associated with the Kelvin-wave spectrum of helical modes that propagate along the vortex line.

We can develop an analytic expression for these zero energy modes that we compare to the numerical solution of the BdG equations. We restrict our attention to a normal-core NSV $(\sigma_{\text{M}},\sigma_{\text{S}})=(0,1)$, which has the form [see Eq. (22)] $\bm{\psi}_{V}=[-\chi_{-1}(\rho)e^{-i\phi},0,\chi_{-1}(\rho)e^{i\phi}]^{\mathrm{T}}$. Considering the vortex to be displaced by a small amount $\text{d}\vec{\rho}$ such that the change in the condensate wavefunction is $\delta\bm{\psi}=\text{d}\vec{\rho}\cdot\vec{\nabla}\bm{\psi}_{V}$, we obtain for the change in the nontrivial $m=\pm 1$ components

where

By inspecting the form of the BdG linearization [Eqs. (36)-(38)] we see that the vortex shift can be mapped to quasiparticle amplitudes with $\eta=\pm 1$. For definiteness we consider a zero energy mode of the BdG solution with $\eta=-1$, which we denote as $\{\tilde{\bm{u}},\tilde{\bm{v}}\}$, with amplitude $\beta$, such that

In comparison to Eq. (45) we observe $\tilde{u}_{\pm 1}\sim\Lambda_{\mp}$, and $\tilde{v}_{\pm 1}\sim\Lambda_{\pm}$, with the complex amplitude $\beta$ determining the vortex displacement. In Fig. 4 we show the numerically calculated BdG zero energy mode result and the functions $\Lambda_{\pm}$ (obtained from the vortex solution), showing they coincide.

It is challenging to numerically calculate the $\eta=\pm 1$ unstable and zero modes accurately. One issue is that our grids are of finite spatial extent $\rho_{\mathrm{max}}=410\xi_{b}$. This is problematic for the zero energy mode which decays rather slowly [noting that $\Lambda_{\pm}\sim\rho^{-1}$, see Eq. (46)]. Additionally, the 2D radial Laplacian is difficult to evaluate accurately with finite differences, particularly near $\rho=0$ Arsoski et al. (2015); Laliena and Campo (2018). In the BdG equations for $\eta=\pm 1$ the unstable and zero energy modes are particularly sensitive to these issues. We have found that second order finite difference schemes (including the improved schemes presented in Refs. Arsoski et al. (2015); Laliena and Campo (2018)) require an impractically large number of grid points to obtain accurate results. This motivated us to implement an 8th order finite difference scheme, which is the basis of the results we present. With this scheme small artefacts are still apparent such as the zero mode having a small imaginary part [e.g. see Fig. 3(b) and inset], causing it to couple to the unstable mode. These artefacts reduce as the range and point density of the numerical grid increase.