Dynamics of a Vortex Dipole in a Holographic Superfluid

We use holography to investigate the dynamics of a vortex-anti-vortex dipole in a strongly coupled superfluid in 2+1 dimensions. The system is evaluated in numerical real-time simulations in order to study the evolution of the vortices as they approach and eventually annihilate each other. A tracking algorithm with sub-plaquette resolution is introduced which permits a high-precision determination of the vortex trajectories. With the increased precision of the trajectories it becomes possible to directly compute the vortex velocities and accelerations. We find that in the holographic superfluid the vortices follow universal trajectories independent of their initial separation, indicating that a vortex-anti-vortex pair is fully characterized by its separation. Subtle non-universal effects in the vortex motion at early times of the evolution can be fully attributed to artifacts due to the numerical initialization of the vortices. We also study the dependence of the dynamics on the temperature of the superfluid.


Introduction
Experimental advances in studying the non-equilibrium dynamics of isolated quantummany-body systems have in recent years led to the quest for a better theoretical understanding of such dynamics. Systems of interest in this context include ultracold Bose gases [1][2][3][4][5][6], semi-conductor exciton-polariton superfluids [7][8][9][10][11] and topological condensed matter systems [12]. The non-equilibrium dynamics of such systems can be observed after a quench into a far-from-equilibrium initial state, as it can be realized experimentally by a sudden change of a parameter of the system, e. g. temperature or the scattering length of the constituent particles. During the relaxation process and before thermal equilibrium is reached the system's dynamics is generically strongly nonlinear and exhibits properties of quantum turbulence, the quantum analogue of classical turbulence. This line of research adds new aspects to the general phenomena of quantum and classical turbulence which have attracted considerable attention since decades and include, among others, universal scaling of correlation functions [13][14][15][16] with characteristic energy cascades in the system, quasi-stationary flow [17], and the observation of non-thermal fixed points [18][19][20][21][22][23] (see [24] for a review).
The strongly nonlinear nature of the dynamics poses a great challenge as it renders a theoretical description very intricate. Even if the system's constituents are weakly coupled, commonly used methods such as perturbation theory break down as a consequence of strong correlations that build up dynamically, often related to the occurrence of highly nonlinear topological excitations in the system. These include in particular quantized vortices, the quantum analogues of classical eddies. In modern experiments with Bose-Einstein-condensed ultracold quantum gases vortex ensembles can even be produced by stirring the condensate with a laser [25,26]. Quantized vortices are known to play a crucial role also in superfluid helium [27] where strong coupling prevails in the system and the theoretical description is even more challenging. In the present work we will be concerned with the dynamics of a vortex-anti-vortex pair which is the most basic configuration that allows one to study the interaction of topological defects.
We will use the theoretical framework of holography or gauge/gravity duality for our study. It offers a promising pathway towards a better understanding of non-equilibrium dynamics and is particularly well suited for the description of strongly coupled quantum systems. In holography, a strongly coupled quantum field theory without gravity in d dimensions is mapped to a theory with classical gravity in d + 1 dimensions. The quantum behavior of the field theory is thus encoded in a classical gravitational theory in a higherdimensional spacetime. A variety of holographic dualities has been devised since such a duality was first proposed in the form of the AdS/CFT correspondence [28][29][30]. AdS/CFT posits that N = 4 super Yang-Mills (SYM) theory, a supersymmetric and conformal gauge theory (CFT) in 3 + 1 dimensions, is dual to type IIB string theory on AdS 5 × S 5 , where AdS 5 stands for an Anti-de Sitter spacetime and S 5 denotes a five-sphere. A key feature is that the AdS/CFT correspondence simplifies in the limit of strong coupling and large number of colors on the field theory side, corresponding to which the string theory side simplifies to a weakly coupled theory of Einstein gravity. Generalizations to dualities between other strongly coupled quantum theories and weakly coupled gravitational theories have become a powerful tool to address many problems in strongly coupled and hence manifestly nonperturbative quantum field theories. The quantum theory in d dimensions is usually referred to as the 'boundary theory' as one can imagine it to 'live' on the boundary of the (d + 1)-dimensional AdS spacetime. The AdS spacetime away from the boundary, on the other hand, is referred to as the 'bulk'.
Holographic dualities can in particular be used to describe quantum field theories at finite temperature. Here, typically, the dual gravitational representation includes a black hole in the bulk of the AdS space, with its Hawking temperature being equal to the temperature in the boundary theory. Similarly, a chemical potential in the field theory corresponds to a charged black hole in the bulk. Holography has been applied to a variety of physical systems (see e. g. [31] for a review), ranging from quark-gluon plasma physics, see e. g. [32], to condensed matter systems, see for instance [33], including the pioneering work on holographic superconductors and superfluids in [34][35][36] (see [37] for a review).
Specifically, it has been found [34][35][36] that an Abelian Higgs model coupled to gravity on a (3+1)-dimensional Anti-de Sitter spacetime with a black hole has a dual interpretation in terms of a superconductor or superfluid (for a review, see [38]). Superfluidity is associated with the breaking of a global U(1) symmetry. In the holographic model, this is incorporated by the breaking of a local U(1) symmetry and the associated condensation of the scalar field in the bulk. Remarkably, the simple Abelian Higgs model with classical gravity encodes the quantum behavior of the superfluid, even if strong correlations are present in the system. It is therefore well suited for studying the real-time dynamics of the superfluid far from equilibrium.
Extensive analytical as well as numerical studies of the holographic superfluid have been performed. Among them are studies of the details of the phase transition into the superfluid phase 1 [36,42], computations of the speed of sound and the investigation of sound waves, see e. g. [36,[43][44][45][46], and studies of topological defects in the system, see e. g. [47][48][49][50][51][52][53][54][55][56][57][58][59][60][61][62][63]. With regard to topological defects, much attention has been devoted to the study of vortices and in particular ensembles of vortices and anti-vortices. While quantized vortices constitute an intense field of research in general, they play a particular role in the holographic superfluid as the main sites of energy and momentum dissipation in this system [51]. To wit, they have a geometric interpretation in that tubes emanating from the vortex cores at the boundary punch holes through the scalar charge cloud in the bulk so that energy can flow from the boundary to the black hole where it is absorbed, see [51].
In the present work, we will study the annihilation process of a single vortex-antivortex pair, or vortex dipole, in the holographic superfluid. By considering just one vortex and one anti-vortex we have a clean environment for investigating their mutual interaction. At the same time, vortex pairs can be viewed as basic building blocks for the interactions in larger ensembles of vortices and are thus relevant in a more general context. Our method builds on previous work in [54] where the vortex dynamics of an ensemble of vortices was studied for certain classes of quench-like initial conditions, and we use the same numerical implementation. While the annihilation process of a vortex dipole in a holographic superfluid has been studied before [59], we are able to study new aspects of this system by using new methods. Specifically, by introducing an advanced tracking technique, we drastically increase the precision of the vortex trajectories that we extract from the numerical simulations. This gain in precision allows us to compute directly the vortex velocities and accelerations. We can then infer universal properties of the motion of the vortex dipole and can cleanly identify a power law for the behavior of the vortexanti-vortex separation. We also compare the vortex trajectories to a general point-vortex model from the literature. Furthermore, we study the temperature dependence of the vortex dynamics. Finally, we discuss some potential further applications of the methods for tracking vortices developed here.
Our paper is organized as follows. In section 2 we review the gravitational model for the holographic superfluid. We discuss its equations of motion and their numerical implementation used for the real-time simulation of the system. Section 3 is concerned with the study of the vortex-anti-vortex pair in the holographic superfluid. We describe the numerical initialization of vortices and introduce a new method to determine their positions during the evolution of the system with high precision. Subsequently, we present our results for the vortex trajectories from the initial separation to the annihilation. Here we also discuss universal properties of the motion of the vortex dipole and compute velocities and accelerations of the vortices. In section 4 we study the dependence of the vortex dynamics on the temperature of the superfluid. We close with a summary and outlook in section 5. Some technical details are presented in the appendices.

Holographic Superfluid
In this section we review the 3 + 1-dimensional gravitational model dual to a superfluid in 2 + 1 dimensions.

Holographic Model
In his seminal work [34], Gubser showed that, if coupled to gravity with a negative cosmological constant, even a simple Abelian Higgs model without self-interaction exhibits spontaneous symmetry breaking. This laid the ground for a holographic description of superconductivity and superfluidity as subsequently developed in [35,36]. 2 The spontaneous breaking of the local U (1) symmetry in the AdS bulk occurs due to the negative contribution of the gauge field to the effective scalar mass in the vicinity of the black-hole horizon which drives the system unstable and leads to the formation of a scalar charge cloud. In the 2 The holographic superfluid was devised in [36] after the dual description of a superconductor was introduced and discussed in [35]. From the field-theory perspective, the main difference between a superfluid and a superconductor is that for a superfluid a global symmetry is spontaneously broken whereas for a superconductor a local symmetry is spontaneously broken. The differences between the holographic superfluid and the holographic superconductor with regard to topological vortex defects are discussed in [52], for example. dual field theory this corresponds to the breaking of a global U (1) symmetry and a phase transition from the normal fluid to a superfluid state. The gravity model holographically describing a (2 + 1)-dimensional superfluid is given by the action [34] where g is the determinant of the metric g M N (M, N = 0, . . . , 3), R its associated Ricci scalar, κ Newton's constant in 4 spacetime dimensions, and Λ = −3/L 2 AdS the cosmological constant of the (3 + 1)-dimensional Anti-de Sitter spacetime with curvature radius L AdS . The matter part L matter of the Lagrangian is given by where the charge q of the scalar field has been pulled out of the Lagrangian L matter by an appropriate rescaling of the fields. Using this convention makes evident that q controls the effect of the matter part of the model on the dynamical gravity, cf. (2.1). The mass of the scalar field Φ is m. We use ∇ M to denote the Levi-Civita connection of the metric. The gravity model is thus given by a simple Abelian Higgs model formulated on a (3 + 1)dimensional gravity background with minimal coupling between the complex scalar field Φ and the gauge connection A M . The kinetic term for the Abelian gauge field is formulated in terms of its field strength tensor Note that no terms of higher order in the scalar field than the quadratic mass term are needed. This is different from the usual Higgs mechanism in the Standard Model where the Φ 4 interaction term is essential for spontaneous symmetry breaking to occur. In contrast to flat Minkowski spacetime, the AdS-type spacetime allows for an instability to occur even if there is no self-interaction of the scalar field. The squared scalar mass m 2 may even be negative without the symmetry being broken as long as it is above the Breitenlohner-Freedman bound [64][65][66]. At low temperatures, the effective scalar mass can be lowered below this bound by the negative contribution of the bulk gauge field, thus giving rise to an instability and symmetry breaking. In this work we choose m 2 = −2/L 2 AdS which is above the Breitenlohner-Freedman bound, and we set L AdS = 1 for convenience.
The interpretation of this model according to the holographic dictionary is the following. The Abelian gauge field in the bulk A M is dual to a global conserved U (1) current operator j µ in the boundary field theory (µ = 0, 1, 2). Similarly, the complex scalar field Φ with mass m and charge q, to which the gauge field A M is dynamically coupled, is dual to a complex scalar field operator Ψ. This operator plays a crucial role as its quantum expectation value ψ = Ψ is the order parameter for the phase transition from a normal fluid to the superfluid. In a usual field-theory context, superfluidity is described by a complex scalar field whose vacuum expectation value assumes non-zero values in the symmetry-broken phase, corresponding to the existence of a superfluid condensate.
Instead of considering the full set of equations of motion corresponding to the action (2.1) we simplify the model by taking the so-called probe approximation. This means that the gravity part of the action is solved independently of the matter fields, i. e. the backreaction of the matter fields on the metric is neglected. This approximation is justified if the scalar charge q is large. The matter part of the model is then solved in the fixed gravity background. For a detailed discussion of the probe limit and its validity see [35] and, in the context of vortex defects, [67].
The pure gravity part which comprises only the first two terms in the action (2.1), i. e. the Ricci scalar and the cosmological constant, can be solved straightforwardly. The vacuum Einstein equations are given by They are solved by a (3 + 1)-dimensional Schwarzschild Anti-de Sitter spacetime that has a planar black hole with horizon situated at the bulk position z = z h . In this work we choose to work with infalling Eddington-Finkelstein coordinates in which the line element is given by with the horizon function The choice of infalling Eddington-Finkelstein coordinates is convenient for two reasons. Firstly, unlike the usual Schwarzschild coordinates they do not exhibit a coordinate singularity at the black-hole horizon. Secondly, all null geodesics falling into the black hole are of radial nature. These two features are particularly useful for numerical simulations in holography [68] due to the ensuing natural implementation of boundary conditions. In these coordinates, z is the holographic (radial) coordinate and (t, x) = (t, x, y) are the time and remaining spatial directions, which also make up the coordinates of the dual field theory that can be thought of as living on the boundary at z = 0 of the Anti-de Sitter spacetime. As a radial coordinate z is taken to be positive and one refers to the spacetime for z > 0 as the bulk. The Hawking temperature associated with the metric (2.4) is and it coincides with the temperature of the boundary theory. Taking the background metric to be fixed one can vary the action (2.1) with respect to the matter fields to obtain their equations of motion, with the current 9) and the general covariant derivative D M = ∇ M − iA M . This is a classical system of coupled partial differential equations on the curved 4-dimensional Anti-de Sitter spacetime (2.4). The Maxwell and Klein-Gordon equations, (2.7) and (2.8) respectively, are coupled through the electromagnetic current (2.9). This classical system of equations of motion together with the appropriate boundary conditions encode the quantum behavior of the dual field theory. We are interested in computing the full quantum evolution of the order parameter field ψ = Ψ , where Ψ is the operator dual to the complex scalar field Φ in the bulk. Finding the time evolution of the expectation value ψ thus reduces to solving the classical equations of motion for Φ. To solve the equations of motion, we first need to impose boundary conditions for the gauge field A M as well as for the scalar field Φ. The boundary condition for the temporal component A t of the gauge field at the conformal boundary z = 0 fixes the chemical potential µ in the dual quantum theory, i. e., where µ is independent of the superfluid's coordinates (t, x). The only parameter of the holographic superfluid, uniquely characterizing it, is the dimensionless ratio µ/T ∼ (µz h ) (cf. (2.6)). One finds that the scalar field Φ forms a charge condensate above a critical value (µz h ) c = 4.06371 (see e. g. [36]) which sets the system in the superfluid phase. It is then straightforward to relate µ/T ∼ (µz h ) to the temperature ratio T /T c , where T c denotes the critical temperature, We fix our units by setting z h = 1. This also fixes the superfluid's temperature to T = 3 4π (see (2.6)). Therefore, µ is the only free parameter that controls the phase transition and hence the temperature ratio T /T c according to (2.11). Henceforth, we will characterize and discuss the superfluid in terms of its temperature ratio T /T c as this appears to be the more intuitive parameter. But we keep in mind that what we actually vary is the chemical potential.
According to the holographic dictionary, the scalar expectation value ψ = Ψ is identified as the next-to-leading order term in the near-boundary expansion of the scalar field Φ, where η(t, x) constitutes a source for the scalar field which has to be set to zero, i. e. η ≡ 0, to not break the U (1)-symmetry explicitly. This condition can be expressed as Similarly, all sources for the spatial gauge field components, which in the holographic interpretation correspond to an external superfluid flow, are switched off, i. e. A x (t, x, z = 0) = A y (t, x, z = 0) = 0. At the black-hole horizon, z = z h , the gauge field components must satisfy regularity conditions, given by Similarly, the scalar field Φ must behave regularly at the horizon which is intrinsically implemented by the choice of infalling Eddington-Finkelstein coordinates. More details on the equations of motion and the boundary conditions are given in appendix A. The scalar order parameter ψ is a complex field whose absolute value squared gives the density of the superfluid, n(t, x) = |ψ(t, x)| 2 . (2.14) We denote the homogeneous background density of the superfluid in thermal equilibrium, where no vortices are present, as n 0 . The phase ϕ(t, x) = arg(ψ(t, x)) of the order parameter encodes the superfluid velocity via v s ( In the probe approximation the black hole constitutes an infinitely sized static heat bath for the superfluid. Due to its presence, the superfluid is equipped with an energy and momentum dissipation mechanism caused by modes falling into the black hole. In this sense, the black hole can loosely be interpreted as a normal component connected to the superfluid condensate in the sense of Tisza's two-fluid model [69]. In the case with back-reaction it has been shown that to non-dissipative order 3 the holographic superfluid is indeed governed by a relativistic version of Tisza's two-fluid model [70]. It appears likely that a similar identification holds, at least approximately, also at dissipative order, cf. [70]. Also for the non-back-reacted system an analogy to the two-fluid model could be expected, but a rigorous identification has not been established so far. Keeping this caveat in mind, it can be useful to view the results of the following sections in light of the two-fluid picture.

Remarks on the Numerical Implementation
We are interested in performing a real-time simulation of the dynamics of the holographic superfluid. The coupled system of equations of motion (see appendix A) cannot be solved analytically and hence must be dealt with numerically. The superfluid is described by the order parameter field ψ which can be extracted from the solution of the equations of motion according to the near-boundary expansion (2.12). For T /T c < 1 the system is in the superfluid phase, i. e. the global U (1) symmetry is spontaneously broken. In section 3 we will keep T /T c fixed at T /T c = 0.68 (corresponding to µ = 6) and vary only the initial condition of our physical setup, that is the initial size of the vortex dipole. In section 4, on the other hand, we will keep the initial condition fixed and only vary the temperature, T /T c . All our numerical simulations are performed on a grid of 512 × 512 points in the (x, y)-plane, with periodic boundary conditions imposed. For a given T /T c we choose the grid constant, a, such that an isolated vortex is resolved by 13 grid points in diameter at 95% of the background density n 0 . For T /T c = 0.68 this yields a = 1/3.5. Identical grid constants are used in x-and y-direction such that an isolated vortex is resolved in a rotationally symmetric fashion. In section 4, where we discuss the dependence of our physical system on the temperature T /T c , we give a table with all temperatures that we consider and the corresponding grid constants. For the holographic z-direction we use 32 collocation points and accordingly a basis of 32 Chebyshev polynomials. The system is propagated in time by an explicit and adaptive Runge-Kutta time-stepping scheme.
In the following sections we often show vortex positions etc. at unit timesteps. One unit timestep (also: unit of time) is composed of an adaptive number of much smaller timesteps τ 1 that we use in the solver for the differential equations. Further details on the numerical implementation are given in appendix B.

Dynamics of Holographic Vortex Dipoles
In this section we discuss the dynamics and the annihilation process of a single vortexanti-vortex pair, or vortex dipole. Our main interest concerns the vortex trajectories and the vortices' local velocities and accelerations. We introduce a high-precision tracking algorithm for the position of the vortices that allows us to extract the trajectories in a quasi-continuous manner. It is this high precision that then enables us to accurately determine the velocity and acceleration fields of the vortices. Knowing the acceleration of the vortices we can in turn infer their (effective) interaction if we assume their (unknown) mass to be constant.
Throughout this section the temperature of the superfluid will be fixed to T /T c = 0.68.

Initial Conditions
We study a system composed of only one vortex and one anti-vortex, separated by an initial distance of u 0 grid points. This is the most basic non-trivial configuration and thus provides a clean environment for studying vortex dynamics. One possible source of numerical uncertainty originates from the periodicity of the grid and thus of the initial configuration. To minimize it, we restrict ourselves to sufficiently small initial dipole sizes.
In other words, we initialize our vortex dipole such that the 'inner' distance, i. e. the u 0 grid points in the initial condition, is much shorter than the 'outer' distance, i. e. the distance from vortex to anti-vortex obtained by crossing the boundary to the next periodic cell of the lattice. (We always assume the dipole to be placed close to the center of a cell of the periodic lattice.) With regard to this problem, we can access somewhat larger initial distances u 0 by placing the two vortices on the diagonal of the quadratic grid. We perform all numerical simulations presented in this work on a grid of 512 × 512 points in the (x, y)directions. We have checked that our results are indeed independent of the grid size by performing one simulation on a 1024 × 1024 grid which gave identical results. We conclude that for all initial dipole sizes u 0 considered here the observed dynamics is not affected by periodicity effects. When we align the vortices along the diagonal it is in general not possible to choose u 0 to be an integer. For simplicity, we then quote the closest integer as the value of u 0 . A vortex is a topological defect of the superfluid and has an intrinsic structure solely determined by the superfluid itself. The main property of a vortex is its phase winding w around the vortex position at which the superfluid density vanishes. The quantized phase winding is given by a contour integral encircling the vortex position, w = 1 2π dϕ for ϕ = arg(ψ), and is by definition positive for vortices and negative for anti-vortices.  Vortices can be built into the system by multiplying the complex scalar expectation value ψ by a localized phase, ψ → ψ · e iwφv , where φ v is the appropriate phase angle in the (x, y)-plane. In our holographic system we multiply the complex scalar bulk field Φ by a localized phase Φ(t, x, z) → Φ(t, x, z) · e iwφv at every z-slice along the holographic direction. At the position of the vortex we in addition set the scalar field Φ to zero, again for all z-values. For the boundary superfluid this corresponds to a zero in the condensate density at the vortex position. As the system is evolved in time, the topological defect imprinted in this way fully builds up within a few unit timesteps and becomes a vortex of characteristic size and shape. 4 In figure 1 we show one of the initial phase and density configurations (after it has fully built up at t = 10) used in this work. The vortex-anti-vortex pair is aligned along the diagonal of the grid, initially separated by a distance u 0 of 80 grid points. The vortices actually start moving immediately after they have been imprinted due to the phase configuration of the respective other vortex. When the vortex shape has fully built up the position of the vortices has therefore already changed and does no longer exactly correspond to the initial distance. The plot on the left shows the superfluid density n = |ψ| 2 scaled to the background condensate density n 0 (the equilibrium density of the superfluid without any topological defects). The plot on the right shows the phase of the order parameter field,

Vortex Dynamics
As the system evolves in time the vortices move on characteristic trajectories, approaching each other while simultaneously performing a motion in the perpendicular direction. We denote the motion along the dipole axis as longitudinal motion, and the one perpendicular to it as transverse motion. In figure 2 we show four snapshots of the superfluid density during the dipole evolution for an exemplary initial separation of u 0 = 80 grid points at different times, from upper left to lower right with proceeding time. It takes t = 1551 unit timesteps for the vortices to annihilate for the chosen initial condition. The first snapshot (a) is taken at time t = 400, i. e. early in the evolution when the vortices are far apart and still perfectly spherical. At this time the width of the vortices is much smaller than the vortex separation and the vortices can thus be interpreted as pointlike defects in the superfluid density. We will come back to this in section 3.4.2 where we compare the trajectories to equations of motion in the point-vortex approximation. Snapshot (b) is taken at time t = 1530, already shortly before the annihilation. At this stage, the superfluid depletions of the vortices start to deform non-elliptically and subsequently merge. This merging process begins only at very late time, here ∼ 25 timesteps before the final annihilation, and proceeds very quickly. The transition from non-overlapping to overlapping density depletions is of course not sharp, but one can still reasonably well estimate the time when it starts. The transition will be relevant for the precise tracking of the vortices, see section 3.3 below. In the following, we will refer to the last 25 timesteps as late times. In snapshot (c), four unit timesteps before the annihilation, the vortices have merged into a peanut-shaped system with strongly overlapping density depletions. At this time, the dipole accelerates most strongly. After the annihilation, we observe a soliton-like excitation propagating along the transverse direction of the vortex-dipole. This is shown in snapshot (d). The soliton-like excitation transitions into sound waves. Each annihilation event of a vortex-anti-vortex system therefore releases energy into the system in the form of sound waves. The sound waves eventually dissipate completely, leaving the system in a homogeneous equilibrium state.

Tracking Algorithm
So far, in all real-time simulations of holographic superfluids vortices were tracked by locating their phase windings on the numerical grid, see e. g. [51,[53][54][55][56][59][60][61][62]. The same holds for many studies of superfluids using the Gross-Pitaevskii equation (GPE) 5 , see for instance [73][74][75]. This procedure yields the vortex positions with plaquette resolution which is sufficient for many applications. We find, however, that to study and resolve the full dynamics of a vortex dipole a sub-grid-spacing resolution is required. Only then can velocities and accelerations be accurately computed. There exists a number of tracking methods in the literature, see e. g. [76][77][78][79][80][81] for applications of different methods, each with different advantages and drawbacks. For the purpose of this work we propose a new tracking algorithm for vortices on a two-dimensional grid and combine it with an existing Newton-Raphson (NR) tracking method, see for example [82]. This ensures that at each stage of the evolution the vortices are most accurately and numerically efficiently tracked.

Gaussian Fitting Routine
In the following, we introduce and discuss a new tracking routine that has the advantage of being numerically very efficient and applicable for most of the evolution of the vortex dipole. At a given time t of the evolution we consider the condensate density n(t, x) = |ψ(t, x)| 2 in a subregion 6 that encloses both vortices. We then use a fitting routine that locates the minima of the density depletions, i. e. the vortex positions, based on an interpolation between grid points. This permits a determination of vortex positions with sub-plaquette precision. As fitting function we use a linear combination of two two-dimensional upsidedown Gaussians, where (x 1 , y 1 ) and (x 2 , y 2 ) are the (x, y)-positions of the vortex and anti-vortex, respectively. The positions are the fitting parameters we are mainly interested in. The widths σ x/y,i/j of the Gaussians in x-and y-direction are also separate fitting parameters to allow for non-spherical deformations of the vortices. Finally, A, B and C are additional parameters in the fit, accounting for the background density and a scaling of the Gaussians. For the fitting routine we use a Levenberg-Marquardt least-squares fitting algorithm. The reason that this fitting routine accurately tracks the vortex positions is the following. Essentially, the method relies on the fact that an extremum of a function of simple and approximately known shape can be found with high precision even when its values are known only on a much coarser discrete lattice. Vortices are characterized by their quantized phase windings which implies that the superfluid density vanishes at their exact location. Consequently, to track the vortices one has to find these minima. For spherical or elliptically deformed vortices the position of the minimum of the density depletions is fully encoded in the flanks. Thus, one only has to accurately describe the flanks (in our method by Gaussians) and interpolate to the minimum. The minimum of the Gaussian then coincides with the vortex location. A Gaussian turns out to be a very good functional approximation in the region of the flanks and can thus be used to determine the vortex positions. We point out that in order to most accurately describe the flanks, the Gaussian fit should not be forced to reach zero density since it does not capture the full vortex shape.
Notably, any other function that describes the flanks of the vortices similarly well could also be used. We have checked this explicitly and found very good agreement between the results. For instance, another choice could be the approximate shape [83,84] of a vortex from the Gross-Pitaevskii equation, given by where n(r) is the density profile with r denoting the radial distance from the vortex center, n 0 the background density, and ξ a free parameter determining the width of the vortex. Using this function in the fitting routine yields a discrepancy with the Gaussian fit of approximately ±0.02 grid points in the vortex positions. We estimate that the error of the vortex position resulting from the Gaussian fits is of the same order of magnitude. This can be reasoned by testing how well a GPE vortex (cf. the same estimate. In particular, we compared the results with the NR method discussed below and found excellent agreement. In all plots shown in this work the error in the vortex positions is much smaller than the size of the plot markers. Figure 3 shows the cross section of a two-dimensional Gaussian fit to a vortex-antivortex system. Instead of showing the full dipole system we representatively focus on just a single vortex but stress that we in fact did a combined fit to both vortex and anti-vortex as described above. The left plot in figure 3 shows the cross section of the two-dimensional Gaussian fit to the superfluid density n(x, y) along the x-direction for fixed y, with the y-position taken as the closest integer y-grid point to the position from the fit. Crucially, the Gaussian fit does not reach zero but captures the vortex flanks to high precision which is sufficient to find the position of the minimum. On the right in figure 3 a contour plot of the same fitted vortex profile is shown. The contour lines of the fit are displayed in grayscale while for the condensate density n we use the same color coding as in figure 1.
This Gaussian fitting algorithm can be used to track the vortices as soon as their shapes have fully developed after initialization, that is after 5 unit timesteps. Only for the last 25 unit timesteps, when the vortices merge and deform non-elliptically, this method breaks down and thus fails to track the vortices as the Gaussians can non longer accurately capture the vortex flanks. For this regime we instead rely on a well-known Newton-Raphson tracking algorithm which we discuss next.

Newton-Raphson Tracking Algorithm
A reliable way to track the vortices also in the late-time regime when the density depletions of the vortices overlap is given by the NR method. As a standard algorithm for numerical root finding, the NR method is optimally suited for vortex tracking on two-dimensional grids, see for example [77,78,85]. In the following we give a brief overview of the method and discuss its advantages and drawbacks as compared to the Gaussian fitting routine. As discussed above, vortices are characterized by ψ(x i ) = 0, with x i the position of the ith vortex, which is a direct consequence of the phase winding. The NR method can therefore be used to solve ψ(x i ) = 0 directly, and no assumptions regarding the shape of the density depletion are needed. To ensure that a zero at x i indeed signals a vortex one subsequently has to check that there is a corresponding phase winding at this position. To solve ψ(x i ) = 0 for x i using the NR method, one performs a Taylor expansion of ψ(x) around an initial guess x g i for the vortex position, with the Jacobian matrix has full rank, (3.4) can be solved for x i , This procedure can then be iterated such that x i approaches the exact vortex position to sub-grid-spacing precision. The precision is mainly limited by the evaluation of the Jacobian at intermesh points. This is done by Fourier interpolating ψ(x) which also guarantees that the periodicity of the order parameter field is retained throughout the iteration procedure. As initial guess for the algorithm we locate the phase winding of the vortices on the numerical grid. This plaquette precision is sufficient to ensure convergence of the NR algorithm. We stop the iteration when the condition ψ(x i ) < 10 −10 n 0 is satisfied. The advantage of this algorithm as compared to the Gaussian fitting routine is that it directly tracks the zeros in the superfluid density. The Gaussian fitting routine, on the other hand, determines the minima of the density depletions only via their flanks. That is very precise as long as the vortices deform only mildly and elliptically but fails as soon as this is no longer the case. For the vortex dynamics under investigation this happens as part of the final annihilation process where the vortices deform strongly. For all earlier times, when the vortices are spherically symmetric or only slightly elliptically deformed, we checked that the two methods agree to a very high precision with a maximum absolute deviation of ±0.02 grid points. This error agrees with the estimate made above for the Gaussian fitting routine. The key advantages of the Gaussian fitting routine are its simplicity and its numerical efficiency which make it much faster than the NR method. Therefore, we use the Gaussian fitting algorithm to track the vortices for all vortex-anti-vortex separations d(t) larger than 25 grid points and the NR method for all smaller separations. As stated above, the two methods yield the same results for the first stage, i. can no longer accurately track the vortices. The trajectories obtained with this method cannot resolve the final annihilation process and have a systematic error, as is evident from the almost parallel final pieces of the trajectories which would, if real, not lead to an annihilation of the vortices. The NR method, in contrast, provides the accurate vortex positions until immediately before the annihilation.

Sub-plaquette Vortex Tracking
Having discussed the tracking algorithms, we now focus on the gain of these methods as compared to the standard plaquette resolution obtained from locating the phase windings on the grid. The position found by locating the phase winding is inherently plagued with an uncertainty of one grid point as the vortex position is always assigned to integer grid points. Consequently the trajectories determined by this method take a stair-like shape, with the vortex position being allocated to one grid point for a certain number of timesteps before jumping to the next one. and with sub-plaquette resolution using our localization method (blue dots). The phasewinding localization is only as precise as one grid point and the trajectory is therefore discrete or stair-like, while our localization method results in quasi-continuous trajectories. resolution). In figure 5 we compare the vortex trajectories determined by our (combined) tracking method to the trajectories determined by using only the phase windings to locate the vortices. Clearly, the gain from the combined use of Gaussian fitting routine and NR method has a strong impact on the extracted trajectories. We note that there is a systematic contribution to the deviation between the two methods coming from the fact that, when locating the vortices via their phase windings on a plaquette, we always attribute the position to the upper left corner of that plaquette. Any other choice would imply a similar systematic deviation that would be stronger in other parts of a given trajectory. Let us give an empirical estimate of the difference in precision between the two methods for our specific vortex-dipole system. At early and intermediate times (to be discussed in more detail below) it takes roughly 50 unit timesteps for the vortex to move as far as one grid point. Locating the vortex solely by its phase winding would thus give the same grid point for all of these 50 unit timesteps. By using our fitting routine, on the other hand, we can clearly locate the vortex at 50 distinct positions. Therefore, our method yields quasi-continuous trajectories for the vortices. As compared to the locations obtained from the phase windings this is a strong increase in precision. Obviously, the gain in precision depends on the velocity of the vortices and, in relative terms, is smaller for fast vortex motion. We will see that the increase in precision due to sub-plaquette tracking is essential for obtaining some of the main results in our analysis of the vortex dipole. For obtaining the data shown in figure 5, as well as for all other results presented in this work, we choose t = 5 as the starting time for our fitting procedure. 7 The trajectories are shown up to the point where the vortices annihilate. Since the vortices accelerate strongly just before annihilation, we increase the frequency of tracking during the last 8 unit timesteps and determine the position 100 times per unit timestep here.

Trajectories of the Dipole System
Let us now investigate in detail the dynamics of a vortex-anti-vortex pair. We first present and discuss the trajectories of the vortex dipole, here for initial separations along the horizontal of the quadratic grid. The results for five different initial separations are shown in figure 6. Note that for conciseness the plot in figure 6 does not show the full numerical grid but only the part containing the actual vortex trajectories. We observe that the trajectories of the vortices at intermediate and late times are independent of the initial separation and follow universal curves. That is, the short-distance behavior is independent of the earlier evolution of the system. We checked this explicitly by shifting the trajectories such that the respective annihilations take place at the same point on the grid. Figure 6 shows the un-shifted vortex trajectories. We refer to all times from which on the trajectories are universal, but prior to the last 25 unit timesteps, as intermediate times. We stress that the distinction between early, intermediate and late times in this sense is only approximate and does not split the evolution into three strictly distinguishable phases. The universality of the trajectories at intermediate and late times is a remarkable feature of the system as it shows that the physics of a vortex-anti-vortex pair in the superfluid is entirely governed by the size of the dipole and independent of its previous evolution. In subsection 3.6 we will discuss this point in more detail in the context of the vortex velocities and accelerations.

Phase Healing
In figure 6 we observe a feature of the trajectories at early times that is, at first sight, very surprising and gives rise to the non-universality of the trajectories in this regime. Both the vortex and the anti-vortex show a slight outward bending in their motion before they approach each other. The effect is discernible only when the vortices are tracked with subplaquette resolution. This is obvious in figure 5 where the quasi-continuous trajectories (blue dots) exhibit a clear outward bending while the discrete trajectories (red dots) limited to plaquette resolution do not show any sign of outward bending. In figure 6, the outward bending is seen in the early-time evolution for all initial separations. It turns out that the effect is not physical but is an artifact of the method used for imprinting the initial configuration of vortices in the superfluid order parameter field. As discussed above, we initialize two single vortices by simply imprinting and linearly superimposing their individual phase configurations on the entire grid. But this does not take into account that a proper vortex-anti-vortex pair is actually one coupled system. Imprinting such a pair would require us to imprint the coupled phase configuration instead of the linear superposition of two separate vortices. However, an analytic formula for the correct phase configuration of a proper vortex dipole is not known. Imprinting the superposition of the phase fields of two separate vortices is therefore the best practical procedure. Fortunately, during the numerical evolution of the system the phase configuration 'heals itself' and transitions into that of a proper vortex dipole. This 'phase healing' process is much slower than the build-up process of the vortex shapes, as can be expected because it is necessarily a non-local process. After the phase healing the dynamics of the dipole becomes independent of its initial configuration and is then solely determined by the distance between vortex and anti-vortex.
We refer to 'intermediate times' of the evolution for times after the phase healing is completed, while we call times during the phase healing process 'early times'. Equivalently, intermediate times refer to the universal regime of the evolution which is independent of the initial separation of the vortex-anti-vortex pair. We find that the phase healing process is mildly dependent on the initial condition, i. e. the initial separation of the vortices, cf. figure 6. However, the effect of the phase healing on the trajectories is rather subtle and fades out slowly. It is therefore difficult to determine the exact time of the transition from the phase-healing regime to the universal regime. The analysis in section 3.6 will show that the vortex velocities permit a reasonably accurate estimate.
For an initial separation of u 0 = 80 grid points we find that it takes about ∆t ∼ 350 unit timesteps for the phase configuration to approach that of a proper vortex dipole. We can understand how the phase healing leads to a slightly outward pointing motion by comparing the phase configurations of two dipoles of identical size, with one of them being in its initial phase configuration and the other one being in a phase configuration after the healing process. The latter can be obtained by evolving a larger dipole through the phase healing. The difference between the two corresponding phase configurations gives the initially 'missing configuration' (as compared to the fully coupled dipole configuration), and we can derive the superfluid flow resulting from it.
To visualize this effect, we initialize a vortex pair as a superposition of the individual phase fields at large distance and propagate it in time for long enough so that its phase field has fully healed. From the corresponding phase configuration we subtract the phase configuration we would have imprinted for a new pair of vortices at this very separation. This gives us the missing phase configuration of the newly initialized dipole as compared to a fully coupled dipole. Its negative gradient is the superfluid flow that acts on the newly imprinted vortices at this separation due to their missing phase configuration. We show the streamlines of this flow field superimposed on the superfluid density of the dipole configuration in figure 7. Here we have chosen the dipole orientation along the diagonal of the grid so that the vortex dipole motion is pointed towards the lower-right corner of the plot. We can clearly see the outward pointing streamlines that lead to the observed outward bending of the trajectories. Furthermore, we observe streamlines pointing in the direction opposite to the vortex motion. In subsection 3.6 we will see that this leads to an initial phase of deceleration in the vortex motion. Understanding the artifacts of the initialization procedure allows us to separate them from actual dynamics of the vortex dipole by concentrating on the universal regime of the evolution.
It is worth pointing out that related effects due to coupled phase fields have been observed in the context of kinks in scalar field theory [86,87].

Comparison to Hall-Vinen-Iordanskii Equations
There exist well-known equations of motion, first introduced by Hall, Vinen and Iordanskii (HVI) [88,89], for the dynamics of a number of vortices in the point-vortex approximation. They can be derived on very general grounds in a hydrodynamical setup, taking into account the interactions of the vortices with the normal and the condensate components of the fluid. The forces that drive the vortex dynamics are the Magnus force, orthogonal to the vortex motion, and drag forces. The equations of motion are derived by setting the sum of all forces in the closed system to zero and solving for the vortex velocities [90,91]. Applied to our setup, there is no velocity of the normal component of the fluid, and each vortex only feels the velocity field created by the respective other vortex. The coupled differential equations for the two-dimensional trajectories x i , where i ∈ {1, 2} labels the vortices, are then given by Here, C and C are phenomenological friction coefficients characterizing the mutual friction between the vortices and the fluid, w i denotes the winding number of the ith vortex, and v s i is the superfluid velocity induced by the respective other vortex. We further denote by e ⊥ a unit vector perpendicular to the (x, y)-plane in a right-handed coordinate system, pointing in an auxiliary direction used only for a convenient notation of this equation and not to be confused with the holographic direction.
For a given set of initial conditions x i (0) = x 0,i , equation (3.7) can straightforwardly be solved. For one vortex (w i = 1) and one anti-vortex (w i = −1) we give a derivation of these solutions in appendix C. It turns out that the vortices move along straight lines and annihilate under a fixed angle, depending on the parameters of the theory.
As we have seen, vortices in the holographic superfluid have a finite extension and in general cannot be considered as point-like defects. However, at large vortex-anti-vortex separation as compared to the vortex width, a point-vortex approximation is possible as far as the effective interaction between the vortices is concerned. We can therefore compare the analytic results from the HVI equations to the vortex trajectories discussed above. Figure 6 clearly shows that the vortices do not follow straight lines for most of the evolution. We note in particular that at early times, when the vortices are furthest apart, the motion is affected by the phase healing and thus not expected to obey the equations of motion for a proper dipole. Furthermore, we find that after the phase healing the vortices quickly accelerate and approach each other, again rendering the point particle approximation invalid. This is clearly visible in the vortex trajectories as the final annihilation process with overlapping density depletions is preceded by a strong bending of the trajectories which sets in not long after the phase healing. Only in between these two regimes a possible straight segment of the trajectories may be found. It turns out that for most of the initial separations we study here, this segment, if it exists, is too short to be properly identified. For an initial separation of u 0 = 120, however, we indeed find that the vortices follow straight lines for a short period of time after the phase healing and prior to the strong bending and ensuing annihilation. We present more details and a plot illustrating this behavior in appendix C. The phase healing ultimately is an artifact of our method of preparing the initial state. Thus we find that for sufficiently large separation the fully developed holographic dipole system indeed obeys the HVI equations of motion.
Comparing in detail the evolution of a vortex dipole in the holographic superfluid, in the Gross-Pitaevskii framework and in the HVI equations, using the tracking methods developed in this work, offers even the possibility to determine the characteristic HVI friction coefficients of the holographic superfluid, see [92].

Size of the Dipole System
Next, we study how the vortex-anti-vortex separation behaves as a function of time. The separation, or size of the vortex dipole, in units of grid points is defined as where x i (t) denotes the position of vortex i ∈ {1, 2} on the numerical grid. In physical units, the distance is given by d(t) = u(t)a, where a denotes the grid constant as discussed above. (We recall that in this section we consider the temperature T /T c = 0.68 for which a grid constant of a = 1/3.5 is chosen.) In figure 8 we show the physical distance as a function of time, again for an initial separation of u 0 = 80 grid points. For all other initial separations discussed in this work d(t) has the same functional behavior, as we discuss in more detail in appendix D. We expect the dipole size at early times to be affected by the phase healing process and thus not to reflect the correct behavior of a proper vortex dipole. Therefore we focus on intermediate and late times in the following. There we find that d(t) follows a power law for all initial separations. Specifically, we perform a fit to our data with the ansatz where t 0 , A and b are free fit parameters, and find that excellent agreement is achieved. According to (3.9) one would expect t 0 to be the time at which the vortices annihilate. However, during the final stages of the annihilation process, the behavior of the vortices is dominated by the deforming and merging density depletions and is not adequately described by the notion of moving vortices in the sense of a formula like (3.9). We therefore exclude the final 5 unit timesteps from the analysis and leave t 0 as a free parameter of the fit. For all initial separations we find for the exponent b = 0.53 ± 0.02. We present more details of the fit and further discussion of the results in appendix D. We compare this result to the behavior expected from the HVI equations (3.7) for point vortices. As discussed in appendix C, the HVI equations predict the vortex-anti-vortex separation d(t) to follow a power law of the form (3.9) with b = 0.5, see equation (C.10). This is very close to what we find numerically for the holographic vortex dipole. However, we note that the holographic vortices in general cannot be considered as point-like, and deviations from the square root behavior are therefore expected. Nevertheless, our results show that also for the case of finite-width vortices their spatial separation obeys a universal power law after the initial phase healing.
The dynamics of a vortex dipole and its annihilation in the holographic superfluid were also addressed in the recent study [59], including in particular a calculation of the time dependence of the dipole size. In the following we comment on how our results compare to those of [59]. We first note that we study exactly the same system at the same temperature. The methods for the numerical simulation are similar, with Fourier decomposition in the (x, y)-plane and a pseudo-spectral method with an expansion in Chebyshev polynomials in the holographic z-direction. A grid of 241×241×25 points in the (x, y, z)-directions is used in [59] while we use a larger domain of 512 × 512 × 32 points. A key difference between our work and [59] is that we use a high-resolution method for tracking vortices while the analysis of [59] relies on tracking by phase winding, thus providing only plaquette resolution. The authors of [59] observe that d(t) splits into two regimes each of which is described by a power law, with a transition at a separation of twice the vortex radius (defined as the radius where 99% of the background density are reached). At large separations the square root behavior predicted by the HVI equations is found, while an exponent of 2/5 is found at small separations. Velocities, accelerations and effective interaction strengths of the vortices are then derived from these power laws, in contrast to our calculation in which velocities and accelerations are computed directly. With improved precision in vortex tracking, we find different results for d(t). For large distances, we find a similar power law, although with a slightly different exponent. For small distances, we do not observe a second, distinct power law. We further observe that the point-vortex approximation breaks down at distances much larger than twice the vortex radius. The likely reason for the different results is the higher precision of our study, but also the larger numerical grid might play a role.

Velocities and Accelerations
The high-precision tracking of the vortices allows us to directly compute their velocities and accelerations. We can use the quasi-continuous trajectories determined above to compute derivatives via finite differences, as we describe in appendix E. Furthermore, we can then determine the dependence of velocities and accelerations on the dipole size which is interesting with regard to the observed universality of the dipole motion.
The results of this section again rely heavily on the improved tracking procedure, in particular for the slow motion of vortices at early and intermediate times. Here, a vortex needs many unit timesteps to traverse a single grid spacing. Accordingly, the velocity calculated from tracking with only plaquette resolution would be zero for most of the time, except for occasional jumps when the next grid point is reached. From quasi-continuous trajectories, on the other hand, we can compute quasi-continuous velocities and accelerations.

Vortex Velocities
In this section we consider physical distances (and velocities), i. e. we use d(t) = au(t), with a the grid constant and u(t) the dipole size in grid point. Consequently, velocities are expressed in units of the speed of light. To ease comparison with the trajectories discussed in the previous sections, we quote the initial separation of the dipole u 0 in units of grid points. We recall that we can split the dipole motion into longitudinal and transverse parts. Analogously, we denote the velocity component along the dipole axis as longitudinal velocity v , and the component perpendicular to the dipole axis as transverse velocity v ⊥ . In figure 9 we show the two velocity components as function of the dipole size for all five different initial separations we considered before. To wit, in these plots time proceeds from right to left as the dipole shrinks with increasing time. Figure 9 shows that after some time the velocities for all different initial separations approach a universal curve. This holds for both the longitudinal and the transverse component. The deviations at early times are caused by the initial phase healing discussed above. At intermediate and late times, the velocities are universal and depend only on the size of the system, confirming the universality observed in the trajectories in 5. We also find that for most of the evolution the total velocity of the dipole system stays well below the speed of light, such that the vortex motion can be considered non-relativistic. This confirms the finding of [51] that the boundary dynamics of the holographic superfluid is non-relativistic although the holographic model has an inherent Poincaré invariance. So far we have only considered the evolution of vortex dipoles of different initial separations. To test the universality of the vortex motion at intermediate and late times, we also performed simulations with a different type of initial condition. For that we imprint the vortex and anti-vortex as before, but in addition perturb the system with random noise 8 at time t = 0. Due to that, the vortices acquire a different initial velocity both in absolute value and in direction. As a result, the trajectories at early times are very different from what we observed in figure 6. At intermediate and late times, however, they can be shown to agree again after an appropriate rotation of the coordinate system (with an angle depending on the initial motion due to the respective noise). Analogously, we find excellent agreement of the velocities and accelerations after the phase healing and after the initial noise has dissipated. This again confirms the universal dynamics of the vortex-dipole in which the system loses information about its initial condition. In the following we will consider again the usual initial dipole configurations without noise.
The phase healing at early times has a marked effect on the velocities. The transverse component of the velocity even decreases initially before it increases again. We note that at this stage the velocity is exclusively induced by the initialized phase configuration which, however, does not match the correct phase field of a coupled dipole. Therefore, the velocities receive a contribution according to the streamlines shown in figure 7. There, we indeed found a component pointing in the direction opposite to the center-of-mass motion of the dipole. Once in the universal regime, on the other hand, the velocities exhibit the expected behavior of a vortex dipole with a purely attractive force due to the strongly dissipative nature of the holographic superfluid in the presence of vortices.
In a superfluid condensate without dissipation, in contrast, vortices are known to exhibit so-called Helmholtz propagation in which vortex and anti-vortex move at fixed distance on straight parallel lines. The velocity of the dipole is then purely dictated by the vortex-anti-vortex distance and follows a simple inverse-distance law due to angularmomentum conservation in a closed system. As our system is dissipative, Helmholtz propagation is suppressed and we expect deviations from the inverse-distance law. To illustrate this, we plot a 1/d curve in figure 9 for the transverse component of the velocity.

Vortex Accelerations
Although the velocity for given initial positions of the vortices already contains the full information about the dipole motion, we still find it instructive to show the corresponding acceleration, which we again split it into its longitudinal and transverse components. Figure  10 shows the longitudinal component a and the modulus of the transverse component |a ⊥ | of the acceleration in physical units. The accelerations exhibit the same universal behavior at intermediate and late times after the phase healing. We note that the dips in |a ⊥ | are due to a change in sign of a ⊥ .
A comment is in order concerning the interpretation of vortex accelerations in terms of force fields. Force fields usually require an interpretation in terms of test particles to be placed at a given position. However, such an interpretation is questionable for a quantized vortex as it always has a strong back-reaction on the fluid. Formally, one could of course compute acceleration vectors and interpret them in terms of a mutual force field acting on the vortex and the anti-vortex. This would furthermore require the knowledge of an effective mass of a vortex. The latter problem has been addressed in the literature, mainly in the context of Gross-Pitaevskii theory, suggesting a definition for the mass of an isolated single vortex, see for instance [93]. Unfortunately, none of these models can be directly applied to our system and the effective vortex mass remains unknown. Nevertheless, we can infer from the accelerations in figure 10 that such an effective force depends only on the spatial separation and is not a central force.  Figure 11: Order parameter as a function of temperature (lower abscissa) or chemical potential (upper abscissa). As the temperature T falls below the critical temperature T c a condensate forms, corresponding to a non-vanishing order parameter that increases with decreasing temperature. On the chemical potential axis we mark the values of µ used in this section.

Vortex Dipole at Varying Superfluid Temperature
So far, we have considered vortex dipole dynamics only at a fixed temperature T /T c = 0.68 of the holographic superfluid. In the following we compare the vortex velocities for a given initial separation while varying T /T c . In general, we do not expect the probe approximation of the holographic superfluid to fully capture all effects of the temperature. However, given that a fully back-reacted computation has not yet been performed, we find it interesting to see how the dynamics changes with temperature. We keep in mind that the results might be different in a fully back-reacted setting.

General Considerations and Initial Conditions
We recall from section 2.1 that the holographic superfluid has only one dimensionless parameter µ/T . Setting units by z h = 1 leaves only µ as free parameter which fixes the temperature T /T c . In figure 11 we show the modulus of the order parameter field in thermal equilibrium, which serves as our background solution, as function of the temperature T /T c or the chemical potential µ. On the µ-axis we mark those values that we use in this section. Changing the temperature of the holographic superfluid implies a change in the superfluid condensate density (see figure 11), a change in the time scales of the system, and a change in the typical size of a vortex. In order to always resolve a vortex by the same number of grid points we account for the changing vortex size by adjusting the grid constant for each fixed temperature. This allows us to avoid any potential artifacts due to velocities and accelerations at different temperatures being determined at different numer-  ical resolution of the vortex structure. Again, all our numerical simulations are performed on a 512 × 512 grid in the (x, y)-direction. We note that choosing different grid constants at different temperatures implies different physical sizes of the total system. We choose five different temperatures covering a broad range of the values in which the system is numerically accessible. Above and below our chosen temperatures the numerical solution turns out to be plagued by instabilities. One of the temperatures is the value chosen above, T /T c = 0.68 in the center of the accessible temperature range. In order to compare the dynamics of the system at different temperatures, we use dimensionless units for the quantities in this section. To account for the dependence of the vortex size on the temperature T /T c we study the vortex velocities for a given T /T c as a function of the dimensionless distance d(t)/ξ, where ξ is the vortex width as determined by fitting (3.3) to the vortex profile for the given temperature. Choosing instead the same initial separation in grid points would imply that the vortices are well separated for some temperature while they would already start to overlap for another. We recall that for given temperature we choose the grid constant a such that an isolated vortex is resolved by 13 grid points in diameter at 95% of the background density n 0 . For the initial separation we adapt u 0 such that d 0 /ξ (with d 0 = au 0 ) it the same for the different temperatures.
In table 1 we summarize our choices for the temperature, the corresponding chemical potentials, and the resulting background densities n 0 and vortex width ξ. The initial separation d 0 is chosen as 39.2 ξ. In addition, we show the corresponding size L of the computational domain in units of ξ.

Velocities of the Vortex Dipole
With these parameters we evolve the systems in time and track the vortices until they annihilate. We find that the initial time it takes the vortices to fully build up their shapes after we imprint them is slightly dependent on the temperature but is always shorter or equal to ∆t = 12 unit timesteps which is well within the regime of the phase healing. As the velocities in this regime are anyway affected by artifacts of the vortex initialization we can disregard any small differences in the build-up process. We start the tracking for all chosen temperatures after the vortices have fully developed. We point out that the tracking based on Gaussian fitting of the vortices works equally well for all temperatures since the vortices are always resolved by the same number of grid points. The longitudinal and transverse velocities, defined as before and again computed via finite differences from the tracked trajectories, are shown in figure 12. The figure shows that the effect of the phase healing is present for all temperatures, but more pronounced at higher temperatures. We further see that after the phase healing process the velocities exhibit a clear ordering according with temperature. Both transverse and longitudinal velocities decrease with increasing temperature.
An interpretation of the ordering of the velocities with temperature in terms of dissipation, for example along the lines of Tisza's two-fluid model, would appear speculative as the results are obtained in the probe-approximation and might be affected when back-reaction is included.

Summary and Outlook
We have studied the dynamics of a single vortex-anti-vortex pair in a holographic superfluid in numerical real-time simulations. A new tracking method has been introduced which allows us to locate the vortices with strongly improved precision as compared to the conventional plaquette-based technique of locating the phase windings in the condensate density. The new tracking method relies on locating the minimum of the density depletion of the vortex rather than only its phase winding. Using a fit with Gaussians, supplemented by a Newton-Raphson method at small vortex separations, we can determine the vortex positions with sub-plaquette resolution in a numerically efficient manner. Fitting both vortices simultaneously with a linear combination of two two-dimensional Gaussians even allows us to capture the mutual deformation at small distances, thus eliminating a further source of error in the vortex positions. Our tracking method can be applied throughout the evolution of the vortex dipole. With this method we have, for the first time for the holographic superfluid, computed quasi-continuous trajectories of the vortices during their approach and eventual annihilation. In addition, we use large numerical grids of 512 × 512 points in the (x, y)-directions such that finite-size (or periodicity) effects are efficiently suppressed, as we have checked on a 1024 × 1024 grid.
The increased precision makes it possible to study many aspects of the dynamics of the quantized vortices which are not accessible with plaquette-resolution. In particular, we are able to derive directly the velocities and accelerations of the vortices during the whole evolution of the dipole. This allows us to improve earlier results of the pioneering study [59] in which the velocities were derived from plaquette-resolution trajectories. We partially confirm the results found there for the power law for the velocity during the dipole evolution. We find a very similar though not exactly equal exponent for large vortex separations, with our result indicating a slight deviation from the behavior expected from the Hall-Vinen-Iordanskii equations for point-like vortices. However, we do not confirm the finding of a second power law for small distances reported in [59].
A key result of our study is the observation of a universal behavior of the vortex dipole in the holographic superfluid. Specifically, we find that its trajectory and velocities are universal independent of the initial separation. Hence their interaction in the superfluid is only determined by their separation and does not depend on the previous motion. This finding holds for all fixed temperatures of the superfluid condensate. We point out that we observe in our simulations some deviations from this universal behavior at early times, visible as a slight initial outward bending of the trajectories. But we have shown that these can be fully attributed to artifacts of the numerical initialization of the phase field of the vortices. Such artifacts occur not only in the holographic superfluid but also in other numerical simulations of superfluids like for example in Gross-Pitaevskii dynamics. It might be interesting to study this problem in more detail so that one could potentially identify better methods for initializing ensembles of vortices in superfluid simulations. We note that our finding of universality in the vortex dynamics could again only be obtained due to the improved precision of the vortex tracking.
We have investigated how the dynamics of the vortex dipole depends on the temperature of the superfluid. Our results indicate an overall decrease in the vortex velocities with increasing temperature if we compare the system for suitably rescaled dimensionless quantities. However, the probe approximation for the holographic superfluid might not be well suited to obtain a full picture of the temperature dependence. It would therefore be interesting to study this question in a fully back-reacted system.
The vortex-anti-vortex system studied here is one of the most basic non-trivial examples of topological defects in a superfluid, and in some sense a basic building block for understanding more complicated ensembles of vortices. Such systems of vortices and antivortices, as well as solitons, have been widely studied for the holographic superfluid. We expect that the improved method for vortex tracking can be very advantageous also for improving our understanding of such systems. It would for example be interesting to see whether the universality of the dynamics of a vortex dipole has an equivalent for systems of many vortices.
The tracking method presented here can also be used for simulations of vortices in Gross-Pitaevskii dynamics. An interesting result of such a study is the quantitative determination of the strong dissipation and of the friction coefficients of vortices in the holographic superfluid [92]. This result has been obtained by a high-precision comparison of vortex trajectories in the holographic superfluid with those of Gross-Pitaevskii dynamics and of the Hall-Vinen-Iordanskii equations. Finally, we expect that our tracking method will also be useful in the study of vortex dynamics in three-dimensional superfluids, including holographic superfluids.
where L matter is the matter Lagrangian, We work in the probe limit and thus neglect the back-reaction of the matter part on the gravity background. The resulting equations of motion are Maxwell's equation and a Klein-Gordon equation for the scalar field, with the current The remaining gauge freedom is fixed by choosing the gauge A z = 0. In our units with L AdS = 1 we fix the squared mass to m 2 = −2. We also fix the position of the black-hole horizon to z h = 1 which leaves us with only one free parameter of the theory, the chemical potential µ. This can be tuned such that the system is in the superfluid phase.
To compute the dynamics of the holographic superfluid we first have to calculate the static and homogeneous background solution. Subsequently we imprint the vortices into this background solution and evolve the system in time. The static background solution is independent of time and of the spatial coordinates x and y of the superfluid but depends non-trivially on the holographic coordinate z. Plugging the AdS-Schwarzschild metric (2.4) into the Maxwell equations (A.3) we obtain where a prime denotes the derivative with respect to the holographic coordinate z and the index i runs over the spatial coordinates x and y.
with A = (A x , A y ). These equations are ordinary differential equations of second order in z. We already discussed the boundary conditions for the fields in section 2.1 but repeat them here for completeness. The chemical potential is given by the boundary value of the temporal component of the gauge field, i. e.
A t (z) = µ + O(z) . (A.10) At the black-hole horizon, z = z h = 1, A t has to be regular, A t (z = z h ) = 0. For the spatial components of the gauge field the boundary conditions at the conformal boundary z = 0 of the AdS spacetime are given by while at the black-hole horizon As discussed in section 2.1, the scalar field has to be regular at the black-hole horizon, which is ensured by the choice of infalling Eddington-Finkelstein coordinates. At the conformal boundary the scalar field Φ has to satisfy which sets the source of the scalar operator in the boundary theory to zero such that the U (1) symmetry is not explicitly broken. After imprinting vortices into the background solution, the system is no longer homogeneous and we have to solve the full set of equations of motion to obtain the real-time dynamics of the superfluid. For the numerical implementation it turns out to be convenient [51] to derive the full set of equations of motion in terms of the rescaled scalar field Φ = Φ/z and to use the 'lightcone derivative', where ∇ = (∂ x , ∂ y ). The full system of equations of motion is then given by The dynamical equation for A z is again used as a constraint equation and is not independent of equations (A.15)-(A.18).

B Numerical Implementation
Here we discuss the numerical implementation of our solver for the system of differential equations given in appendix A. For the holographic z-direction we employ a collocation method with 32 Gauß-Lobatto grid points and expand the fields in a basis of 32 Chebychev polynomials. After setting z h = 1, we have to use the rescaled coordinatez ∈ [−1, 1], defined implicitly via to implementz-derivatives via simple matrix multiplication. After using a Newton-Kantorovich algorithm to linearize the equations of motion for the static and homogeneous background, we solve them with an LU decomposition. The set of equations of motion for the time evolution (A.15)-(A.18) is solved on a quadratic grid with 512 × 512 grid points in the (x -y)-plane and periodic boundary conditions. Derivatives with respect to x and y are implemented via Fourier transforms. The equations are then solved as follows. At every unit of time, we compute the right hand sides of (A.15)-(A.18). At each (x, y)-position we then integrate them with respect to z, implementing the boundary conditions (A.10), (A.11) and (A.13) at z = 0. This allows us to undo the shifts of the lightcone derivative (A.14) which yields ∂ t X, for X again given by X = A x , A y ,Φ. These fields can then be propagated in time using a fourth-order/fifthorder Runge-Kutta-Fehlberg algorithm with adaptive timesteps. The adaptive timesteps vary between 0.001 and 0.1 in units defined by z h = 1. One unit timestep corresponds to ∆t = 1, and is therefore composed of a varying number of numerical timesteps. Knowing the fields A x , A y ,Φ at the new time we can plug them into equation (A.15) and solve the boundary value problem to also obtain A t at the new time. FromΦ we can then derive ψ according to equation (2.12).

C Vortex Trajectories from the Hall-Vinen-Iordanskii Equations
The Hall-Vinen-Iordanskii equations [88,89] for the velocities of a system of point vortices for vanishing velocity of the normal fluid component are given by where C, C are phenomenological friction coefficients and v s i is the superfluid velocity at the position of the ith vortex of winding number w i created by all other vortices. We use the unit vector e ⊥ perpendicular, in a right-handed coordinate system, to the (x, y)-plane.
(We point out that e ⊥ is not the holographic direction but a third spatial direction only used to write the HVI equations in a concise way.) All external forces are set to zero. The superfluid velocity at x i is given by with winding number w i of the ith vortex. We further use I to denote the unit matrix in twodimensions and introduce J = 0 −1 1 0 . Then equation (C.1) can be solved straightforwardly for the case of one vortex (w 1 = 1) and one anti-vortex (w 2 = −1). Taking the gradient of the Kirchhoff-Onsager functional we obtain Inserting this into (C.1) leads to with i ∈ {1, 2}, and we have used w i w 3−i = −1. We define the vortex-anti-vortex separation d = x 1 − x 2 with d = de d and the center of mass D = 1 2 (x 1 + x 2 ). For the separation we find dd dt = − 4 d 2 From this the solution for x i can be obtained by changing back coordinates. Irrespective of the time dependence of the dipole size we can infer directly from the HVI equation (C.1) that, in the point-vortex limit and for vanishing velocity of the normal component, the vortices trace out straight lines. In short, one expects straight spatial trajectories of the vortices when the conditions for the applicability of the HVI equations are met.
As we discussed in the main text, a comparison of the holographic vortex trajectories to the HVI equations has two main caveats. First, there is the initial phase healing (see section 3.4.1) which renders the dynamics unphysical in the sense that the system does not exactly follow the expected behavior of a vortex dipole system for some time after initialization of the two vortices. Second, the point-particle approximation of the HVI equation is no longer valid when the vortices come too close to each other. Only for an initial separation of u 0 = 120 grid points (or larger) we can clearly identify a segment of the trajectories for which the vortices follow straight lines, as shown in figure 13. For all other initial separations discussed in the main text the two regimes excluding the application of the HVI equations seem to overlap and thus prohibit the characteristic behavior to emerge. We checked explicitly that for initial separations larger than u 0 = 120 grid points (but small enough for finite-size effects to be negligible) the interval for which the trajectories are straight lines increases. This is of course expected due to the universality of the dynamics.

D Dipole Size Power Law
According to the HVI equations discussed in appendix C, the vortex-anti-vortex separation for point-like vortices in a superfluid decreases with a power law. In the following we check if such a behavior is also observed for the holographic superfluid. We again exclude from the analysis the early regime in which the phase of the system has not transitioned (or healed) into that of a proper vortex dipole. The phase healing is a dynamical process during the evolution of the system which is continuous and relatively slow, so that it is difficult to determine its exact endpoint. However, using the universality of the dipole dynamics as a criterion, it is possible to estimate reasonably well the time after which the system behaves like a proper dipole. For that, we can read off from figure 9 the time when the vortices of a given initial separation follow the universal (velocity) curve. This time is a good estimate, with an error of about ∆t = ±30 timesteps, for the endpoint of the phase healing.
For all times after that endpoint of phase healing we analyze whether the separation d(t) exhibits power-law behavior. We make the ansatz for the functional dependence, where A, t 0 and b are free parameters, and test how well it describes d(t) in the considered interval. In practice we use a Levenberg-Marquardt least-squares fitting algorithm to find the free parameters and an error estimate. For the fitting routine we ensure that all points are weighted equally. We find that d(t) is indeed well described by the power law (D.1) at intermediate and late times. For all initial separations we find the same exponent, given by b = 0.53 ± 0.02, but the parameters A and t 0 differ. The fitting routine itself yields an error for b of ∆b = ±10 −4 . However, this is not the true error for the exponent in our analysis as the uncertainty in the initial time (i. e. the endpoint of the phase healing process) used in the fitting routine is much more severe. Thus, we vary the initial time by ∆t = ±30 to find ∆b = ±0.02 which appears to be a more realistic estimate for the error.
In figure 14 we show d(t) on a logarithmic scale for all initial separations. Together with the data we also plot as the black line the power law (D.1) with b = 0.53.

E Finite Difference Methods
To compute the velocities and accelerations from the trajectories of the vortex-anti-vortex system we use finite difference methods. Finite differences allow for a good approximation of the derivatives of our quasi-continuous vortex trajectories.
As we use an adaptive time-stepping scheme, our timesteps are unevenly spaced. This has to be taken into account when computing finite differences. Nevertheless the calculation is straightforward. For simplicity, we restrict ourselves here to one spatial dimension in the following, the generalization to two dimensions being obvious. We start out by Taylor expanding the position of one of the vortices at times t i+1 and t i−1 , x(t i+1 ) = x(t i + ∆t) = x(t i ) +ẋ(t i ) ∆t +ẍ(t i ) (∆t) 2 2 + . . . , (E.1) where a dot denotes the derivative with respect to the time t and we introduced both ∆t = t i+1 − t i and ∆t = t i − t i−1 to account for the in general not evenly spaced timesteps. Computing the difference (E.1) − (E.2) and neglecting all terms of second or higher order in the derivative gives the local velocity of the vortex at time t i , v(t i ) =ẋ(t i ) = 1 2 We are also interested in the vortex acceleration. For that, we compute the sum (E.1) + (E.2) and neglect all terms of third or higher order in the derivative. We obtain which, after some algebra, gives us the central second finite differences derivative used in this work,