Vortex pair dynamics in three-dimensional homogeneous dipolar superfluids

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 axisa 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 mag-netic 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 investigationsincluding 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 3D 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 µ d which is polarised uniformly by an applied magnetic field parallel to d.For this dipolar BEC, the dGPE is given by [31,44,45], Here, the two-body short-ranged interaction strength is g = 4π ̵ h 2 a s /m with a s the scattering length of the atom-atom scattering potential, n = |ψ| 2 is the atomic density of the condensate, and µ the chemical potential of the condensate which fixes the normalisation of ψ and satisfies the self-consistency relation, (2) The parameter ε dd = mµ 0 µ 2 /(12π ̵ h 2 a s ), where µ 0 is the permeability of free space, serves as an effective dipolar-tocontact 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 dd (r) = 1 4π [ and has a Fourier transform given by 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, 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 µ g = gn 0 {1 + ε dd lim q→0 [3( d ⋅ q) 2 − 1]}.Thus, in free space the solutions to Eq. ( 1) are unstable when ε 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 < ε 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 µ g , length in units of the healing length ξ = ̵ h/ √ mµ g , time in units of ̵ h/µ g , and ψ in units of √ n 0 [54].
Using these natural units, Eq. ( 1) can be written in dimensionless variables as [24,43] where 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 ψ in the Madelung form ψ = √ ne iS and interpreting the superfluid phase S as a scaled velocity potential for the superfluid velocity v such that v = ∇S .Under these transformations, the dGPE may be re-expressed as a pair of superfluid hydrodynamic equations [44], from which we shall extract considerable intuition about the expected form of numerical solutions to the dGPE.For instance, the quantization of circulation of ψ [1], viz.
is inherently a manifestation of the single-valued nature of ψ while the fact that the superfluid density vanishes at the core of a quantum vortex, where q ≠ 0, is a logical consequence of finding a solution of Eq. (7) for n(r) where v(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 righthand-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 ψ 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 3D dipolar BECs have received little attention when compared to those in trapped 3D systems or those that are uniform in quasi-2D 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, ε dd , and the dipole alignment axis 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 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 Γ, we employ Neumann (reflecting) boundary conditions in the x-y plane [56,57], where the normal derivative of ψ 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, ∆x = ∆y = ∆z = 25/256, along each direction [58].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 ∆k i = π/[(N i − 1)∆r i ] in the k x -k y plane and ∆k z = 2π/(N z ∆z).Then, the spatial derivatives may be computed using the discrete cosine transform [59] (DCT) [60] 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 ñ(k 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 ε 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 ∈ 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 ψ 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 ↦ it and renormalizing ψ at each timestep such that the mean value of the density obeys ⟨n⟩ = ⟨|ψ| 2 ⟩ = 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(ρ) is proportional to ρ 2 as ρ, the distance from the centre of the vortex core, tends to zero, and approaches µ/µ g when ρ ≫ ξ [1].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, θ d .Since V dd (r) ∝ 1 − 3 cos 2 θ d , it is energetically preferable for the two dipole moments to minimise the angle θ, 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 FIG. 1. Cross-sections of the density profiles, n(r), in the x-y plane for a vortex of circulation 2π where the magnetic field is polarised along the x-axis and ε dd = 0 (a), 0.3 (b), 0.6 (c), and 0.9 (d).The density n(ρ) ∼ ρ 2 as ρ → 0 and approaches a constant value far from the vortex core.For nonzero ε dd , the core is elongated along the dipolar axis.
naturally leads to higher values of the density along this axis in stationary states of the Euler equation, Eq. ( 8) [22,23,61].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] 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 π/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 [27,62] and quasi-two-dimensional uniform [35,36] 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 d = 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 ε 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 ε 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 [35,36]; 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 the discrepancy of the density at a given point from the angleaveraged density at that radius, and δS (r) = S (r)−S 0 (x+iy), the discrepancy of the phase from the incompressible solution, satisfying Neumann boundary conditions, given by Eq. (A1) in Appendix A 1. Furthermore, we evaluate ∇ ⋅ v, a measure of the compressibility of the superfluid velocity fields of the stationary solutions, and Φ dd (r) = ⟨V dd (r)⟩ ≡ n(r) ⊛ V dd (r), the local expectation value of the dipolar interaction potential.Each of these quantities are presented, below, in Fig. 2 for the values ε dd = {0.3,0.6, 0.9}.From Fig. 2, it is immediately clear that the density profile anisotropy δn(r) in the first row is inversely correlated with Φ(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 Φ 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 δS in the second row demonstrate that the vortex phase is modified by dipolar interaction.Interestingly, both δn and δS take the form of lobes reminscent of d-wave orbitals, a symmetry obeyed by the DDI, with the lobes corresponding to δn being out of phase compared to those of δS by π/4.We also observe that the fact that the maximum magnitude of δn increases with ε dd is accompanied by a similar behaviour exhibited by δ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 FIG. 2. Cross-sections of the axial density anisotropy δn (first row), phase anisotropy δS (second row), local dipolar interaction energy Φ dd (third row), and compressibility ∇ ⋅ v (fourth row) as defined in the main text for ε dd = 0.3 (first column, a-d), 0.6 (second column, e-h), and 0.9 (third column, i-l).An inverse correlation between δn and Φ dd is evident and the maximum magnitude of both δn and δS increases with ε dd .N.B. for the sake of clarity, we have set ∇ ⋅ v -which diverges at the phase singularity -to zero at r = 0. vortex core in the fourth row of Fig. 2. In particular, this behaviour exhibited by δ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 Φ, 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.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 Γ = 2π at the positions (x, y) = (0, ±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 ∆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πq, q ∈ 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 r i , the positions of each vortex evolve according to the following coupled set of equations, For a pair of same-signed vortices of circulation Γ i = Γ j = 2π 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 ω = 2ẑ/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 r = 0 at a radius s/2 and period 4π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ξ (a), s = 6.25ξ (b), and s = 12.5ξ (c), with ε dd ∈ {0, 0.3, 0.6, 0.9} and 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 π/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 ε 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 ε dd .When ε 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 ŷ.A larger degree of ellipticity is seen for ε dd = 0.6, though, and for ε 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, Φ 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 ε dd .
It is important to note that in previous theoretical studies of vortex dynamics in quasi-two-dimensional systems, likesigned vortex pairs were observed to undergo a similar transition in orbit from circular, through elliptical, to oval-like with increasing ε dd when the dipole polarisation had a nonzero projection in the x-y plane [36].The same studies in the quasi-2D 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 Γ 1 = −Γ 2 = 2π one finds that ṙ12 is now conserved and that Ṙ12 = (ẑ × r12 )/r 12 .Thus, the vortex-antivortex pair undergoes translation along the binormal axis parallel to ẑ × r12 at the constant speed v = 1/r 12 .For instance, for a pair of vortices with circulations ±2π at the initial positions (x, y) = [x(0), y(0) ± s/2], In each plot, we have ε 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 ε dd (s) is apparent.
the trajectories of the vortices as predicted by the point-vortex model are 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 Γ = ±2π at the positions (x, y) = (0, ±s/2) with s = 6.25ξ 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 ∆x = ∆y = ∆z = 25/256, and the spatial derivatives are computed via FFTs; the wavenumber grids have spacings ∆k i = 2π/(N i ∆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 ε dd from 0 to 0.9 and allow d to be parallel to x, ŷ and ẑ.The results of these simulations are presented as a plot of v x = ẋi as a function of the full range of ε dd in Fig. 4.
We have found that for s ≫ ξ, 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, ⟨∆y⟩, after an equilibrium has been attained; we note that, in the point-vortex limit, the quantity v x ⟨∆y⟩ should be equal to unity.Instead, we find that a strong dipolar dependence exists.For the full range of ε dd that we probe, v x ⟨∆y⟩ increases monotonically with ε 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 ε dd and undergoes a maximum (which exceeds the corresponding value of v x ⟨∆y⟩ with the dipole moments polarised along the y-axis) before decreasing for larger ε 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 ε 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-twodimensional 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 pre-dicted 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 ∼ 4 healing lengths, the relationship between the translational velocities, dipole orientation and dipolar interaction strength is found to be essentially universal with vortexantivortex 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 vortexantivortex 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 [67-69] [70].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 semianalytical 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 pointvortex 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.In a 3D domain with Neumann boundary conditions in the x-y plane and periodicity in the z-direction, the infinitedomain 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 2D Poissonian Green function in this domain which, in turn, represents the superfluid equivalent of the stream function.To wit, in a rectangular grid x ∈ [0, L x ], y ∈ [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]  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 v, the superfluid velocity.Thus, we can only consider ensembles of vortexantivortex 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 ±2π is located at the position (x ± , y ± ), 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 ∞ (r) = arctan(y − y + , x − x + ) − arctan(y − y − , x − x − ) with (x ± , y ± ) replaced by the coordinates of a given periodic copy of the original vortex pair.For a grid with x ∈ [0, L x ) and y ∈ [0, L y ), the solution for the superfluid phase has been provided in the literature by [79], where X ± = x − x ± and Y ± = y − y ± .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 ∑ ∞ p=−∞ with a finite series ∑ 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

FIG. 3 .
FIG. 3. Trajectories of a vortex of circulation Γ = 2π 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 ε 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 ε dd (s) is apparent.

FIG. 4 .
FIG. 4. Velocities along the x-axis of a vortex-antivortex pair of circulation Γ = ±2π with an initial separation of d = 3.125ξ (a), 6.25ξ (b), and 12.5ξ (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 = vx(ε dd )/vx((ε dd = 0) and multiply them by the scaled vortex-antivortex separation along the y-axis, ⟨∆y⟩ ′ = ⟨∆y(ε dd )⟩/⟨∆y(ε dd = 0)⟩.