A damped point-vortex model for polar-core spin vortices in a ferromagnetic spin-1 Bose-Einstein condensate

Ferromagnetic spin-1 Bose-Einstein condensates in the broken-axisymmetric phase support polar-core spin vortices (PCVs), which are intimately linked to the nonequilibrium dynamics of the system. For a purely transversely magnetized system, the Turner point-vortex model predicts that PCVs behave like massive charged particles interacting via a two-dimensional Coulomb potential. We test the accuracy of the Turner model for two oppositely charged PCVs, via comparisons with numerical simulations. While the bare Turner model shows discrepancies with our numerical results, we find that a simple rescaling of the PCV mass gives much better agreement. This can be explained via a phenomenological damping arising from coupling to modes extrinsic to the point-vortex phase space. We also identify the excitations produced following PCV annihilation, which help elucidate recent phase ordering results. We extend the Turner model to cases where the system is magnetized both transversally and axially, identifying a crossover to scalar vortex dynamics for increasing external Zeeman field.


I. INTRODUCTION
Spinor Bose-Einstein condensates (BECs) can exhibit both ferromagnetic and antiferromagnetic order, and possess a rich array of phases with distinct symmetry properties. Associated with these phases are a variety of topological defects and spin textures [1,2], which play an essential role in nonequilibrium processes such as symmetry breaking and the Kibble-Zurek mechanism [3,4], phase ordering dynamics [5] and quantum turbulence [6,7]. Most work on defects and spin textures in spinor condensates have focused on topological and stability aspects , with less exploration into defect-defect interactions and consequent dynamics. Studies on the interactions and resulting dynamics of half-quantum vortices in antiferromagnetic spin-1 condensates [34][35][36][37] and the collisional dynamics of non-abelian vortices in spin-2 condensates [38,39], as well as the dynamics of vortex dipoles across distinct magnetic phases [40], are notable exceptions, and reveal the rich dynamics possible due to the multicomponent nature of defects in spinor condensates.
The ground state manifold of a ferromagnetic spinor BEC is SO(3), supporting both nonsingular and singular defects [21]. Through the use of a quadratic Zeeman field, the spin vector can be constrained to point in a plane transverse to the field direction, and hence has SO (2) symmetry. This phase is termed the broken-axisymmetric (or easy-plane) phase. In this work we explore the dynamics of polar-core spin vortices (PCVs), which are point defects arising in the broken-axisymmetric phase. These vortices consist of a phase winding of transverse spin density around a polar (i.e. unmagnetised) core, and have been observed in situ in experiments [41]. PCVs have been shown to play a crucial role in the symmetry breaking following a quench to the easy-plane phase [41,42], leading to Kibble-Zurek scaling [43][44][45] and universal long-time phase ordering dynamics [46][47][48][49].
In scalar condensates, the point-vortex model [50] has been shown to describe accurately a plethora of nonequilibrium dynamical effects [51][52][53][54][55][56][57][58][59][60][61][62][63]. An analogous point-vortex model of PCV dynamics was introduced by Turner [64], where it was argued that the constituent circulations that make up a PCV interact like scalar vortices, but are confined due to spinexchange interactions. This confining energy manifests itself as a mass, resulting in PCVs behaving like massive charged particles interacting via a two-dimensional Coulomb force. We term this model the "Turner model". The presence of a vortex mass makes PCV dynamics vastly different from scalar vortices. With scalar vortices, two like charged vortices circulate around their centroid, while two oppositely charged vortices move in parallel lines (in the absence of damping) [65][66][67]. In contrast, the Turner model predicts that two like charged PCVs should repel, while two oppositely charged PCVs should attract. The qualitative features of the Turner model have been confirmed for two oppositely charged PCVs via simulations of the spin-1 Gross-Pitaevskii equations [68]. However this work also identified quantitative discrepancies between simulations and the Turner model, and hence left questions open regarding the precise quantitative details of the process.
In this work we carry out a detailed comparison between the Turner model of PCV dynamics and numerical simulations of two oppositely charged PCVs. While the bare Turner model shows deviations from numerical simulations, we find that the agreement can be vastly improved by introducing a phenomenological damping of the modes that give rise to the PCV mass, resulting in a simple rescaling of the bare PCV mass. The damping arises due to a coupling of the PCV coordinates to extrinsic modes, resulting in a loss of energy from the point-vortex phase space. We also identify the excitations produced following PCV annihilation, which likely play an important role in the anomalous phase ordering identified in [69]. A linear Zeeman field results in the ground state condensate magnetization rotating out of the transverse plane to partially align with the Zeeman field [70]. We extend the Turner model to PCVs in this phase, finding good agreement with our model and numerical simulations.
The paper is organized as follows. In Sec. II we present background material to understand PCVs and the Turner model. In Sec. III we compare the Turner model with numerical simulations for two oppositely charged PCVs, and identify a crucial damping processes missing from the original Turner model. In Sec. IV we study the excitations produced following PCV annihilation and in Sec. V we extend the Turner model to describe systems with both axial and transverse magnetization. In Sec. VI we conclude, with various ideas for future work.  The atoms in a spin-1 BEC have access to three spin levels m = −1, 0, 1 and can be described by a spinor of three classical fields Ψ(x) = (ψ 1 (x), ψ 0 (x), ψ −1 (x)) T . In a flat-bottomed quasi-2D trap, the Gross-Pitaevskii equations (GPEs) are [1],

II. BACKGROUND
with [71][72][73] The spin-1 atoms interact via spin-independent (strength g n ) and spin-dependent (strength g s ) interactions, with n(x) = Ψ † (x)Ψ(x) the areal density and F(x) = µ=x,y,z Ψ † (x)f µ Ψ(x)ŝ µ the areal spin density for spin-1 Pauli matrices f µ and spin directionsŝ µ . We consider ferromagnetic (g s < 0) spin interactions, arising for example in 87 Rb [74,75] and 7 Li [76], which favours a nonzero magnetization density. A magnetic field alongŝ z , along with microwave dressing techniques [77,78], results in a Zeeman splitting of the spin levels, and has both a linear p and quadratic q contribution. The linear Zeeman shift can be combined with the Lagrange multiplier λ that enforces conservation of total axial magnetization and hence p = p + λ [79]. For 0 < q < q 0 = 2|g s |n 0 (n 0 is the mean condensate density) and |p| < q the ground state is in the brokenaxisymmetry (BA) phase, whereby the quadratic Zeeman energy favours a transverse magnetization that breaks the axial symmetry of the Hamiltonian [70,79]. Here q 0 is a quantum critical point separating the p = 0 BA phase from from the polar [Ψ = (0, ψ 0 , 0) T ] phase.

B. Polar-core spin vortices
For p = 0 in the BA phase, the quadratic Zeeman field confines the spin to point entirely in the (ŝ x ,ŝ y ) plane and the ground state has only transverse magnetization. The ground state can be parameterized as [70] Ψ = e iθ e ifzϕ n 0 2 with U(1) symmetries due to global phase rotations e iθ and spin rotations e −ifzθ aroundŝ z . Here cos(2β) = q/q 0 . The resulting spin density is A 2πκ (κ ∈ Z \ {0}) phase winding of the transverse spin F ⊥ = F x + iF y results when ϕ = κφ(x), with φ(x) = phase(x + iy) the phase of the complex number x + iy. When there is no additional global phase rotation 1 , the state is a PCV [1], This PCV state is approximate as in reality a core will form, to avoid the divergences of ∇ψ ±1 . Hence the density and potentially the phase profiles will be modified for distances |x| ξ s , with ξ s = / √ M q 0 the spin healing length. Within the core there is a peak in the occupation of the ψ 0 (polar) component, hence the terminology "polar-core spin vortex". As we will see, the core structure plays an integral role in the PCV motion, drastically changing the dynamics compared to the scalar case. Since the ψ ±1 components of a single PCV circulate in opposite directions, the velocity field generated by a collection of PCVs consists of counter-flowing ψ 1 and ψ −1 currents. Hence the flow field around a PCV transports no mass but will transport axial magnetization F z = |ψ 1 | 2 − |ψ −1 | 2 [80].

C. Scalar point-vortex model
In a scalar BEC, the point-vortex model is derived by assuming that the superfluid can be described by a collection of vortex coordinates X k that give the points of δ-function divergence of the vorticity of the superfluid velocity field. One then stipulates a vortex ansatz for the condensate field, are the vortex charges), which evolves according to the scalar Gross-Pitaevskii equation. Deriving equations of motion for the vortex positions X k (t) using either variational Lagrangian [81] or hydrodynamic [82] techniques gives the point-vortex equations of motion, Equation (6) shows that a scalar vortex follows the velocity field generated by the remaining vortices. This is necessary to conserve the fluid momentum in the presence of no external forces, in analogy to the Magnus effect in classical fluid dynamics [83][84][85][86]. Two fundamental assumptions of the scalar point-vortex model are that interactions with sound waves can be ignored and that the precise core features of the vortices do not affect the mesoscopic dynamics of the system, hence the use of δ-function cores suffice. In contrast, we will show that the core features of a PCV have a radical effect on the mesoscopic dynamics of the system and that coupling to spin waves is non-negligible. Despite this, a point-vortex model is still extremely useful at describing the bulk dynamics of PCVs. transverse spin, but will otherwise have different properties to the PCVs explored here.

D. PCV point-vortex model
The variational Lagrangian method to derive the pointvortex model for scalar vortices, Eq. (6), can be adapted to PCVs [64,68]. We start with an ansatz consisting of a product of PCVs [68], with X mk (t) the centre of a vortex in spin component m, κ k the PCV charge and g mk (x, t) the amplitude profiles giving rise to the PCV core. The X 1k and X −1k coordinates that constitute a single PCV will tend to move in opposite directions, due to the counter-flowing ψ ±1 fields generated by other PCVs and Eq. (6). This "stretching" is restricted, ultimately by the spin exchange term in g s |F| 2 , which favours the phase profiles of ψ 1 and ψ −1 to overlap. In addition, other energy terms arising from overlapping vortex cores in the ψ ±1 components may restrict the stretching. Hence the additional energy term H int [Ψ, Ψ † ] must be considered in the PCV dynamics.
The ansatz (7), as it currently stands, allows for coupling between phase and amplitude excitations via coupling between X ±1k and excitations of the amplitude profiles g mk . To obtain a point-vortex description of PCV dynamics, we make the approximation that the core profiles g mk are stationary, depending only on t through the vortex coordinates X mk , depends only on the vortex coordinates X mk . Assuming also that the contribution of each PCV to H int depends only on the magnitude of the PCV stretch |X 1k − X −1k | allows one to further decompose H int as for individual PCV stretch energies ∝ u(|X 1k − X −1k |) (the prefactor has been chosen to simplify the equations of motion below). Adding this to the scalar vortex interaction from Eq. (6), assuming that the PCV stretching remains small compared to the distance between PCVs, and neglecting the kinetic energy arising from the vortex cores, gives [64,68] with ∇ u =x∂/∂u x +ŷ∂/∂u y the gradient derivative with respect to a vector u = (u x , u y ). We see that the fluid momentum in component m is now balanced by the force ∝ −∇ X mk u [64]. We introduce position R k and stretch r k coordinates for each PCV, and hence rewrite Eqs. (9) aŝ Hence Eqs. (11) both take the form of a Magnus force balanced by the gradient of a potential.
The transverse spin will be zero at x = R k to avoid a singularity in the phase winding of ϕ. Equation (11b) predicts that the ψ ±1 flow fields will cause the PCV to stretch. This results in a force on the ψ ±1 components, arising from the stretch energy (8), with the net effect that two PCVs of opposite(same) charge attract(repel), Eq. (11a). The equations of motion (11) are derived from a time translationally invariant and hence energy conserving Hamiltonian, with a conserved dimensionless PCV energy [68], with (The condensate energy is related to Eq. (12) by scaling by the prefactor in Eq. (8), which includes the factor of π often included in the scalar point-vortex energy.) The assumption that the profiles g mk depend only on time through the coordinates X mk removes the possibility of coupling to dynamical degrees of freedom outside of the {X 1k , X −1k } phase space (or equivalently the {R k , r k } phase space). For scalar vortices, this decoupling is a reasonable approximation in many cases [81]. As we will later show, it is in fact unreasonable for PCVs. However, a simple linear damping term in Eq. (11b) allows for the effects of additional modes to be included, leading to an accurate pointvortex model for the coordinates R k .

III. TESTING THE TURNER MODEL OF PCV DYNAMICS
We analyze the dynamics of two oppositely charged PCVs with centres of circulation at positions R 1 (0) = (−10ξ s , 0) (charge κ 1 = −1) and R 2 (0) = (10ξ s , 0) (charge κ 2 = 1), see Fig. 1(a). The initial two-PCV state is, which we evolve using the spin-1 GPEs (1). For times −10t s ≤ t < 0 (t s = /q 0 ) we include in Eq. (1) an energy damping term of strength γ 2 . This allows a core structure to develop, while having only a small affect on |R 1 − R 2 | as long as γ 1. For times t ≥ 0 we set γ = 0. An example case of PCV dynamics obtained from our simulations is shown in Fig. 1. As predicted by Eqs. (11), the oppositely charged PCVs attract. Within each PCV, the ψ ±1 vortex cores separate, along a line orthogonal to R 2 −R 1 , see Fig. 1 [In all simulations we use n 0 = 10 4 ξ s and g n = 10|g s |. We solve Eq. (1) on a N = 1024 × 1024 grid with physical size 100ξ s × 100ξ s using a fourth order Runge-Kutta method (time step = 0.002t s ) with periodic boundary conditions and kinetic energy operator evaluated to spectral accuracy. The PCVs remain far from the boundary, which mitigates the effect of image charges and an initial phase discontinuity along x = (x, ±100ξ s ). To detect vortices, we interpolate the wavefunction onto a 100×denser grid in a region around each vor-tex and find points of diverging vorticity ∇ × ∇ phase(ψ ±1 ).]

A. PCV position dynamics
Solving the Turner model (11) for the PCV coordinates R k and r k requires stipulating a form for the potential u(r k ). We find that a harmonic potential describes the dynamics very well, with a a fitted "spring constant". Taking a second time derivative of Eq. (11a) and using Eqs. (11b) and (15) gives with, The prefactors in Eq. (16) are such that the right-hand side is the gradient (with respect to R k ) of the condensate kinetic energy. Equation (16) predicts that PCVs behave like massive charged particles moving under the influence of the twodimensional Coulomb interaction 3 . In Fig. 2(a) we show simulation results for the PCV separation D(t) = |R 2 (t) − R 1 (t)| for different q values. We test Eq. (11a), assuming a potential (15), by comparing D(t) with the integral 2a t 0 dτ r(τ )/t s . Here r = |r 1 | = |r 2 | is the stretching of either PCV, which we obtain from the GPE simulations. Equation (11a) is satisfied accurately for q 0.1q 0 , with a as a fitting parameter. For small q, the full SO(3) spin manifold will become accessible and hence we expect deviation from the Turner model 4 .
Equation (16) can be solved analytically for the two-PCV setup in Fig. 1 to give, Here 3 PCVs of charge |κ k | > 1 decay into PCVs of lower charge, which repel due to the repulsion of like-charged PCVs predicted by Eq. (16). 4 The fit for 0 < q 0.1q 0 can be improved by adding a term ∝ r k in the potential (15).
is the analytic prediction for the PCV collision time t coll , with c s = ξ s /t s the characteristic spin-wave speed, and allows for a nonzero initial PCV velocityḊ(0), which is present in our simulation due to the initial damped evolution, see enlarged region in Fig. 2(a). We can easily extract the PCV collision time t coll from our GPE simulations. Rescaling time by t coll results in the curves for different q values from Fig. 2(a) reducing to the same functional form, see Fig. 2(b). Hence the entire q dependence of the dynamics of D is contained in the parameter t coll , as predicted by Eq. (18). The agreement with Eq. (18) is excellent with t ann coll replaced by t coll and b ≈ 0.17. We find that the predicted scaling of (t coll ) 2 ∝ 1/a from Eq. (19) (to zeroth order inḊ(0)) is also very well satisfied, see Fig. 2(c). In addition, we find that t coll ∝ D(0) (not shown), as predicted by Eq. (19). The initial PCV velocity, obtained from Eq. (20) √ a, and decreases as q is increased, see enlarged region in Fig. 2 This can be interpreted as an enhancement of the PCV mass to the effective mass The necessity of using an effective mass, rather than the bare PCV mass, is the first evidence of a coupling to extrinsic degrees of freedom, outside of the {R k , r k } phase space. We will provide more evidence of this shortly, where we will see that the scaling factor A can be written in terms of a damping rate of the coordinates r k . The notion of an effective mass arising from a coupling to extrinsic modes is reminiscent of the notion of electron effective mass in solids [87]. Note that we see no q dependence on the enhancement factor in Eq. (22), hence the q dependence is entirely contained within the bare mass m v .

B. PCV stretch dynamics
We now turn to the PCV stretch coordinate r. In Fig. 3(a) we show simulation results for r(t) for different q values. We test Eq. (11b) by comparing r(t) with the integral t 0 dτ 2/D(τ ) for the two extreme cases q = 0.1q 0 and q = 0.9q 0 (the remaining t 0 dτ 2/D(τ ) results fall between these two curves). Evidently, the stretch coordinate r is much smaller than the prediction of Eq. (11b). Note that the oscillation of r observed in [68] is not present in Fig. 3(a). This is due to the initial damped evolution, which removes the energy liberated when the core forms and the resulting long-lived oscillation. The analytic prediction for r, obtained from solving Eqs. (11) with the potential (15), is Results for √ ar(t) versus t/t coll are shown in Fig. 3(b), where we see that the different √ ar(t) show poor collapse onto a single functional form (in contrary to the prediction of Eq. (23)). Furthermore the analytic prediction (23) (with t ann coll → t coll ) overestimates the stretching. This suggests that energy from the Coulomb potential H R is not being entirely transferred into the energy term H r , in contrary to the conservation law Eq. (12).
In Fig. 3(c),(d) we plot the two energy terms from Eq. (12) using the simulation results for R k and r k , with the potential (15) for H r and a obtained from the fits in Fig. 2(a). The PCV energy H R decreases as the PCVs collide, Fig. 3(c) while the energy H r increases, Fig. 3(d). However, the loss in H R is approximately four times greater than the gain in H r and hence most of the liberated H R energy is lost to other excitations in the system, rather than coupling to the PCV dynamical variable r.

C. The damped Turner model
We now present a simple modification to Eq. (11b) that gives rise to the effective mass, Eq. (22), and which also improves the agreement in Fig. 3(b). Motivated by the loss of energy identified in Fig. 3(c),(d), we introduce a phenomenologically damping term −Γẑ ×ṙ k on the right-hand side of Eq. (11b). Solving Eq. (11a) and the modified Eq. (11b) gives the same solution (18) but with the modified PCV collision time Eq. (21). The scaling factor A is related to the damping rate Γ via, Meanwhile Eq. (23) is modified to Using the effective mass Eq. (22) gives Γ ≈ 6, which when used in Eq. (25) gives good agreement with the r curve of the smallest q value, see Fig. 3(b). Note Fig. 3(b) suggests a q dependence on the scaling (25) that is not accounted for by our simple damping model. A more accurate determination of this damping process is an interesting topic for future work. The effective mass, however, appears insensitive to this q dependence and hence our phenomenological damping theory describes very well the behaviour in Fig. 2. Note also that despite the addition of this damping and the new parameter Γ, the effect on the dynamics of D is to simply modify the PCV collision time.

D. Probing the PCV stretch energy
In [64] it was proposed that the stretch energy (13b) arises from the spin exchange interaction ∝ g s d 2 x Re(ψ 2 0 ψ * 1 ψ * −1 ), which favours phase coherence between ψ 1 and ψ * −1 . The spin exchange energy will increase as a PCV stretches and the circulating ψ ±1 phases separate. We expect this energy to dominate for sufficiently large stretching where the cores of the ψ ±1 vortices do not overlap and the |ψ m | 2 profiles are independent of the stretching. However, for smaller stretching, energy terms involving |ψ m | 2 may also be important. In Fig. 4(a) we plot the density profiles δ|ψ m | 2 = |ψ m | 2 − |ψ gs m | 2 of the three spin components ψ m through the centre of the rightward moving (negatively charged) PCV along y (the stretch axis). Here |ψ gs m | 2 are the PCV-free ground-state densities computed from Eq. (3). Clearly the stretching remains small compared to the PCV core size and the |ψ m | 2 profiles change as the PCV stretches. Hence there will be a coupling between the stretch coordinate and energy terms involving |ψ m | 2 . Note that as the cores in the ψ ±1 components separate, each core is partially filled by the other ψ ∓1 component, resulting in the dipole of F z magnetization in Fig. 1(b). The density δn = n − n 0 is also plotted in Fig. 4(a).
In Fig. 4(b) we plot different contributions to the condensate energy, computed across a 12ξ s × 12ξ s area S centred The total mass density δn = n − n0 is also shown (red curve). (b) Within a 12ξs ×12ξs area around either PCV, both the spin exchange energy Ese (blue dots) and the quadratic Zeeman energy Eq (yellow triangles) increase as the PCV stretches. The increase in Ese is compensated by a decrease in the remaining terms in EF ⊥ , such that EF ⊥ (green stars) is close to constant. The energy EF z (red squares) decreases as the PCV stretches and the Fz dipole forms. The energy En (orange asterisks) remains very close to constant and hence excitations of the total mass are decoupled from the PCV dynamics. The zero point of each energy has been shifted so that the smallest r value plotted is at zero energy. For (a),(b) q = 0.3q0. (c) The spring constant a from Fig. 2(a) increases as a ≈ q/q0. at either PCV, for q = 0.3q 0 . The spin interaction energy consists of terms E Fz = (g s /2) S d 2 x F z (x) 2 and E F ⊥ = (g s /2) S d 2 x |F ⊥ | 2 . The latter term contains the spin exchange energy E se = 2g s S d 2 x Re(ψ 2 0 ψ * 1 ψ * −1 ). Changes in |ψ m | 2 may also change the quadratic Zeeman energy E q = q d 2 x (|ψ 1 | 2 + |ψ −1 | 2 ). Both the quadratic Zeeman energy E q = q d 2 x (|ψ 1 | 2 + |ψ −1 | 2 ) and the spin exchange energy E se = 2g s S d 2 x Re(ψ 2 0 ψ * 1 ψ * −1 ) increase, roughly in proportion, as the PCV stretches. Therefore it is not clear whether the stretch energy arises from E q or E se . The energy E Fz = 2g s S d 2 x F z (x) 2 becomes increasingly more negative as the PCV stretches, as expected due to the formation of the F z dipole. The energy E F ⊥ remains approximately constant and hence the increase in E se must coincide with a decrease of the remaining terms in E F ⊥ . We also plot the density interaction term E n = (g n /2) S d 2 x n(x) 2 , which remains very close to constant. Hence although the mass density variation across the PCV core is not negligible, see Fig. 4(a), the mass density modes do not couple strongly to the PCV dynamics, as expected since g n |g s |. Although our analysis of Fig. 4(b) provides some clues into the origin of the stretch energy, in general it is difficult to identify what excitations are contributing to H r and what excitations are contributing to the the loss of the PCV energy H [ Fig. 3(b),(c)]. Furthermore, unlike in Fig. 4(b), for q 0.5q 0 we find no clear trend in E q and E se with increasing r. This may relate to the qualitative change in the instability of a uniform polar state at q = 0.5q 0 [43][44][45]. A more detailed analysis of the PCV core excitation spectrum, as has been done recently for nematic vortices in an antiferromagnetic spin-1 condensate [32], may help elucidate more clearly what modes contribute to H r and what modes contribute to the loss of H. The stretch energy spring constant a, obtained from the fits to the results in Fig. 2(a), increases close to linearly with increasing q, see Fig. 4(c), with a ≈ q/q 0 5 . It is difficult, however, to associate this behaviour with E q or E se , as the density and phase profiles within the PCV core are not uniform. A detailed study of a stretched PCV excitation spectrum may reveal modes that increase in energy as the PCV stretches and that scale with q as in Fig. 4(c), which would provide more insight into the basis of the PCV stretch energy.

IV. PCV ANNIHILATION
Two oppositely charged PCVs attract until they collide and annihilate. In this section we examine the dynamics of the excitation produced following PCV annihilation. Figure 5 shows the evolving (a) transverse spin and (b) F z spin density, following the PCV annihilation imminent in Fig. 1(a). The ringshaped excitation that propagates out from the PCV annihilation point consists of variation in both the transverse spin direction and the F z spin density.
In Fig 5(c) we plot F z (x) along the line x = (0, y), extracted from Fig. 5(b). The excitation is initially dominated by a single peak, but this disperses as the excitation propagates out. We estimate the speed of the propagating waves from the propagation speed of the minimum of F z (0, y). The position y ex of this minimum travels at a constant speed v ex , see Fig. 5(d), which increases for increasing q. The speed is closely matched by the speed of long-wavelength Bogoliubov excitations of F z , which dynamically couple to quadratic excitations of phase(F ⊥ ), see inset to Fig. 5(d). For g n |g s | these excitations both have the dispersion relation ω k = (q 0 ξ s k/2) ξ 2 s k 2 + 2q/q 0 [73], which travel at a speed v Bog = ω k /k. Fitting v Bog to v ex , we extract a small but non-zero wavevector k ≈ 2π × 0.1ξ −1 s , which shifts the speed slightly from the sound (k = 0) speed. In [69], which studied the phase ordering dynamics following a quench to the easy-plane phase, it was found that out-of-equilibrium F z and phase(F ⊥ ) excitations remain after all PCVs have annihilated, leading to anomalously slow thermalization. However the origin of the excitations was not identified. Following the results from Fig. 5, we expect that the excitations are produced from the PCV annihilation events during the PCVdriven phase ordering.

(27)
In addition to the transverse magnetization F ⊥ , the state (26) has a non-zero axial magnetization, A phase winding of the transverse spin gives rise to a PCV, with the state (7) modified to Carrying out an identical procedure that leads to the derivation of Eq. (9) (see [68]), we obtain According to Eqs. (30), the Coulomb potential between two oppositely charged PCVs initially separated along the x-axis causes the ψ ±1 components to stretch along the y-axis, like in the F z = 0 case. With non-zero F z magnetization, however, the prefactor of the balancing stretch force ∝ ∇ X mk u(|r k |) is different for each of the ψ ±1 components. As a result, the centres of circulation of the ψ ±1 components move along the x-axis at different velocities. This results in an additional stretching along the x-axis, which needs to be balanced by a Magnus force along the y-axis. Hence in addition to the motion along the x-axis, the PCVs also move along the y-axis.
Continuing as for Eqs. (11), we rewrite Eqs. (30) in terms of the coordinates R k and r k . Including a damping term discussed in Sec. III C giveŝ B. PCV dynamics with small and large axial magnetization Figure 6 shows the PCV dynamics from GPE simulations as in Sec. III, but with the initial state (14) modified according to Eq. (29) such that F z = 0.1n 0 . Here F z = A −1 d 2 x F z is the spatially averaged F z density (A is the system area). As predicted, the PCVs not only attract but also move transverse to their separation, see Fig. 6(a). The dynamics of R k follows closely the integral [( κ k )/(2M sin 2 2α)] t 0 dτẑ × ∇ r k u(|r k (τ )|) using the potential (15), see Fig. 6(b). The spring constant a is fitted using GPE values for |R k | and |r k | (noteẑ × ∇ r k r 2 k = 2ẑ × r k ). The dynamics of r 2 is shown in Fig. 6(c), compared with the prediction of Eq. (31b) using simulation results for R 2 . Equation (31b) describes the dynamics of r 2 reasonably well with damping rate Γ ≈ 3. Figure 7 shows the PCV dynamics from GPE simulations with F z = 0.2n 0 . For this higher magnetization, we find that the stretch coordinate r k rotates until it is parallel R k , see Fig. 7(a). Hence the attractive force between the PCVs becomes negligible and instead the two PCVs move parallel, analogous to scalar vortices. (This effect becomes more pronounced as | F z | is increased.) The dynamics of R k follows closely the integral [( κ k )/(2M sin 2 2α)] t 0 dτẑ × ∇ r k u(|r k (τ )|) using the potential (15), see Fig. 7(b). We find that a/ sin 2 (2β) has little dependence on F z for | F z |/n 0 0.2, with a/ sin 2 (2β) = 0.46, 0.45, 0.48 for F z /n 0 = 0, 0.1, 0.2 respectively. As observed in Fig. 7(a), r 2 rotates until y 2 ≈ 0, see Fig. 7(c). The small stretch x 2 also approximates a steady state. For this larger magnetization, the prediction of Eq. (31b), using simulation results for R 2 , is less accurate, but still provides a qualitative description of the dynamics with damping rate Γ ≈ 3. The precise nature of the oscillations in Fig. 7(c) are sensitive to the length of the initial damped evolution, but the overall trend of the curves remains unchanged.
The crossover to scalar vortex dynamics is interesting and warrants further discussion. Settingṙ k = 0 in Eq. (31b) gives, identical in form to the point-vortex equations for scalar vortices, Eq. (6). For the two-PCV setup in Fig. 6, 7, Eq. (32) allows for constantṘ k and hence a constant r k is compatible with Eq. (31a). We expect that this crossover from PCV dynamics to scalar vortex dynamics will have interesting consequences in PCV-driven phase ordering, following quenches to the BA phase. For zero F z , the phase ordering is governed by a dynamic critical exponent z = 1 [46,47]. For sufficiently large | F z |, we expect the attractive force between PCVs to become negligible compared to scalar-vortex damping [61], and the PCVs will behave like scalar vortices. In this scalarvortex regime we expect critical behaviour that follows that of a scalar condensate, with z = 2 [88][89][90].
Considering an arbitrary | F z | > 0, we envisage two possible phase-ordering scenarios, dependent on whether or not the scalar-vortex regime is reached prior to PCV annihilation. Scenario 1: The scalar vortex regime is always reached prior to PCV annihilation, given a sufficiently large PCV separation, and hence the long-time limit of phase ordering with nonzero F z would be that of a scalar system. This is analogous to the small-q phase ordering around the isotropic ferromagnetic phase [91]. Scenario 2: There exists some critical value of | F z | below which the scalar-vortex regime is not reached prior to PCV annihilation, for any initial PCV separation, and hence the BA phase would be divided into two dynamic universality classes either side of this critical value of | F z |. A many-PCV simulation that allows for soundinduced damping [61] would be needed to definitively determine which scenario occurs, however identifying regimes where Eq. (32) is reached from Eq. (31) may give insight into this question. For example, for the two-PCV setup in Figs. 6, 7, one could compare t coll with the time for r k to rotate parallel to R k , as a function of F z and initial PCV separation.

VI. SUMMARY AND OUTLOOK
In this work we have systematically studied the accuracy of the Turner point-vortex model at describing the dynamics of two oppositely charged PCVs, via comparisons with GPE simulations. While the bare Turner model, as introduced in [64], shows discrepancies with GPE results, the agreement can be vastly improved by introducing a phenomenological damping of the PCV stretch coordinate. The damping arises from coupling to other excitations not accounted for in the point-vortex phase space, and reduces the PCV stretching. The net effect on the mesoscopic dynamics, however, is to simply increase the PCVs bare mass to give an effective mass that increases the PCV collision time. In addition, we identified the excitations produced following PCV annihilation and extended the Turner model to describe PCV dynamics with axial magnetization, both of which have applications to the phase ordering of spin-1 BECs. Some important questions raised by this work remain to be answered. Firstly, can the modes that contribute to the PCV stretch energy be identified; and secondly, can the modes that are responsible for the damping of the stretch coordinate be identified. In addition, a derivation of the Turner point-vortex model from the hydrodynamic formulation of a ferromagnetic spin-1 condensate [80] may offer additional insights into the fluid dynamics of the Turner model, analogous to the derivation of the scalar point-vortex model from the Euler equation [92]. For example, the scalar point-vortex model, which describes mass currents, reflects conservation of momentum. As the circulating transverse spin around a PCV corresponds to an F z magnetization current, Eqs. (11) and their modifications Eqs. (31) may arise from an analogous conservation law for F z "spin momentum".
A useful property of PCVs is that they can be observed in situ in experiments [41] using non-destructive phase contrast imaging [93]. To observe the attraction of two PCVs, the axial magnetization should be kept small, | F z | 0.2n 0 . The effective mass of two PCVs could be obtained from measuring the PCV collision time. The spring constant a could be measured with highly resolved absorption imaging [94,95] of the F z dipole, which forms as the PCV stretches. Combined with the effective mass, this would allow extraction of the damping rate Γ.
A host of possible applications of PCVs could also be explored, for example few-body PCV dynamics, PCV driven turbulence, emergent hydrodynamics [63,96,97], and thermodynamic properties of PCV systems. Also, the coupling to other dynamical modes through the PCV stretching and the associated modification of the PCV mass raises the question of whether PCV dynamics could be manipulated via the use of additional interactions or external fields. The success of our damped Turner model provides a useful starting point for these explorations.