# Vortex pair dynamics in three-dimensional homogeneous dipolar superfluids

Srivatsa B. Prasad<EMAIL_ADDRESS>Joint Quantum Centre
Durham–Newcastle, School of Mathematics, Statistics and Physics, Newcastle
University, Newcastle upon Tyne, NE1 7RU, United Kingdom  Nick G. Parker
<EMAIL_ADDRESS>Joint Quantum Centre Durham–Newcastle, School of
Mathematics, Statistics and Physics, Newcastle University, Newcastle upon
Tyne, NE1 7RU, United Kingdom  Andrew W. Baggaley
<EMAIL_ADDRESS>Joint Quantum Centre Durham–Newcastle, School
of Mathematics, Statistics and Physics, Newcastle University, Newcastle upon
Tyne, NE1 7RU, United Kingdom

###### Abstract

The static and dynamic properties of vortices in dipolar Bose-Einstein
condensates (dBECs) can be considerably modified relative to their nondipolar
counterparts by the anisotropic and long-ranged nature of the dipole-dipole
interaction. Working in a uniform dBEC, we analyze the structure of single
vortices and the dynamics of vortex pairs, investigating the deviations from
the nondipolar paradigm. For a straight vortex line, we find that the induced
dipolar interaction potential is axially anisotropic when the dipole moments
have a nonzero projection orthogonal to the vortex line. This results in a
corresponding elongation of the vortex core along this projection as well as
an anisotropic superfluid phase and enhanced compressibility in the vicinity
of the vortex core. Consequently, the trajectories of like-signed vortex pairs
are described by a family of elliptical and oval-like curves rather than the
familiar circular orbits. Similarly for opposite-signed vortex pairs their
translation speeds along the binormal are found to be dipole interaction-
dependent. We expect that these findings will shed light on the underlying
mechanisms of many-vortex phenomena in dBECs such as quantum turbulence,
vortex reconnections, and vortex lattices.

##  I Introduction

The quantisation of circulation and discretization of vorticity is a striking
manifestation of superfluidity in interacting atomic Bose-Einstein condensates
(BECs). Whereas classical vortices are characterised by a continuous vorticity
field concentrated at the cores of individual vortices, the phase coherence of
superfluids results in a vanishing vorticity everywhere but along discrete
topological defect lines where the superfluid density must necessarily vanish.
While superfluid vortices undergo reconnections and annihilations, thereby
releasing compressible energy in the process in the form of phonons, the
overall vorticity of a non-dissipative superfluid is conserved [1]. In the
absence of such processes the dynamics of ensembles of quantum vortices in a
uniform superfluid background are well-described theoretically by the Biot-
Savart law [2, 3]. Furthermore, if all of the vortices are polarised
(anti-)parallel to a given axis and are restricted to move in the plane normal
to this axis, they can be modelled as incompressible point vortices [4, 5, 6].
Ensembles of quantum vortices have been the subject of extensive theoretical
and experimental studies. Notably, the ground states of rotating superfluids
are generally triangular vortex lattices, accompanied by a spectrum of lattice
vibrations in response to external perturbations [7, 8, 9, 10]. Meanwhile,
highly agitated superfluids exist in a state of quantum turbulence, a
disordered state of vortex lines and phonon excitations of the fluid [11, 12,
13]. In these macroscopic systems of vortices, the interaction between
vortices is key, for example, setting the equilibrium distance between
vortices in lattices and driving vortex-vortex reconnections in turbulent
systems, which are key to changing the topology of the flow and transfering
energy across scales [14, 15, 16]. This motivates the importance of
understanding the microscopic detail of the vortex-vortex interaction.

These phenomena manifest themselves in strikingly contrasting ways in dipolar
Bose-Einstein condensates (dBECs) which are composed of lanthanide atoms with
large, permanent dipole moments, such as chromium [17], dysprosium [18, 19],
erbium [20], and europium [21]. Here, assuming that all of the atoms are
uniformly polarised along an applied magnetic field, the superfluid properties
of the system are strongly modified by the dipole-dipole interaction (DDI)
between the atoms. Since the DDI is responsible for magnetostriction [22, 23]
– the energetic preference for dipolar atoms to mutually align themselves
along the dipolar axis – a directional dependence manifests itself in the
dynamics of a dBEC including a direction-dependent speed of sound and Landau
critical velocity [24, 25]. Thus, an intricate interplay exists between the
DDI and ‘environmental effects’ such as external trapping or rotation. Long
before the first experimental detection of vortex structures in dipolar BECs
(dBECs) occurred in 2022 [26, 27], a considerable body of research has been
conducted on single vortex ground states, vortex lattices and the dynamics of
vortex ensembles in dBECs subjected to harmonic confinement either in all
three dimensions or with a dominant confinement parallel to the axis of
vorticity, viz. the quasi-two-dimensional regime. Notably, the geometry of
vortex lattices has been predicted to be strongly dependent on the direction
of dipole alignment and strength of the dipole moments, with striped density
phases and the possibility of square lattices when the vorticity and the
applied magnetic field are not (anti-)parallel [28, 29, 30, 31, 32].
Additionally, in these trapped dBECs the cores of individual vortices have
been found to be perturbed by the dipolar interaction, leading to elongated
cores along the dipolar axis or density rippling associated with the
excitation of rotons [33, 34, 35]. Furthermore, the dynamical behaviour of
small ensembles of vortices in trapped dBECs has been demonstrated
theoretically to differ qualitatively from those of the equivalent nondipolar
systems [36, 37, 38, 39]. We also note that the competing influences of the
DDI and externally imposed confinement allows for vortices to be seeded in a
trapped dBEC through direct rotation of the magnetic dipole moments by
rotating the applied magnetic field [40], a method which was responsible for
the first unambiguous observation of vortices in dBECs in 2021 [26].

However, relatively little work has been conducted on the physics of vortices
in dipolar superfluids in domains that are (roughly) uniform over a
significantly sized volume in all three dimensions and which can be
effectively realised in optically generated ‘box traps’ [41, 42]. These would
serve as experimental platforms to investigate vortex-vortex interactions in a
relatively clean system where the interactions of the vortices with boundaries
and density inhomogeneities are minimised [13], allowing for the study of
quantum turbulence in a regime where the vast majority of theoretical
investigations – including in a dipolar setting [43] – have been carried out.
The possibility of future experiments in this direction motivates the
addressing of outstanding questions regarding the elementary properties of
vortices in this regime. For instance, it is not presently known whether the
anisotropic vortex cores and dipole-mediated dynamics previously predicted in
the presence of at least some external confinement are inherent in any dipolar
superfluid or if these phenomena are the result of instabilities arising from
the interplay between the DDI and the trapping. The work we present here aims
to resolve some of these questions and thereby shed light on the properties of
single vortex lines and pairs of vortex lines in a uniform $3$D dipolar BEC.

This article is structured as follows. In Section II we provide an overview of
the mean-field model employed to study vortices in a three-dimensional,
uniform dBEC while, in Sec. III, we search for stationary states of the system
containing a single vortex and examine the ways in which its properties
diverge from those of a vortex in a nondipolar BEC. These single-vortex
stationary solutions provide insight into the dynamics of vortex pairs, which
we first examine for the case where both vortices boast the same circulation
in Sec. IV before proceeding to analyse pairs of opposite circulation in Sec.
V. These findings are summarised in Sec. VI along with an outlook to future
lines of enquiry in this field.

##  II Formalism

In this work, we study the structure and dynamics of quantum vortices in
dipolar superfluids using the dipolar Gross-Pitaevskii equation (dGPE), which
adequately models the behaviour of interacting BECs in the mean-field
approximation. Let us consider a single atomic species of mass $m$ and
magnetic dipole moment $\mu_{\mathrm{d}}$ which is polarised uniformly by an
applied magnetic field parallel to $d$. For this dipolar BEC, the dGPE is
given by [44, 31, 45],

$i\hbar\partial_{t}\psi=\left\\{-\frac{\hbar^{2}}{2m}\nabla^{2}+g\left[n+3\varepsilon_{\mathrm{dd}}V_{\mathrm{dd}}(\mathbf{r})\circledast
n(\mathbf{r})\right]-\mu\right\\}\psi.$ (1)

Here, the two-body short-ranged interaction strength is
$g=4\pi\hbar^{2}a_{\mathrm{s}}/m$ with $a_{\mathrm{s}}$ the scattering length
of the atom-atom scattering potential, $n=|\psi|^{2}$ is the atomic density of
the condensate, and $\mu$ the chemical potential of the condensate which fixes
the normalisation of $\psi$ and satisfies the self-consistency relation,

$\mu=\int\mathrm{d}^{3}r\,\left\\{\frac{\hbar^{2}}{2m}|\nabla\psi|^{2}+gn^{2}+3g\varepsilon_{\mathrm{dd}}n(\mathbf{r})\left[V_{\mathrm{dd}}(\mathbf{r})\circledast
n(\mathbf{r})\right]\right\\}.$ (2)

The parameter
$\varepsilon_{\mathrm{dd}}=m\mu_{0}\mu^{2}/(12\pi\hbar^{2}a_{\mathrm{s}})$,
where $\mu_{0}$ is the permeability of free space, serves as an effective
dipolar-to-contact interaction ratio that encapsulates the degree to which the
dynamics of the superfluid are governed by the dipolar interaction.
Furthermore, the two-body dipole-dipole interaction is defined as [44]

$V_{\mathrm{dd}}(\mathbf{r})=\frac{1}{4\pi}\left[\frac{1-3(\hat{\mathbf{d}}\cdot\hat{\mathbf{r}})^{2}}{r^{3}}\right],$
(3)

and has a Fourier transform given by

$\widetilde{V}_{\mathrm{dd}}(\mathbf{q})=(\hat{\mathbf{d}}\cdot\hat{\mathbf{q}})^{2}-\frac{1}{3}.$
(4)

Note that while the global isotropy of the dGPE is broken in the presence of a
dipolar interaction, a rotational symmetry is still preserved about the
dipolar axis, $\hat{\mathbf{d}}$.

Since no external confining potential enters Eq. (1), its ground state
solution is a constant density $n_{0}$ and the corresponding chemical
potential is given by
$\mu_{\mathrm{g}}=gn_{0}\left\\{1+\varepsilon_{\mathrm{dd}}\lim_{q\rightarrow
0}\left[3(\hat{\mathbf{d}}\cdot\hat{\mathbf{q}})^{2}-1\right]\right\\}$. Thus,
in free space the solutions to Eq. (1) are unstable when
$\varepsilon_{\mathrm{dd}}>1$. While superfluids can still remain stable
beyond this threshold, as seen by the experimental realisation of self-
stabilised quantum droplets [46, 47] and supersolids [48, 49, 50] in dipolar
BECs in recent years, a theoretical description of this regime requires an
extension of the mean-field dGPE to incorporate quantum corrections to the
interaction [51, 52, 53]. In this work, we explore only the regime
$0<\varepsilon_{\mathrm{dd}}<1$ for which the mean-field description is
sufficient [45].

We choose to work in natural units where energy is expressed in units of
$\mu_{\mathrm{g}}$, length in units of the healing length
$\xi=\hbar/\sqrt{m\mu_{\mathrm{g}}}$, time in units of
$\hbar/\mu_{\mathrm{g}}$, and $\psi$ in units of $\sqrt{n_{0}}$ [54]. Using
these natural units, Eq. (1) can be written in dimensionless variables as [24,
43]

$i\partial_{t}\psi=\left[-\frac{1}{2}\nabla^{2}+V_{\mathrm{int}}(\mathbf{r})\circledast
n(\mathbf{r})-\mu\right]\psi,$ (5)

where

$\widetilde{V}_{\mathrm{int}}(\mathbf{q})=\frac{1+3\varepsilon_{\mathrm{dd}}\widetilde{V}_{\mathrm{dd}}(\mathbf{q})}{1+3\varepsilon_{\mathrm{dd}}\lim_{q\rightarrow
0}\widetilde{V}_{\mathrm{dd}}(\mathbf{q})}$ (6)

is the Fourier transform of the full scaled nonlocal two-body interaction.

It is also helpful to recast the dGPE in a hydrodynamic form by writing $\psi$
in the Madelung form $\psi=\sqrt{n}e^{iS}$ and interpreting the superfluid
phase $S$ as a scaled velocity potential for the superfluid velocity
$\mathbf{v}$ such that $\mathbf{v}=\nabla S$. Under these transformations, the
dGPE may be re-expressed as a pair of superfluid hydrodynamic equations [44],

$\displaystyle\partial_{t}n=\nabla\cdot(n\mathbf{v}),$ (7) $\displaystyle
m\partial_{t}\mathbf{v}=-\nabla\cdot\left\\{\frac{1}{2}m\mathbf{v}^{2}+V_{\mathrm{int}}(\mathbf{r})\circledast
n(\mathbf{r})-\mu\right\\},$ (8)

from which we shall extract considerable intuition about the expected form of
numerical solutions to the dGPE. For instance, the quantization of circulation
of $\psi$ [1], viz.

$\Gamma=\oint\,\mathrm{d}\mathbf{s}\cdot\mathbf{v}=2\pi q\,:\,q\in\mathbb{Z},$
(9)

is inherently a manifestation of the single-valued nature of $\psi$ while the
fact that the superfluid density vanishes at the core of a quantum vortex,
where $q\neq 0$, is a logical consequence of finding a solution of Eq. (7) for
$n(\mathbf{r})$ where $\mathbf{v}(\mathbf{r})$ is divergent. Hereafter, we
consider only vortices with a single quanta of circulation, i.e. $|q|=1$.

Our results are based on numerical solutions of the dGPE obtained by
propagation via the split-step pseudospectral method [55], in which the
spatial derivatives on the right-hand-side of the dGPE are computed using
Fourier-based methods and are propagated separately to the remaining terms.
Given that we wish to simulate a system where, far from the vortex cores, the
density is approximately constant, the choice of boundary conditions must
reflect this. This necessitates the use of distinct spatial grids and
corresponding spectral methods for computing the spatial derivatives of $\psi$
depending on the total circulation generated by vortices in the system. These
shall be described in greater detail in the relevant sections of this article.

##  III Single Vortex Structures for Orthogonal Dipole Alignments

In order to best understand vortex pairs, it is important to first examine the
properties of individual vortex lines. Surprisingly, vortex lines in uniform
$3$D dipolar BECs have received little attention when compared to those in
trapped $3$D systems or those that are uniform in quasi-$2$D limit, a regime
where the external confinement along the axis of vorticity is strong enough to
render the dynamics of the condensate effectively two-dimensional. We proceed
thus to probe the structure of a single vortex in a uniform dipolar BEC and
explore the effects of varying the dipolar interaction ratio,
$\varepsilon_{\mathrm{dd}}$, and the dipole alignment axis $\hat{\mathbf{d}}$.
Without loss of generality we take the vortex to be parallel to the $z$-axis
and focus on the two most disparate scenarios, viz. dipole alignments along
either the $z$\- or $x$-axes, which represent alignments parallel or
orthogonal to the vorticity, respectively. If the circulation of $\mathbf{v}$
around the boundary of a planar cross-section of a computational cell is
nonzero, the boundary conditions in this plane cannot be periodic if the
density profile of the condensate is to remain finite since a phase jump
between opposite sides of the cell is inevitable. Indeed, any attempt to solve
the dGPE in a periodic box with an initial state of nonzero circulation will
result in phase defects artificially emerging at the boundaries in order to
yield a vortex-neutral state inside the box. Thus, whenever there exists an
uneven number of vortices of positive and negative $\Gamma$, we employ Neumann
(reflecting) boundary conditions in the $x$-$y$ plane [56, 57], where the
normal derivative of $\psi$ vanishes at the boundaries, and periodic boundary
conditions along the $z$-axis. In what follows, the computational grid is of
size $(N_{x},N_{y},N_{z})=(513,513,2)$ with with equal grid spacings, $\Delta
x=\Delta y=\Delta z=25/256$, along each direction 111We use only 2 gridpoints
along $z$ as uniformity is assumed along this dimension.. To use
pseudospectral methods for computing the spatial derivatives, we also require
wavenumber grids of the same size as the position grids but with grid spacings
$\Delta k_{i}=\pi/[(N_{i}-1)\Delta r_{i}]$ in the $k_{x}$-$k_{y}$ plane and
$\Delta k_{z}=2\pi/(N_{z}\Delta z)$. Then, the spatial derivatives may be
computed using the discrete cosine transform [59] (DCT) 222More specifically,
to ensure that the same parity conditions are observed at both boundaries
along a given axis, we employ the DCT-I transform [82]. while the derivative
along the $z$-axis is computed via the discrete Fourier transform (DFT)/fast
Fourier transform [55] (FFT).

Let us first consider a dipole alignment parallel to the $z$-axis, for which
the condensate solution will be uniform along $z$. Writing
$\tilde{n}(\mathbf{k}_{\perp})=\mathcal{F}\left\\{n(\mathbf{r}_{\perp});\,\mathbf{r}_{\perp}\rightarrow\mathbf{k}_{\perp}\right\\}$
and evaluating $\left[V_{\mathrm{dd}}\circledast n\right](\mathbf{r})$
explicitly yields

$\displaystyle\left[V_{\mathrm{dd}}\circledast n\right](\mathbf{r})$
$\displaystyle=\frac{1}{(2\pi)^{2}}\int\mathrm{d}^{3}k\,e^{i\mathbf{k}_{\perp}\cdot\mathbf{r}_{\perp}}e^{ik_{z}z}\tilde{n}(\mathbf{k}_{\perp})\delta(k_{z})\left(k_{z}^{2}-\frac{1}{3}\right)$
$\displaystyle=-n(\mathbf{r})/3,$ (10)

implying that the only effect of the dipolar interaction in this regime is to
shift the value of the chemical potential depending on the value of
$\varepsilon_{\mathrm{dd}}$. Thus, neither the density or phase profile of a
straight vortex line are affected by the dipolar interaction if the atomic
dipole moments are polarised parallel to the vortex line. In the limit of an
infinite domain, the phase of such a vortex at the position
$(x,y)=(x_{v},y_{v})$ is thus given by the nondipolar limit
$S(x,y)=q\arctan(y-y_{v},x-x_{v}):\,n\in\mathbb{Z}$. On a grid with reflecting
boundaries in the $x$-$y$ plane, the phase is modified at the edges of the
grid to account for the boundary condition and it is this modified phase that
we initially impose upon the trial state $\psi$ prior to an attempt to obtain
a stationary state; a fuller discussion of the nondipolar vortex phase in a
reflecting box is provided in Appendix A.1. We obtain this vortex solution by
propagating the dGPE via the normalised gradient method, viz. evolving Eq. (5)
with $t\mapsto it$ and renormalizing $\psi$ at each timestep such that the
mean value of the density obeys $\langle n\rangle=\langle|\psi|^{2}\rangle=1$
as well as reapplying the vortex phase at each timestep. The corresponding
density profile is shown in Fig. 1 (a). As expected, the density is axially
isotropic due to the corresponding isotropy of the nondipolar GPE.
Furthermore, matching the properties of a $q=1$ vortex in a uniform nondipolar
BEC, $n(\boldsymbol{\rho})$ is proportional to $\rho^{2}$ as $\rho$, the
distance from the centre of the vortex core, tends to zero, and approaches
$\mu/\mu_{g}$ when $\rho\gg\xi$ [1].

Figure 1: Cross-sections of the density profiles, $n(\mathbf{r})$, in the
$x$-$y$ plane for a vortex of circulation $2\pi$ where the magnetic field is
polarised along the $x$-axis and $\varepsilon_{\mathrm{dd}}$ $=$ $0$ (a),
$0.3$ (b), $0.6$ (c), and $0.9$ (d). The density
$n(\boldsymbol{\rho})\sim\rho^{2}$ as $\rho\rightarrow 0$ and approaches a
constant value far from the vortex core. For nonzero
$\varepsilon_{\mathrm{dd}}$, the core is elongated along the dipolar axis.

By contrast, whenever the dipole polarisation has a nonzero projection
orthogonal to the vorticity, we must look to the superfluid hydrodynamic
equations to provide some insight about the form of the stationary solutions.
Let us consider a dipole moment located at the origin and polarised along the
$x$-axis. From the form of the DDI in real space, we observe that its long-
ranged interaction of this dipole moment with a second moment is dependent on
the angle between the dipole polarisation and the separation axis,
$\theta_{d}$. Since $V_{\mathrm{dd}}(\mathbf{r})\propto
1-3\cos^{2}\theta_{d}$, it is energetically preferable for the two dipole
moments to minimise the angle $\theta$, an effect referred to as
magnetostriction and which is responsible for the elongation of trapped
dipolar BEC clouds along the polarisation axis. From a superfluid hydrodynamic
perspective, a lowered dipolar interaction potential along the axis of the
polarisation naturally leads to higher values of the density along this axis
in stationary states of the Euler equation, Eq. (8)[61, 22, 23]. Via a similar
line of reasoning, the anisotropy of the density profile is diminished in the
far-field due to the $1/r^{3}$ dependence of the DDI.

In the particular case of a quantum vortex, its qualitative features can be
elucidated through a variation of the preceding arguments. Given that a vortex
line is characterised by a localised absence of dipolar superfluid atoms in
its core, it may be interpreted as a localised line of dipolar holes polarised
antiparallel to the atomic dipoles [33] in a manner analogous to the
counterparts of electrons in solid-state systems. Given that the generalised
form of the DDI for two dipole moments not necessarily parallel to each other
is of the form [44, 45]

$V^{(12)}_{\mathrm{dd}}(\mathbf{r})=\frac{3}{4\pi}\left[\frac{1-3(\hat{\mathbf{d}}_{1}\cdot\hat{r})(\hat{\mathbf{d}}_{2}\cdot\hat{r})}{r^{3}}\right],$
(11)

it is clear that the energetically preferential alignment for an atomic
dipole-vortex hole dipole pair is such that the angle between the applied
magnetic field and the mutual separation is $\pi/2$. Thus, the vortex core
becomes elongated along the direction of the applied magnetic field as the
virtual vortex dipole moments align preferentially along the magnetic field,
thereby expelling superfluid atoms in a larger domain along this axis relative
to the orthogonal axis. This, we note, is a novel form of the same
magnetostrictive effect that induces an elongation of the condensate density
envelope of a trapped dipolar BEC along the magnetic axis but is mediated by
virtual, rather than real, dipole moments [62].

While this core elongation has been discussed extensively in the literature
for both three-dimensional trapped [62, 27] and quasi-two-dimensional uniform
[36, 35] dipolar BECs, little attention has been paid to its consequence on
the phase profile of the vortex. Focussing on the continuity equation, Eq. (7)
rather than the Euler equation, Eq. (8) we argue that a logical consequence of
an anisotropy of the stationary state density profile must be a corresponding
anisotropy of the superfluid phase and velocity field. Thus, the normalised
gradient descent method is again used to find stationary states of the dGPE
with $\hat{\mathbf{d}}=\hat{x}$ where a single vortex is located at the centre
of the numerical grid but, while the density is renormalised throughout the
propagation, both the (normalised) density and phase are allowed to evolve
freely during the search for a stationary state of Eq. (5). The results of
these simulations are depicted as cross-sections of the density for
$\varepsilon_{\mathrm{dd}}=\\{0.3,0.6,0.9\\}$ in Fig. 1 (b), (c), and (d),
respectively. Here, it is clear that the vortex core becomes elongated along
the dipolar axis for nonzero $\varepsilon_{\mathrm{dd}}$ and that the degree
of elongation increases with the dipolar interaction strength. Note that no
global density rippling is present in the system, unlike in the quasi-two-
dimensional limit where rotons can be excited in the vicinity of the vortex
core [36, 35]; these rotonic ripples are a consequence of external confinement
which is absent in our analysis.

To gain further insight into these deviations from the nondipolar limit, we
also compute

$\delta
n(\mathbf{r})=n(\mathbf{r})-\frac{1}{2\pi}\int_{0}^{2\pi}\mathrm{d}\varphi\,n(r,\varphi),$
(12)

the discrepancy of the density at a given point from the angle-averaged
density at that radius, and $\delta S(\mathbf{r})=S(\mathbf{r})-S_{0}(x+iy)$,
the discrepancy of the phase from the incompressible solution, satisfying
Neumann boundary conditions, given by Eq. (15) in Appendix A.1. Furthermore,
we evaluate $\nabla\cdot\mathbf{v}$, a measure of the compressibility of the
superfluid velocity fields of the stationary solutions, and
$\Phi_{\mathrm{dd}}(\mathbf{r})=\langle
V_{\mathrm{dd}}(\mathbf{r})\rangle\equiv n(\mathbf{r})\circledast
V_{\mathrm{dd}}(\mathbf{r})$, the local expectation value of the dipolar
interaction potential. Each of these quantities are presented, below, in Fig.
2 for the values $\varepsilon_{\mathrm{dd}}=\\{0.3,0.6,0.9\\}$.

Figure 2: Cross-sections of the axial density anisotropy $\delta n$ (first
row), phase anisotropy $\delta S$ (second row), local dipolar interaction
energy $\Phi_{\mathrm{dd}}$ (third row), and compressibility
$\nabla\cdot\mathbf{v}$ (fourth row) as defined in the main text for
$\varepsilon_{\mathrm{dd}}$ $=$ $0.3$ (first column, a–d), $0.6$ (second
column, e–h), and $0.9$ (third column, i–l). An inverse correlation between
$\delta n$ and $\Phi_{\mathrm{dd}}$ is evident and the maximum magnitude of
both $\delta n$ and $\delta S$ increases with $\varepsilon_{\mathrm{dd}}$.
N.B. for the sake of clarity, we have set $\nabla\cdot\mathbf{v}$ – which
diverges at the phase singularity – to zero at $\mathbf{r}=0$.

From Fig. 2, it is immediately clear that the density profile anisotropy
$\delta n(\mathbf{r})$ in the first row is inversely correlated with
$\Phi(\mathbf{r})$ in the third row. In other words, the relative elongation
of the vortex cores along the $x$-axis is reflected in the larger values of
$\Phi$ along the $x$-axis as the superfluid atoms are ejected from regions in
which the local dipolar interaction energy is higher. As we had predicted
earlier, the plots of $\delta S$ in the second row demonstrate that the vortex
phase is modified by dipolar interaction. Interestingly, both $\delta n$ and
$\delta S$ take the form of lobes reminscent of $d$-wave orbitals, a symmetry
obeyed by the DDI, with the lobes corresponding to $\delta n$ being out of
phase compared to those of $\delta S$ by $\pi/4$. We also observe that the
fact that the maximum magnitude of $\delta n$ increases with
$\varepsilon_{\text{dd}}$ is accompanied by a similar behaviour exhibited by
$\delta S$, which we note is a consequence of $n$ and $S$ mutually satisfying
the continuity equation, Eq. (7). It is also evident that the phase
discrepancy manifests itself in the superfluid velocity field in the form of a
nonneglible compressibility in the vicinity of the vortex core in the fourth
row of Fig. 2. In particular, this behaviour exhibited by $\delta S$ suggests
that the velocity of test particles around the vortex is not uniform for a
given radius and that this non-uniformity is enhanced for more strongly
dipolar condensates. This, in conjunction with the anisotropy of $\Phi$, bears
intriguing consequences for the trajectory of a second vortex with respect to
this vortex, which we explore for the remainder of this article.

##  IV Same-signed Vortex Pair Dynamics

In Section III, it was observed that orthogonal magnetic dipole polarisations
had a considerable effect upon the density structure of a vortex and a small,
but non-negligible effect upon its phase structure. Just as the necessity of
satisfying the superfluid continuity equation in the stationary regime
provides an explanation for the interplay between the axial anisotropies of
the density and phase, we would expect divergent time-dependent dynamics of
small numbers of vortices with orthogonal polarisations. Conversely, as the
structure of a straight vortex line was found to be unmodified relative to its
nondipolar counterpart when the vorticity is parallel to the dipole
polarisation, it would be expected that the dynamics of ensembles of such
vortices would likewise be unaffected by a nonzero dipolar interaction
strength.

Figure 3: Trajectories of a vortex of circulation $\Gamma=2\pi$ initially
imprinted at the position $x(t=0)=0$, $y(t=0)=-s/2$, where $s$ $=$ $3.125$
(a), $6.25$ (b), or $12.5$ (c), and accompanied by a vortex of the same
circulation located initially at $x(t=0)=0$, $y(t=0)=s/2$. In each plot, we
have $\varepsilon_{\mathrm{dd}}$ $=$ $0.3$ (blue), $0.6$ (red), $0.9$ (black).
Due to the redistribution of energy at the start of each simulation the
vortices initially increase their mutual separation and we present the
trajectories of one closed orbit after a new equilibrium has been attained.
Here, an increasing degree of divergence from perfectly circular orbits for
larger (smaller) values of $\varepsilon_{\mathrm{dd}}$ ($s$) is apparent.

Thus, we employ the same numerical grid parameters and mixed Neumann-periodic
boundary conditions as in Sec. III and imprint the phases of two vortices,
both of circulation $\Gamma=2\pi$ at the positions $(x,y)=(0,\pm s/2)$,
continuously during a normalised gradient descent propagation of the dGPE
until convergence of the stationary solution is obtained. Subsequently, these
stationary solutions are propagated in real time with a timestep of $\Delta
t=0.001$; the positions of each vortex are tracked throughout the simulation
by identifying points on the grid around which the circulation is
approximately $2\pi q,\,q\in\mathbb{Z}$ [63, 64] and then employing a subgrid
least-squares interpolation method to further refine the estimated positions
of the vortices [65]. In the nondipolar, incompressible limit, the dynamics of
an ensemble of straight, (anti-)parallel vortex lines are governed by the
point-vortex model. More specifically, assuming a vorticity about the $z$-axis
and indexing each vortex by its position $\mathbf{r}_{i}$, the positions of
each vortex evolve according to the following coupled set of equations,

$\dot{\mathbf{r}}_{i}=\frac{-1}{2\pi}\sum_{j\neq
i}\,\Gamma_{j}\frac{\hat{z}\times(\mathbf{r}_{i}-\mathbf{r}_{j})}{|\mathbf{r}_{i}-\mathbf{r}_{j}|^{2}}.$
(13)

For a pair of same-signed vortices of circulation $\Gamma_{i}=\Gamma_{j}=2\pi$
and initial separation $r$, it is easy to show that the solution of Eq. (13)
is the vortex pair orbiting their centre of mass at a fixed radius, i.e. in a
circular orbit, and moving anticlockwise with an angular velocity
$\boldsymbol{\omega}=2\hat{z}/r^{2}$. Thus, the solution of the point-vortex
model for the two imprinted vortices in our dGPE simulations is that the two
vortices are locked in two circular orbits around the centre of mass
$\mathbf{r}=0$ at a radius $s/2$ and period $4\pi s^{2}$.

However, our numerical simulations demonstrate the existence of orbits whose
profiles are distinctly noncircular. In Fig. 3, we present plots of the same-
signed vortex trajectories with an initial displacement parameter $s=3.125\xi$
(a), $s=6.25\xi$ (b), and $s=12.5\xi$ (c), with
$\varepsilon_{\mathrm{dd}}\in\\{0,0.3,0.6,0.9\\}$ and $\hat{\mathbf{d}}$
aligned along the $x$-axis throughout the simulations. For the sake of
clarity, only one of the two vortex trajectories is plotted for each parameter
set; the phase shift of the second vortex with respect to the depicted vortex
is always approximately $\pi/2$. We also note that there is an initial
transient of phonons released from the vortices at $t=0$ as the solution
adjusts, and that this leads to a small initial increase in vortex separation;
in Fig. 3 we plot only the trajectories after the vortices have relaxed into
the equilibrium orbits that they occupy thereafter.

In Fig. 3 the qualitative shape of the trajectory of the vortex at fixed
$\varepsilon_{\text{dd}}$ appears to be largely independent of the initial
separation as the plots in each panel of the figure are very similar to each
other. However, striking differences are evident when focussing on any one of
the three panels, thereby fixing the initial vortex-vortex separation, and
instead varying $\varepsilon_{\text{dd}}$. When
$\varepsilon_{\mathrm{dd}}=0.3$, the trajectory of the vortex is almost
perfectly circular and is characterised by only a small ellipticity with a
semi-major axis along $\hat{y}$. A larger degree of ellipticity is seen for
$\varepsilon_{\mathrm{dd}}=0.6$, though, and for $\varepsilon_{\text{dd}}=0.9$
the vortex trajectories are clearly no longer elliptical but are instead
reminiscent of Cassini ovals [66].

We interpret these qualitative differences in the trajectories of the same-
signed vortices as being due to the anisotropy of the local dipolar
interaction energy, $\Phi_{\mathrm{int}}$, experienced by either vortex due to
the presence of the other, cf. Figs. 2 (c), (g), and (k). Whereas the dipolar
interaction energy for a dipolar atom located along the $x$-axis relative to a
vortex (and thus separated from the vortex line parallel to the dipolar axis)
is higher than that for an atom along the $y$-axis, the presence of a second
vortex line at a separation parallel to the dipolar axis is energetically
favourable compared to the orthogonal case because the second vortex is, as
argued in Sec. III, effectively a line of virtual dipole moments polarised
antiparallel to the atomic dipole moments. Thus, the higher kinetic energy
arising from the dipole moments being closer together when their separation is
parallel to the $x$-axis is counterbalanced by the lower dipole-dipole
interaction energy, thereby arising in elliptical and oval-like orbits with a
semi-major axis orthogonal to the dipole polarisation. Since the dipole-dipole
interaction energy decays with the distance from the vortex core, as evident
in the third row of Fig. 2, it is less favourable for orbits to be anisotropic
as the mean intervortex distance increases and this is evident in comparing
the trajectories in each panel of Fig. 3 for a given value of
$\varepsilon_{\mathrm{dd}}$.

It is important to note that in previous theoretical studies of vortex
dynamics in quasi-two-dimensional systems, like-signed vortex pairs were
observed to undergo a similar transition in orbit from circular, through
elliptical, to oval-like with increasing $\varepsilon_{\mathrm{dd}}$ when the
dipole polarisation had a nonzero projection in the $x$-$y$ plane [36]. The
same studies in the quasi-$2$D limit found a similar dependence of the dipolar
energy cost associated with the existence of a vortex pair on the orientation
of the intervortex separation axis, thereby explaining the deviations from the
nondipolar point-vortex paradigm. Similarly, a study of same-signed vortex
pairs in a harmonically trapped dBEC subject to dissipation predicted
different decay rates for the vortex trajectories depending on the initial
angle between the vortex separation and the dipolar axis [37]. Given our
results, we conclude that these discrepancies are an intrinsic feature of
same-signed vortex pairs in dBECs and do not arise necessarily from the
interplay of external confinement and the DDI.

##  V Vortex Dipole Dynamics

We now proceed to explore the dynamics of straight, antiparallel vortex lines,
i.e. vortex dipoles, in an orthogonally polarised dipolar BEC. Returning to
the point vortex dynamics of the incompressible limit and setting
$\Gamma_{1}=-\Gamma_{2}=2\pi$ one finds that $\dot{\mathbf{r}}_{12}$ is now
conserved and that
$\dot{\mathbf{R}}_{12}=(\hat{z}\times\hat{\mathbf{r}}_{12})/r_{12}$. Thus, the
vortex-antivortex pair undergoes translation along the binormal axis parallel
to $\hat{z}\times\hat{\mathbf{r}}_{12}$ at the constant speed $v=1/r_{12}$.
For instance, for a pair of vortices with circulations $\pm 2\pi$ at the
initial positions $(x,y)=[x(0),y(0)\pm s/2]$, the trajectories of the vortices
as predicted by the point-vortex model are

$[x_{i}(t),y_{i}(t)]=\left[x(0)+\frac{t}{s},\,y(0)\pm\frac{s}{2}\right].$ (14)

In the preceding section, however, it was shown that the strongly enhanced
compressibility of the superfluid flow and, more crucially, the vortex-vortex
dipolar interaction cause deviations from the point vortex model that are
observable even when the vortex cores are not overlapping. Such deviations
would then be detectable in future experiments in box-trapped dipolar BECs
[41, 42] where vortex-antivortex pairs are generated by dragging a potential
barrier through the system. Hence, we imprint a pair of vortices of
circulation $\Gamma=\pm 2\pi$ at the positions $(x,y)=(0,\pm s/2)$ with
$s=6.25\xi$ on a grid with periodic boundary conditions in all three
directions since the neutral overall circulation of the system allows for a
stationary state without the formation of phase defects at the boundaries.
Specifically, in this case the dGPE is solved on a
$(N_{x},N_{y},N_{z})=(512,512,2)$ grid with equal grid spacings in each
direction, given by $\Delta x=\Delta y=\Delta z=25/256$, and the spatial
derivatives are computed via FFTs; the wavenumber grids have spacings $\Delta
k_{i}=2\pi/(N_{i}\Delta r_{i})$. As in the previous section, the true phase of
the vortex-antivortex pair is not simply a composition of the phases of the
two vortices in the limit of an infinite domain,
$\arctan(y,x-s/2)-\arctan(y,x+s/2)$, as we need to account for the effect of
the boundaries; the modified phase that we specify in the initial conditions
is provided in Appendix A.2. Similar to our analysis of the same-signed vortex
pair, we choose to vary $\varepsilon_{\mathrm{dd}}$ from $0$ to $0.9$ and
allow $\hat{\mathbf{d}}$ to be parallel to $\hat{x}$, $\hat{y}$ and $\hat{z}$.
The results of these simulations are presented as a plot of
$v_{x}=\dot{x}_{i}$ as a function of the full range of
$\varepsilon_{\mathrm{dd}}$ in Fig. 4.

Figure 4: Velocities along the $x$-axis of a vortex-antivortex pair of
circulation $\Gamma=\pm 2\pi$ with an initial separation of $d=$ $3.125\xi$
(a), $6.25\xi$ (b), and $12.5\xi$ (c). The vortices initially approach each
other slightly before relaxing into a new equilibrium separation due to the
effects of the transient energy and thus we ignore the first $25$ units of
time. We scale the velocities after this time as
$v_{x}^{\prime}=v_{x}(\varepsilon_{\mathrm{dd}})/v_{x}((\varepsilon_{\mathrm{dd}}=0)$
and multiply them by the scaled vortex-antivortex separation along the
$y$-axis, $\langle\Delta y\rangle^{\prime}=\langle\Delta
y(\varepsilon_{\mathrm{dd}})\rangle/\langle\Delta
y(\varepsilon_{\mathrm{dd}}=0)\rangle$.

We have found that for $s\gg\xi$, the vortex-antivortex trajectories are still
straight lines along the binormal axis except at the earliest times when $s$
adjusts slightly due to the initial transient rearrangement of energy in the
simulations. Thus, the plots of the translational velocity in Fig. 4 are
scaled by the separation of the two vortices, $\langle\Delta y\rangle$, after
an equilibrium has been attained; we note that, in the point-vortex limit, the
quantity $v_{x}\langle\Delta y\rangle$ should be equal to unity. Instead, we
find that a strong dipolar dependence exists. For the full range of
$\varepsilon_{\mathrm{dd}}$ that we probe, $v_{x}\langle\Delta y\rangle$
increases monotonically with $\varepsilon_{\mathrm{dd}}$ when the dipole
moments are polarised along the $y$-axis, viz. parallel to the intervortex
separation.

The corresponding behaviour of the vortices when the dipole moments are
polarised along the $x$-axis, i.e. translational axis, is less simple. In Fig.
4 (a), where the vortex separation is initially $3.125$ healing lengths, the
translational velocity initially increases with $\varepsilon_{\mathrm{dd}}$
and undergoes a maximum (which exceeds the corresponding value of
$v_{x}\langle\Delta y\rangle$ with the dipole moments polarised along the
$y$-axis) before decreasing for larger $\varepsilon_{\mathrm{dd}}$. This is in
contrast to the behaviour observed in Figs. 4 (b) and (c), where the initial
separations are $6.25$ and $12.5$ healing lengths, respectively, which we
attribute to the effects of a slight overlap of the vortex cores in Fig. 4
(a). Instead, for the larger separations the scaled translational velocities
decrease for larger $\varepsilon_{\mathrm{dd}}$ before attaining a minimum and
increasing slightly. The dipolar dependence of the translational velocities
are likely due to an interplay between the dipolar correction to the
superfluid velocity induced by a single vortex, as elucidated in Sec. III, and
the dipole-dipole interaction between the vortices themselves. We note that
our solution scheme approximately preserves the total energy of the system
and, due to the point-vortex prediction of constant intervortex separation
being upheld by our simulations after the initial mutual approach, the vortex-
vortex dipolar energy is also conserved after this initial time. Thus, in
comparison to the trajectories of same-signed pairs in Sec. IV, a simple
explanation of the vortex-antivortex translation velocity discrepancy in terms
of conserving the dipolar interaction energy is less forthcoming. Further
investigations of this regime are necessary to clarify these questions.

It is also pertinent to note that previous theoretical studies of vortex-
antivortex pairs in dipolar superfluids have noted a dependence of the
critical velocity of vortex-antivortex formation when dragging an obstacle
through a dipolar superfluid on the angle between the obstacle velocity and
the dipole polarisation, an effect attributed to the anisotropy of the roton
dispersion at finite wavenumber resulting from the obstacle’s presence [24,
67]. Other studies in the quasi-two-dimensional regime have suggested that the
presence of a transverse dipole-dipole interaction can even suppress the
annihilation of a vortex-antivortex pair initially separated closely enough
that the equivalent pair in a nondipolar condensate would be annihilated [36].
However, to the best of our knowledge, no prior investigations have uncovered
a relation between the vortex-antivortex propagation velocity – as opposed to
the critical velocity of an obstacle dragged through a superfluid that is
required to nucleate such a pair – and the direction of the dipole
polarisation.

##  VI Conclusion

In this work, we have demonstrated that the anisotropy of the magnetic dipolar
interaction can fundamentally alter the properties of both the stationary
states and dynamics of quantum vortices in a superfluid as compared to the
nondipolar paradigm. Whereas prior studies had focussed on scenarios that were
then more amenable to experimental probes, such as harmonically trapped
systems and systems strongly confined along the axis of vorticity, we have
shown that the qualitative deviations that these investigations observed are
evident even in a wholly uniform, three-dimensional system. Firstly, the
stationary states of a single, straight vortex line are shown to exhibit axial
anisotropy in both the density and phase profiles. This is seen to be due to
the effective dipolar interaction near the vortex core between the dipolar
atoms of the bulk and the virtual oppositely-polarised dipolar holes of the
vortex line. While a modification of the phase had been noted in a previous
study of single vortices in harmonically trapped dBECs [62], this phenomenon
had not previously been predicted in the uniform limit. For systems of pairs
of vortices, both the like-signed and opposite-signed cases exhibit deviations
from the nondipolar regime. In the case of both vortices having the same
circulation, we have found that their orbits are described by a family of
curves from circles to ovals depending on the strength of the dipole-dipole
interaction and their initial separation, a phenomenon we attribute to the
effective dipolar interaction between the vortex lines. When the two vortices
comprise a vortex dipole, their trajectories remain a translation at constant
velocity but these velocities are found to be dipole-dependent; for
intervortex separations larger than $\sim 4$ healing lengths, the relationship
between the translational velocities, dipole orientation and dipolar
interaction strength is found to be essentially universal with vortex-
antivortex pairs in a dBEC polarised parallel to the mutual separation moving
at a higher velocity than those in a dBEC polarised along the translational
axis.

While the direct relevance of our study is for straight vortex lines in a
three-dimensional uniform dBEC, it is pertinent to note that the qualitative
ramifications of our findings extend beyond this regime to any system with
phase defects in a uniform background. For instance, we expect that the
nonneglible dependence of the translational velocity of a vortex-antivortex
dipole on the direction and strength of the magnetic dipoles is a general
feature of soliton-like excitations in dipolar BECs, regardless of the
effective dimension of the system, such as Jones-Roberts solitons, vortex
rings and rarefaction pulses [68, 69, 67] 333A. Villois and D. Proment,
private communication. Indeed, to the best of our knowledge, the properties of
vortex rings have not been studied systematically in a dipolar BEC.
Furthermore, given our demonstration that the phase of a vortex is modified by
the dipolar interaction when the dipole orientation is not fully parallel to
the vortex line, it must be assumed that the superfluid velocity field induced
by any two- or three-dimensional phase defect is similarly modified.

These results also raise further questions about the deeper nature of a vortex
in a dipolar BEC. For instance, the presence of lobes of compressibility about
the vortex core when the dipole moments are orthogonal to the vortex warrants
a semi-analytical treatment of the radial and angular dependence of the
discrepancy of the phase profile relative to the nondipolar paradigm. We also
believe that the strong divergence of the same-signed vortex pair trajectories
from the familiar point-vortex predictions warrants an investigation into
extending the point-vortex model to account for the dipolar interaction
heuristically, thereby allowing for computationally inexpensive predictions of
the trajectories of larger ensembles of vortices; we note that such modified
point-vortex models have been employed in the studies of BECs subjected to
effects such as external trapping, dissipation or the presence of multiple
species [71, 72, 73]. Additionally, we note that while the effects of
transversely polarised magnetic dipoles on the Kelvin wave excitations of
single vortex lines have been studied in trapped dBECs [33, 74, 75], the phase
profile modification observed in our study would bear ramifications for the
Kelvin wave spectrum of straight vortex lines in uniform dBECs as well. It is
also possible that the dipolar dependence of the vortex-antivortex
translational velocities bears consequences for reconnection dynamics in
three-dimensional dBECs where the translational motion of the vortices induces
growth of Kelvin waves; our results would suggest the existence of a
directional dependence in quantities such as the reconnection timescale.
Finally, while the Tkachenko oscillations of a vortex lattice in a
transversely polarised dipolar BEC have been studied theoretically in
harmonically trapped systems [38], our predictions of a large degree of
modification to the dynamics of same-signed vortex pairs suggests that the
Tkachenko oscillation spectrum would be modified in the uniform background
case as well due to the influence of the vortex-vortex dipolar interaction on
the lattice dynamics.

###### Acknowledgements.

This work was funded by Grant No. RPG-2021-108 from the Leverhulme Trust. S.
B. P. acknowledges fruitful discussions with Thomas Billam, Brendan Mulkerin
and Alberto Villois. This research made use of the Rocket High Performance
Computing service at Newcastle University.

##  Appendix A Vortex Phase Profiles in the Nondipolar Limit

###  A.1 Neumann Boundary Conditions in the x-y Plane

In a $3$D domain with Neumann boundary conditions in the $x$-$y$ plane and
periodicity in the $z$-direction, the infinite-domain expression for the phase
of a nondipolar quantum vortex at $(x,y)=(x_{v},y_{v})$,
$S(x,y)=n\arctan(y-y_{v},x-x_{v})$, does not respect the boundary conditions.
Previous studies of vortex dynamics in systems with these boundary conditions
have avoided any resulting issues by imprinting additional image vortices
outside the computational domain in such a way that the velocity projections
normal to the grid boundaries are roughly zero [56, 76].

In this work we adopt an alternative approach. In a sufficiently large but
finite domain, we assume incompressibility and take the phase to be the
harmonic conjugate of the $2$D Poissonian Green function in this domain which,
in turn, represents the superfluid equivalent of the stream function. To wit,
in a rectangular grid $x\in[0,L_{x}],\,y\in[0,L_{y}]$ with reflecting
boundaries, the appropriate incompressible phase for a vortex at
$(x,y)=(x_{v},y_{v})$ is given by [77, 78]

$S_{0}(w)=\frac{\Gamma}{2\pi}\operatorname{Arg}\left[\frac{\sigma(w-w_{v};\omega)\sigma(w+w_{v};\omega)}{\sigma(w-\overline{w}_{v};\omega)\sigma(w+\overline{w}_{v};\omega)}\right],$
(15)

where $w=(x+L_{x}/2)+i(y+L_{x}/2)$, $w_{v}=(x_{v}+L_{x}/2)+i(y_{v}+L_{x}/2)$
$\omega=(Lx,iL_{y})$ represents the pair of half-periods of the Weierstrass
$\sigma$ function,

$\displaystyle\sigma[w;(\omega_{1},\omega_{2})]=\prod_{\begin{subarray}{c}m,n=-\infty\\\
(m,n)\neq(0,0)\end{subarray}}^{\infty}\left[\left(1-\frac{w}{\omega_{mn}}\right)\left(\frac{w}{\omega_{mn}}+\frac{w^{2}}{2\omega_{mn}^{2}}\right)\right],$
$\displaystyle\omega_{mn}=m\omega_{1}+n\omega_{2}.$ (16)

###  A.2 Periodic Boundary Conditions in the x-y Plane

By contrast, the imposition of periodic boundary conditions in the $x$-$y$
plane precludes the existence of solutions to the dGPE with a nonzero
circulation of $\mathbf{v}$, the superfluid velocity. Thus, we can only
consider ensembles of vortex-antivortex pairs and the fundamental building
block of their incompressible velocity profiles – and their corresponding
superfluid phases – is that of a single vortex-antivortex pair for the
specified computational grid.

For a single vortex-antivortex dipole where the vortex with circulation $\pm
2\pi$ is located at the position $(x_{\pm},y_{\pm})$, the effects of the
periodic domain boundaries on the corresponding superfluid can be evaluated by
considering not only the given pair but an infinite series of image copies
comprising a tiling of the $x$-$y$ plane. The solution for the phase is thus
the summation of the solution in the infinite-cell limit,
$S_{\infty}(\mathbf{r})=\arctan(y-y_{+},x-x_{+})-\arctan(y-y_{-},x-x_{-})$
with $(x_{\pm},y_{\pm})$ replaced by the coordinates of a given periodic copy
of the original vortex pair. For a grid with $x\in[0,L_{x})$ and
$y\in[0,L_{y})$, the solution for the superfluid phase has been provided in
the literature by [79],

$\displaystyle S(x,y)=\sum_{p=-\infty}^{\infty}\biggl{\\{}$
$\displaystyle\mp\operatorname{atan}\left[\tanh\left(\frac{\pi
Y_{\pm}}{L_{y}}+p\pi\right)\tan\left(\frac{\pi
X_{\pm}}{L_{y}}-\frac{\pi}{2}\right)\right]$
$\displaystyle+\pi\left[\Theta(X_{+})-\Theta(X_{-})\right]\biggr{\\}}-\frac{2\pi(x_{+}-x_{-})y}{L_{x}L_{y}},$
(17)

where $X_{\pm}=x-x_{\pm}$ and $Y_{\pm}=y-y_{\pm}$. Taking the gradient of this
expression yields an equivalent quantity, up to a constant shift, to the well-
established result for the superfluid velocity of a vortex-antivortex dipole
in a periodic grid [80, 81]. Following the procedure of Ref. [79], we replace
the infinite series $\sum_{p=-\infty}^{\infty}$ with a finite series
$\sum_{p=-P}^{P}$ for the sake of computational tractability; we have found
that no significant suppression of initial ‘transient’ in the real time
propagation of the dGPE occurs if $P$ is increased beyond $O(10^{1})$. Thus,
in our work, we have set $P=11$.

## References

  * Fetter [2009] A. L. Fetter, Rotating Trapped Bose-Einstein Condensates, Reviews of Modern Physics 81, 647 (2009).
  * Schwarz [1985] K. W. Schwarz, Three-Dimensional Vortex Dynamics in Superfluid 4He: Line-Line and Line-Boundary Interactions, Physical Review B 31, 5782 (1985).
  * Schwarz [1988] K. W. Schwarz, Three-Dimensional Vortex Dynamics in Superfluid 4He: Homogeneous Superfluid Turbulence, Physical Review B 38, 2398 (1988).
  * Lin [1941a] C.-C. Lin, On the Motion of Vortices in Two Dimensions-I. Existence of the Kirchhoff-Routh Function, Proceedings of the National Academy of Sciences 27, 570 (1941a).
  * Lin [1941b] C.-C. Lin, On the Motion of Vortices in Two Dimensions-II Some Further Investigations on the Kirchhoff-Routh Function, Proceedings of the National Academy of Sciences 27, 575 (1941b).
  * Onsager [1949] L. Onsager, Statistical Hydrodynamics, Il Nuovo Cimento 6, 279 (1949).
  * Ho [2001] T.-L. Ho, Bose-Einstein Condensates with Large Number of Vortices, Physical Review Letters 87, 060403 (2001).
  * Cooper  _et al._ [2001] N. R. Cooper, N. K. Wilkin, and J. M. F. Gunn, Quantum Phases of Vortices in Rotating Bose-Einstein Condensates, Physical Review Letters 87, 120405 (2001).
  * Baym [2003] G. Baym, Tkachenko Modes of Vortex Lattices in Rapidly Rotating Bose-Einstein Condensates, Physical Review Letters 91, 110402 (2003).
  * Watanabe  _et al._ [2004] G. Watanabe, G. Baym, and C. J. Pethick, Landau Levels and the Thomas-Fermi Structure of Rapidly Rotating Bose-Einstein Condensates, Physical Review Letters 93, 190401 (2004).
  * Henn  _et al._ [2009] E. A. L. Henn, J. A. Seman, G. Roati, K. M. F. Magalhães, and V. S. Bagnato, Emergence of Turbulence in an Oscillating Bose-Einstein Condensate, Physical Review Letters 103, 045301 (2009).
  * Tsatos  _et al._ [2016] M. C. Tsatos, P. E. S. Tavares, A. Cidrim, A. R. Fritsch, M. A. Caracanhas, F. E. A. dos Santos, C. F. Barenghi, and V. S. Bagnato, Quantum Turbulence in Trapped Atomic Bose-Einstein Condensates, Physics Reports 622, 1 (2016).
  * Navon  _et al._ [2016] N. Navon, A. L. Gaunt, R. P. Smith, and Z. Hadzibabic, Emergence of a Turbulent Cascade in a Quantum Gas, Nature 539, 72 (2016).
  * Vinen and Niemela [2002] W. F. Vinen and J. J. Niemela, Quantum Turbulence, Journal of Low Temperature Physics 128, 167 (2002).
  * Madeira  _et al._ [2020] L. Madeira, M. A. Caracanhas, F. E. A. dos Santos, and V. S. Bagnato, Quantum Turbulence in Quantum Gases, Annual Review of Condensed Matter Physics 11, 37 (2020).
  * Barenghi  _et al._ [2023] C. F. Barenghi, H. A. J. Middleton-Spencer, L. Galantucci, and N. G. Parker, Types of Quantum Turbulence, AVS Quantum Science 5, 025601 (2023).
  * Griesmaier  _et al._ [2005] A. Griesmaier, J. Werner, S. Hensler, J. Stuhler, and T. Pfau, Bose-Einstein Condensation of Chromium, Physical Review Letters 94, 160401 (2005).
  * Lu  _et al._ [2011] M. Lu, N. Q. Burdick, S. H. Youn, and B. L. Lev, Strongly Dipolar Bose-Einstein Condensate of Dysprosium, Physical Review Letters 107, 190401 (2011).
  * Tang  _et al._ [2015] Y. Tang, N. Q. Burdick, K. Baumann, and B. L. Lev, Bose-Einstein Condensation of 162Dy and 160Dy, New Journal of Physics 17, 045006 (2015).
  * Aikawa  _et al._ [2012] K. Aikawa, A. Frisch, M. Mark, S. Baier, A. Rietzler, R. Grimm, and F. Ferlaino, Bose-Einstein Condensation of Erbium, Physical Review Letters 108, 210401 (2012).
  * Miyazawa  _et al._ [2022] Y. Miyazawa, R. Inoue, H. Matsui, G. Nomura, and M. Kozuma, Bose-Einstein Condensation of Europium, Physical Review Letters 129, 223401 (2022).
  * Eberlein  _et al._ [2005] C. Eberlein, S. Giovanazzi, and D. H. J. O’Dell, Exact Solution of the Thomas-Fermi Equation for a Trapped Bose-Einstein Condensate with Dipole-Dipole Interactions, Physical Review A 71, 033618 (2005).
  * Stuhler  _et al._ [2005] J. Stuhler, A. Griesmaier, T. Koch, M. Fattori, T. Pfau, S. Giovanazzi, P. Pedri, and L. Santos, Observation of Dipole-Dipole Interaction in a Degenerate Quantum Gas, Physical Review Letters 95, 150406 (2005).
  * Ticknor  _et al._ [2011] C. Ticknor, R. M. Wilson, and J. L. Bohn, Anisotropic Superfluidity in a Dipolar Bose Gas, Physical Review Letters 106, 065301 (2011).
  * Wenzel  _et al._ [2018] M. Wenzel, F. Böttcher, J.-N. Schmidt, M. Eisenmann, T. Langen, T. Pfau, and I. Ferrier-Barbut, Anisotropic Superfluid Behaviour of a Dipolar Bose-Einstein Condensate, Physical Review Letters 121, 030401 (2018).
  * Klaus  _et al._ [2022] L. Klaus, T. Bland, E. Poli, C. Politi, G. Lamporesi, E. Casotti, R. N. Bisset, M. J. Mark, and F. Ferlaino, Observation of Vortices and Vortex Stripes in a Dipolar Condensate, Nature Physics 18, 1453 (2022).
  * Bland  _et al._ [2023] T. Bland, G. Lamporesi, M. J. Mark, and F. Ferlaino, Vortices in Dipolar Bose-Einstein Condensates, Comptes Rendus Physique 24, 1 (2023).
  * Cooper  _et al._ [2005] N. R. Cooper, E. H. Rezayi, and S. H. Simon, Vortex Lattices in Rotating Atomic Bose Gases with Dipolar Interactions, Physical Review Letters 95, 200402 (2005).
  * Komineas and Cooper [2007] S. Komineas and N. R. Cooper, Vortex Lattices in Bose-Einstein Condensates with Dipolar Interactions Beyond the Weak-Interaction Limit, Physical Review A 75, 023623 (2007).
  * Zhang and Zhai [2005] J. Zhang and H. Zhai, Vortex Lattices in Planar Bose-Einstein Condensates with Dipolar Interactions, Physical Review Letters 95, 200403 (2005).
  * Martin  _et al._ [2017] A. M. Martin, N. G. Marchant, D. H. J. O’Dell, and N. G. Parker, Vortices and Vortex Lattices in Quantum Ferrofluids, Journal of Physics: Condensed Matter 29, 103004 (2017).
  * Cai  _et al._ [2018] Y. Cai, Y. Yuan, M. Rosenkranz, H. Pu, and W. Bao, Vortex Patterns and the Critical Rotation Frequency in Rotating Dipolar Bose-Einstein Condensates, Physical Review A 98, 023610 (2018).
  * Klawunn  _et al._ [2008] M. Klawunn, R. Nath, P. Pedri, and L. Santos, Transverse Instability of Straight Vortex Lines in Dipolar Bose-Einstein Condensates, Physical Review Letters 100, 240403 (2008).
  * Abad  _et al._ [2009] M. Abad, M. Guilleumas, R. Mayol, M. Pi, and D. M. Jezek, Vortices in Bose-Einstein Condensates with Dominant Dipole-Dipole Interactions, Physical Review A 79, 063622 (2009).
  * Mulkerin  _et al._ [2014] B. C. Mulkerin, R. M. W. van Bijnen, D. H. J. O’Dell, A. M. Martin, and N. G. Parker, Vortices in the Two-Dimensional Dipolar Bose Gas, Journal of Physics: Conference Series 497, 012025 (2014).
  * Mulkerin  _et al._ [2013] B. C. Mulkerin, R. M. W. van Bijnen, D. H. J. O’Dell, A. M. Martin, and N. G. Parker, Anisotropic and Long-Range Vortex Interactions in Two-Dimensional Dipolar Bose Gases, Physical Review Letters 111, 170402 (2013).
  * Gautam [2014] S. Gautam, Dynamics of the Corotating Vortices in Dipolar Bose–Einstein Condensates in the Presence of Dissipation, Journal of Physics B: Atomic, Molecular and Optical Physics 47, 165301 (2014).
  * Jia  _et al._ [2018] L. Jia, A.-B. Wang, and S. Yi, Low-Lying Excitations of Vortex Lattices in Condensates with Anisotropic Dipole-Dipole Interaction, Physical Review. A 97, 043614 (2018).
  * Zhao [2021] Q. Zhao, Effects of Dipole–Dipole Interaction on Vortex Motion in Bose–Einstein Condensates, Journal of Low Temperature Physics 204, 1 (2021).
  * Prasad  _et al._ [2019] S. B. Prasad, T. Bland, B. C. Mulkerin, N. G. Parker, and A. M. Martin, Vortex Lattice Formation in Dipolar Bose-Einstein Condensates via Rotation of the Polarization, Physical Review A 100, 023625 (2019).
  * Navon  _et al._ [2021] N. Navon, R. P. Smith, and Z. Hadzibabic, Quantum Gases in Optical Boxes, Nature Physics 14, 1334 (2021).
  * Juhász  _et al._ [2022] P. Juhász, M. Krstajić, D. Strachan, E. Gandar, and R. P. Smith, How to Realize a Homogeneous Dipolar Bose Gas in the Roton Regime, Physical Review A 105, L061301 (2022).
  * Bland  _et al._ [2018] T. Bland, G. W. Stagg, L. Galantucci, A. W. Baggaley, and N. G. Parker, Quantum Ferrofluid Turbulence, Physical Review Letters 121, 174501 (2018).
  * Lahaye  _et al._ [2009] T. Lahaye, C. Menotti, L. Santos, M. Lewenstein, and T. Pfau, The Physics of Dipolar Bosonic Quantum Gases, Reports on Progress in Physics 72, 126401 (2009).
  * Chomaz  _et al._ [2009] L. Chomaz, I. Ferrier-Barbut, F. Ferlaino, B. Laburthe-Tolra, B. L. Lev, and T. Pfau, Dipolar Physics: A Review of Experiments with Magnetic Quantum Gases, Reports on Progress in Physics 86, 026401 (2009).
  * Schmitt  _et al._ [2016] M. Schmitt, M. Wenzel, F. Böttcher, I. Ferrier-Barbut, and T. Pfau, Self-Bound Droplets of a Dilute Magnetic Quantum Liquid, Nature 539, 259 (2016).
  * Chomaz  _et al._ [2016] L. Chomaz, S. Baier, D. Petter, M. J. Mark, F. Wächtler, L. Santos, and F. Ferlaino, Quantum-Fluctuation-Driven Crossover from a Dilute Bose-Einstein Condensate to a Macrodroplet in a Dipolar Quantum Fluid, Physical Review X 6, 041039 (2016).
  * Böttcher  _et al._ [2019] F. Böttcher, J.-N. Schmidt, M. Wenzel, J. Hertkorn, M. Guo, T. Langen, and T. Pfau, Transient Supersolid Properties in an Array of Dipolar Quantum Droplets, Physical Review X 9, 011051 (2019).
  * Chomaz  _et al._ [2019] L. Chomaz, D. Petter, P. Ilzhöfer, G. Natale, A. Trautmann, C. Politi, G. Durastante, R. M. W. van Bijnen, A. Patscheider, M. Sohmen, M. J. Mark, and F. Ferlaino, Long-Lived and Transient Supersolid Behaviors in Dipolar Quantum Gases, Physical Review X 9, 021012 (2019).
  * Tanzi  _et al._ [2019] L. Tanzi, E. Lucioni, F. Famà, J. Catani, A. Fioretti, C. Gabbanini, R. N. Bisset, L. Santos, and G. Modugno, Observation of a Dipolar Quantum Gas with Metastable Supersolid Properties, Physical Review Letters 122, 130405 (2019).
  * Schützhold  _et al._ [2006] R. Schützhold, M. Uhlmann, Y. Xu, and U. R. Fischer, Mean-Field Expansion in Bose–Einstein Condensates with Finite-Range Interactions, International Journal of Modern Physics B 20, 3555 (2006).
  * Lima and Pelster [2012] A. R. P. Lima and A. Pelster, Beyond Mean-Field Low-Lying Excitations of Dipolar Gases, Physical Review A 86, 063609 (2012).
  * Bisset  _et al._ [2016] R. N. Bisset, R. M. Wilson, D. Baillie, and P. B. Blakie, Ground-State Phase Diagram of a Dipolar Condensate with Quantum Fluctuations, Physical Review A 94, 033619 (2016).
  * Barenghi and Parker [2016] C. F. Barenghi and N. G. Parker,  _A Primer on Quantum Fluids_ , Springer Briefs in Physics (Springer, 2016).
  * Bao and Cai [2013] W. Bao and Y. Cai, Mathematical Theory and Numerical Methods for Bose-Einstein Condensation, Kinetic and Related Models 6, 1 (2013).
  * Kerr [2011] R. M. Kerr, Vortex Stretching as a Mechanism for Quantum Kinetic Energy Decay, Physical Review Letters 106, 224501 (2011).
  * Rorai  _et al._ [2013] C. Rorai, K. R. Sreenivasan, and M. E. Fisher, Propagating and Annihilating Vortex Dipoles in the Gross-Pitaevskii Equation, Physical Review B 88, 134522 (2013).
  * Note [1] We use only 2 gridpoints along $z$ as uniformity is assumed along this dimension.
  * Lee  _et al._ [2021] A.-C. Lee, D. Baillie, and P. B. Blakie, Numerical Calculation of Dipolar-Quantum-Droplet Stationary States, Physical Review Research 3, 013283 (2021).
  * Note [2] More specifically, to ensure that the same parity conditions are observed at both boundaries along a given axis, we employ the DCT-I transform [82].
  * Yi and You [2000] S. Yi and L. You, Trapped Atomic Condensates with Anisotropic Interactions, Physical Review A 61, 041604(R) (2000).
  * Yi and Pu [2006] S. Yi and H. Pu, Vortex Structures in Dipolar Condensates, Physical Review A 73, 061602(R) (2006).
  * Foster  _et al._ [2010] C. J. Foster, P. B. Blakie, and M. J. Davis, Vortex Pairing in Two-Dimensional Bose Gases, Physical Review A 81, 023623 (2010).
  * White  _et al._ [2012] A. C. White, C. F. Barenghi, and N. P. Proukakis, Creation and Characterization of Vortex Clusters in Atomic Bose-Einstein Condensates, Physical Review A 86, 013635 (2012).
  * O’Riordan [2017] L. J. O’Riordan,  _Non-Equilibrium Vortex Dynamics in Rapidly Rotating Bose–Einstein Condensates_ , Ph.D. thesis, Okinawa Institute of Science and Technology (2017).
  * Moon and Spencer [1971] P. Moon and D. E. Spencer,  _Field Theory Handbook_, 2nd ed. (Springer-Verlag, 1971).
  * Sabari and Kumar [2018] S. Sabari and R. K. Kumar, Effect of an Oscillating Gaussian Obstacle in a Dipolar Bose-Einstein Condensate, European Physical Journal D 72, 48 (2018).
  * Tsuchiya  _et al._ [2008] S. Tsuchiya, F. Dalfovo, and L. Pitaevskii, Solitons in Two-Dimensional Bose-Einstein Condensates, Physical Review A 77, 045601 (2008).
  * Tikhonenkov  _et al._ [2008] I. Tikhonenkov, B. A. Malomed, and A. Vardi, Anisotropic Solitons in Dipolar Bose-Einstein Condensates, Physical Review Letters 100, 090406 (2008).
  * Note [3] A. Villois and D. Proment, private communication.
  * Groszek  _et al._ [2018] A. J. Groszek, D. M. Paganin, K. Helmerson, and T. P. Simula, Motion of Vortices in Inhomogeneous Bose-Einstein Condensates, Physical Review A 97, 023617 (2018).
  * Richaud  _et al._ [2022] A. Richaud, P. Massignan, V. Penna, and A. L. Fetter, Dynamics of a Massive Superfluid Vortex in $r^{k}$ Confining Potentials, Physical Review A 106, 063307 (2022).
  * Richaud  _et al._ [2023] A. Richaud, G. Lamporesi, M. Capone, and A. Recati, Mass-Driven Vortex Collisions in Flat Superfluids, Physical Review A 107, 053317 (2023).
  * Klawunn and Santos [2009] M. Klawunn and L. Santos, Phase Transition from Straight into Twisted Vortex Lines in Dipolar Bose–Einstein Condensates, New Journal of Physics 11, 055612 (2009).
  * Wilson  _et al._ [2009] R. M. Wilson, S. Ronen, and J. L. Bohn, Stability and Excitations of a Dipolar Bose-Einstein Condensate with a Vortex, Physical Review A 79, 013621 (2009).
  * Rorai  _et al._ [2016] C. Rorai, J. Skipper, R. M. Kerr, and K. R. Sreenivasan, Approach and Separation of Quantised Vortices with Balanced Cores, Journal of Fluid Mechanics 808, 641 (2016).
  * Smirnov [1964] V. I. Smirnov,  _Part 2: Complex Variables – Special Functions_ , 1st ed., A Course of Higher Mathematics, Vol. III (Pergamon Press, 1964).
  * Caselle  _et al._ [1996] M. Caselle, F. Gliozzi, U. Magnea, and S. Vinti, Width of Long Colour Flux Tubes in Lattice Gauge Systems, Nuclear Physics B 460, 397 (1996).
  * Billam  _et al._ [2014] T. P. Billam, M. T. Reeves, B. P. Anderson, and A. S. Bradley, Onsager-Kraichnan Condensation in Decaying Two-Dimensional Quantum Turbulence, Physical Review Letters 112, 145301 (2014).
  * Weiss and McWilliams [1991] J. B. Weiss and J. C. McWilliams, Nonergodicity of Point Vortices, Physics of Fluids A: Fluid Dynamics 3, 835 (1991).
  * Griffin  _et al._ [2020] A. Griffin, V. Shukla, M.-E. Brachet, and S. Nazarenko, Magnus-Force model for Active Particles Trapped on Superfluid Vortices, Physical Review A 101, 053601 (2020).
  * Britanak  _et al._ [2007] V. Britanak, P. C. Yip, and K. R. Rao,  _Discrete Cosine and Sine Transforms: General Properties, Fast Algorithms and Integer Approximations_ , 1st ed. (Academic Press, 2007).