Exactly solvable model for a solitonic vortex in a compressible superfluid

Vortex motion is a complex problem due to the interplay between the short-range physics at the vortex core level and the long-range hydrodynamical effects. Here we show that the hydrodynamic equations of vortex motion in a compressible superfluid can be solved exactly in a model"slab"geometry. Starting from an exact solution for an incompressible fluid, the hydrodynamic equations are solved with a series expansion in a small tunable parameter provided by the ratio of the healing length, characterizing the vortex cores, to the slab width. The key dynamical properties of the vortex, the inertial and physical masses, are well defined and renormalizable. They are calculated at leading order beyond the logarithmic accuracy that has limited previous approaches. Our results provide a solid framework for further detailed study of the vortex mass and vortex forces in strongly-correlated and exotic superfluids. The proposed geometry can be realised in quantum-gas experiments where high-precision measurements of vortex mass parameters are feasible.


I. INTRODUCTION
The quantisation of vortices in superfluids is one of the most prominent quantum phenomena of condensed matter physics. Yet the dynamical properties of vortices are poorly understood, with the mass [1][2][3] and applicable forces [4][5][6] still the subject of controversies after decades of research [7][8][9][10][11][12]. Recent progress in ultra-cold atomic gases has demonstrated that the dynamics of individual vortices can be observed in defect-free superfluids near zero temperature [13][14][15][16][17]. In particular the determination of the inertial to physical mass ratio from the observed oscillation frequency of a solitonic vortex in a trapped unitary Fermi gas in Ref. [15] has demonstrated that, in principle, high precision measurements of the dynamical properties are feasible [18]. The ability to derive a deeper understanding of vortex motion from such experiments is currently impeded by the limited precision of available theoretical predictions. For harmonically trapped quantum gases [16,[19][20][21] these are typically performed at logarithmic accuracy, i.e. up to an undetermined factor that varies only slowly (logarithmically) with system parameters. In addition to the longrange hydrodynamic laws governing the vortex motion, one can also expect the short-range, mesoscopic structure of the vortex core and its interaction with the fluid to influence the precise value of the vortex mass and other dynamical properties. This is particularly interesting for strongly-correlated quantum fluids such as superfluid helium [22] where no quantitative microscopic theory exists. In Fermi superfluids such as neutron matter [23] and the unitary Fermi gas [24] one may further expect the Kopnin [6] and Iordanskii forces [4] to become relevant even though their effects may easily be dwarfed by dominant contributions from long-range hydrodynamics including the Magnus force [25]. Hence it is important to have precise control and understanding of the long-range hydrodynamic effects on the vortex motion. In this work we propose a simple slab geometry that provides this precise control and solve the hydrodynamic Euler equations for the vortex motion by a series expansion, where each term can be computed exactly. Thouless and Anglin concluded from a theoretical study that the vortex mass is poorly defined as its value depends on the process of measurement [7]. While this remains true in general, we demonstrate that the vortex mass can in fact be properly defined if an appropriate context is provided. Here this is achieved by specifying the geometry of a long but transversely confined channel, or slab. A vortex line perpendicular to the channel's long axis is known as a solitonic vortex [26][27][28]. The concept was introduced in [26] in order to describe localised nonlinear wave excitations with vortex structure that travel at sub-sonic velocities along the channel [28]. Solitonic vortices have recently been observed as long-lived products from the decay of dark solitons [16,29,30] and as spontaneously formed defects during the Bose-Einstein condensation transition [31,32]. The solitonic vortex appears as a decay product of other families of non-linear localised excitations and is, in fact, the only known stable and slow-moving solitary wave excitation in channels that are sufficiently wide to host a vortex [33].
By virtue of the transverse confinement, the vortex becomes a localised object. Indeed, the exact solutions obtained in this work and visualised in Fig. 1 show that the associated velocity field and the corresponding excitation energy density are exponentially localised in the longitudinal direction on the length scale of the transverse dimension. Thus the vortex loses its long-range nature beyond the transverse length scale of the channel. The localisation makes it possible to consider the solitonic vortex as a quasiparticle with a well-defined inertial mass M * through the framework of Landau quasi- 1st order phase S1/(q 3 ξ 2 /D 2 ) (d) 1st order phase Numerical evaluation of the renormalised first-order compressible phase correction S1 with h/D = 0.7. (d) Same as b, but with the second-order density correction n2, which depends on S0 and S1. For n1 + n2, we have chosen q = 1, ξ/D = 0.05, and γ = 3/2, representing the unitary Fermi gas (UFG). Note that n1 is independent of γ, the polytropic index.
particle dynamics of solitary waves [34]. It allows us to derive a Newtonian equation of motion under the condition that the localised quasiparticle moves sufficiently slowly between regions of different background density, to suppress radiation and to conserve the quasiparticle's energy. The effects of a trapping potential are hereby encapsulated in the chemical potential of local equilibrium µ(X) = µ 0 − V trap (X), which is assumed to vary by negligible amounts over the length scale of the solitonic vortex. If E s (µ, V ) denotes the excitation energy associated with a solitonic vortex at chemical potential µ moving with velocity V =Ẋ, then requiring conservation of energy, dE s (µ(X),Ẋ)/dt = 0, leads to where is the inertial or effective mass. Since the vortex depletes particles from its core, we also have F = M ph g, a buoyancy-type force with the physical mass where g = m −1 dµ/dX is a buoyancy parameter with dimension of acceleration, and m the mass of the elementary bosons. While g characterises the environment of the solitonic vortex, the mass parameters M * and M ph are measurable and well-defined intrinsic properties that depend both on the microscopic details of the many-body physics, e.g. the vortex core structure, as well as the transverse confinement. In this way the Landau quasiparticle picture of the solitonic vortex provides a framework for defining vortex properties, where additional effects like non-conservative and fluctuating forces can be added later.
Here, we will calculate the mass parameters M * and M ph that arise as a consequence of the superfluid hydrodynamics given by the Euler and continuity equations. These equations govern the fluid dynamics at length scales large compared to the healing length ξ and other microscopic length scales. The ratio ξ/D of the healing length to the channel width D thus naturally emerges as a small parameter of the theory. The key to finding exact solutions is to choose a geometry where simple analytic solutions are available for the simplified model of an incompressible fluid. From this starting point, compressible corrections are later calculated as a perturbation series.
The paper is organised as follows: The main hydrodynamic theory is laid out in Sec. II, starting with the definition of the model and the analytic solution for the solitonic vortex in an incompressible fluid followed by the perturbation theory for the effects of compressibility in the Euler equations. Section III then reports on a number of results for the properties of the solitonic vortex at low orders of perturbation theory. The consequences of these results for the vortex dynamics under the influence of weak harmonic trapping in the quasiparticle dynamics framework are reported in Sec. IV and the relevance of the results for further experimental study are discussed in Sec. V. Two appendices discuss the derivation of the Euler equation in the form used in this work and an analytic approximation for the first order correction to the superfluid phase, respectively.

A. Solitonic vortex in an incompressible fluid
We consider a three-dimensional uniform slab geometry as sketched in Fig. 1a. Hard-wall boundaries confine the superfluid to a rectangular region in the yz-plane where the potential is uniform. Such a geometry can be realised experimentally with a flat-bottom trap for atomic gases [35]. For simplicity, we assume an infinite extent along the longitudinal x-direction [36], a width D in the y-direction and additional confinement < D along the third direction, with which the vortex is aligned. Since the third dimension is irrelevant for the further discussion we will work in the projected two-dimensional xy-plane. For the geometry of a two-dimensional channel the velocity field v 0 (x, y) of a vortex in an incompressible inviscid fluid can be determined with the method of images [25].
The velocity field is related by v 0 = m ∇S 0 to the superfluid phase S 0 (x, y) shown in Fig 1b. The phase S 0 of the incompressible superfluid with a vortex at position ρ and hard wall boundaries can be constructed by superimposing the known solution for a vortex without boundaries with those of image vortices, such as to satisfy the condition of no flow normal to the boundary in the resulting velocity field. While a single straight wall requires a single image vortex, the two parallel walls of the channel generate a (doubly) infinite array of image vortices. The problem can be solved in an elegant way on the complex plane by introducing the meromorphic velocity potential w(z) with poles as the vortex locations, where S 0 (x, y) = {w(x + iy)}. The solution for the channel can then be found with the help of a conformal transformation between the channel and the half plane, where a single image vortex suffices [25,37]: Here, q ∈ Z is the quantized charge with q = 1 for a right-handed singly-charged vortex, the location of the vortex is ρ = (0, h), and the two-dimensional coordinate space is now represented by the complex plane z = x+iy. The x-axis is taken parallel to the walls of the channel, the origin is at the bottom wall, and the y-axis passes through the vortex. S 0 is given by the real part of w and can be written as For practical purposes, it is useful to introduce the stream function χ(x, y) = {w(x + iy)}, which has the important property of being single-valued, as opposed to the phase S 0 . By virtue of the Cauchy-Riemann equations ∂ x S 0 = ∂ y χ and ∂ y S 0 = −∂ x χ any expression involving derivatives of S 0 are easily rewritten in terms of χ. Even though the phase singularity at the position of the vortex results in a locally divergent velocity field, meaning that for example the energy E s (Fig. 2a) obtained as an integral over the slab of the kinetic energy density is formally divergent, it turns out that our model can be renormalised. Importantly, the compressible corrections to the density and phase fields that we calculate here are renormalisable. This means that we can meaningfully assign several interesting properties of the solitonic vortex; in particular, we find that the vortex moves along the channel in the x-direction with a velocity V 0 = πq 2Dm cot(hπ/D) (Fig. 2b), where q is the vortex charge and h the distance from the left channel wall as in Fig. 1a. With the non-divergent canonical momentum P 0 = 2π qn 0 h (Fig. 2c), the inertial mass can be calculated from M * 0 = dP 0 /dV 0 as where n 0 is the two-dimensional density of the incompressible fluid. When the vortex is located at the centre of the channel (h = D/2), it is at rest (V 0 = 0), and the inertial mass M * 0 = −4mn 0 D 2 /π is simply proportional to the mass of the fluid contained in the square of the channel width D. While the inertia is thus determined by the fluid inside the localisation volume, the negative sign indicates that the solitonic vortex accelerates in the direction opposite to an applied force according to Eq. (1). The incompressible fluid model is directly relevant to real superfluids where it corresponds to the Thomas-Fermi approximation (density follows the local chemical potential) in a flat-bottom trap. Significantly, Eq. (6) for the inertial mass presents an exact result for strongly and weakly-correlated superfluids in the Thomas-Fermi limit ξ/D → 0, in contrast to previous results obtained for harmonically trapped superfluids [16,[19][20][21]. A (finite) healing length ξ, then, defined by dµ/dn = 2 /(2nmξ 2 ), is associated with the compressibility of the fluid. All physical fluids are always compressible.

B. Perturbation expansion for compressible corrections
Limitations of the incompressible model appear when we try to evaluate the force term of the equation of motion (1), wherein the physical mass formally vanishes. A suitable model that overcomes this limitation is that of inviscid isentropic flow described by the Euler equation, which in the co-moving reference frame of the vortex be- Here µ 0 is the bulk chemical potential, and µ[n] is the chemical potential in a local density approximation evaluated at n(r). For definiteness, we assume a polytropic equation of state n ∝ µ γ , where γ is the polytropic index; for example, for a weakly-interacting Bose-Einstein condensate γ = 1, while for the unitary Fermi gas γ = 3/2 [38]. The dimensionless factor is used for collecting orders in a perturbation expansion, and will be set to 1 at the end. We note at this point that the Gross-Pitaevskii equation for Bose-Einstein condensates (BECs) can be rewrit-ten as an Euler equation with an additional term known as quantum pressure [39] where we have defined the quantum pressure term (first term on the left) to be of order . Together with the usual continuity equation (see below), Eq. (8) is fully equivalent to Gross-Pitaevskii theory and provides a quantitative description of dilute-gas BECs. The quantum pressure term governs the detailed structure of solitons and vortex cores. Whether a similar term exists for the crossover superfluid Fermi gas, and the unitary Fermi gas in particular, is not known.
Complemented with the continuity equation the Euler equation (7) describes a compressible superfluid at zero temperature. Setting = 0 recovers the incompressible problem already solved with the method of images. The more general problem can be solved with a perturbation expansion for the density n = n 0 + n 1 + 2 n 2 +· · · and phase S = S 0 + S 1 + 2 S 2 +· · · . Requiring that n and S solve Eqs. (7) and (9), sorting orders of reveals that the density n k at order k can be expressed through derivatives of the lower-order phase fields, while the phase S k at order k can be determined from a Poisson equation with a source term containing derivatives of n k and lower-order phase fields. Through the perturbation expansion, we have essentially linearised the problem by writing the non-linear Euler and continuity equations in a different form as an infinite set of coupled, but linear Poisson equations. By integrating the Poisson equations it is thus possible to determine the phase and density fields to arbitrary order. Specifically, the density can be obtained using the Euler equation (7) and inverting the equation of state µ[n]. For definiteness we assume a power-law relation n 3D = αµ γ for the three-dimensional density [40]. Further assuming uniform confinement perpendicular to the xy plane with a length scale B < D, we obtain for the two-dimensional density The power-law equation of state covers in particular the case of a weakly-interacting BEC, where γ = 1 and µ = gn 3D , where g is the Gross-Pitaevskii coupling so that α = 1/g. For the unitary Fermi gas we have γ = 3/2, and where ξ mb ≈ 0.37 is known as the Bertsch parameter [41]. An important length scale related to the equation of state is the healing length ξ with µ 0 = γ 2 /2mξ 2 .
We expand n = n 0 + n 1 + 2 n 2 + . . ., S = S 0 + S 1 + 2 S 2 + . . ., where n k and S k are kth order terms in in the perturbation expansion. Comparing coefficients of in Eq. (10) gives the density corrections Expressions at arbitrary order can be generated by writing n as a generalised Cauchy product using Mertens' Theorem: . . .
Similarly, the phase corrections can be deduced using the continuity equation under the condition of stationary flow: ∇ · (nv) = − ∂n ∂t = 0. We get Clearly ∇n 0 = 0, ∇ 2 S 0 = 0, and so Equations (11) - (14) were derived from the classical Euler equation (7). If we include the quantum pressure term for BECs from Eq. (8), the results are modified by an additive correction to the density expressions of Eq. (11) at each order in with n i → n i + n (qp) i but otherwise unchanged (Note that γ = 1 for BECs). The quantum pressure corrections can be calculated order by order and read In Fig. 1, we show the results for the phase up to first order and the density up to second order. After renormalisation, every order of the phase and all compressible corrections to the density formally diverge only at the vortex position (Appendix B). We note that only S 0 has a logarithmic branch point. Even though our counting device is not a small parameter itself, it turns out that successive orders in also accumulate factors of ξ 2 , which is small compared to the relevant length scale D 2 . The addition of a quantum pressure term to the Euler equation (7), which would formally make it equivalent to the Gross-Pitaevskii equation (with γ = 1), will affect the density and phase correction terms at order O ξ 4 /D 4 .

III. SOLITONIC VORTEX IN A COMPRESSIBLE SUPERFLUID
In this section, we obtain several physical properties of the vortex from the compressible hydrodynamics developed in the previous section.

A. Velocity and phase step
The vortex speed is given by the superfluid flow field at the vortex location after subtracting the free vortex q/r divergence (and any divergences of the velocity corrections), wherer = |r − ρ| is the distance from the vortex. The vortex velocity is parallel to the channel walls: where θ is an angle measured at the vortex and a series expansion of V = V 0 + V 1 + · · · is defined through the corresponding expansion of the phase S. From the incompressible phase S 0 of Eq. (5) we easily obtain It turns out that ∇S 1 diverges at the vortex. Expanding the analytical approximation around the vortex, where we have taken only the leading term for S approx 1 .
We can see that once the divergingr −2 term is removed, taking the limitr → 0 gives the following first-order correction to the vortex velocity (in the x-direction): We show how to obtain the full V 1 (Fig. 2b) numerically below in Sec. III C. The phase step across the solitonic vortex can be obtained by evaluating the difference of the limits (x → ∞) − (x → −∞) of the phase field S 0 + S 1 . Using the analytical approximation: A direct numerical evaluation of Eq. (B3) yields an additional term of order ξ 2 /D 2 with the same sign.

B. Canonical momentum
Similar to the case of a dark soliton [42] or a vortex ring in a cylindrical wave guide [21], we need to distinguish between the physical momentum P ph and the canonical momentum P = P ph + n 0 D(2π − ∆S). The physical momentum is the integral of the momentum density of the solitonic vortex solution with open boundary conditions for x → ±∞ and a phase step ∆S. The canonical momentum is relevant in the context of Hamiltonian dynamics and is obtained by adding the momentum of a back flow current to compensate for the phase step. Indeed the canonical momentum can be obtained directly by integrating the x component of the velocity field where a background phase gradient has been added such that the solution fulfills periodic boundary condition (pbc) in x direction with box length L D: where S pbc = S + (2π − ∆S)x/L. The first term of the right hand side of Eq. (21) is easy to evaluate since ∂ x S pbc 0 dx has the result 2π for y < h and 0 for y > h. The result is (at all orders of when expanding S) There are no branch cuts or logarithmic branch points in S 1 . In this sense, the vortex is already contained within S 0 . Hence, the integrations become straightforward as we do not have to worry about crossing any branch cuts. At first order in we have where the second term on the right hand side vanishes in the limit L → ∞ and thus can be discarded. In the region Ω \ B R ρ , we can write the integrand of the first term I 1 as a total derivative: The integral is solved by means of the divergence theorem and gives Finally taking R = ξ we obtain for the momentum

C. Effective mass
The effective mass M * of the vortex can be calculated as where the canonical momentum P is given in Eq. (26). We write the velocity derivative in terms of an expansion β(h) ≡ dV /dh = β 0 + β 1 + . . ., where β i corresponds to a term of order (ξ 2 /D 2 ) i . The first term β 0 is easily found from Eq.
For the effective mass of a stationary solitonic vortex at h = D/2, this yields where |Ξ 1 | 10 −4 , as determined by our numerical accuracy. We note that the first-order correction to M * | h= D 2 has the same sign as the leading-order term provided that ξ D < 2 π e − 2+π 2 2π 2 ≈ 0.349, making the effective mass even more negative.
By definition, β(h) = dV /dh. For β 1 , we need V 1 = ( /m)∇S 1 at the vortex with the divergent 'free vortex' terms removed as in Appendix B and Sec. III A. Using Eq. (B4), we start by writing As before for S 1 , we renormalise β 1 by expanding the integrand about the vortex position as a power series in the distancer from the vortex. The diverging terms are analytic expressions that scale with inverse powers of the distancer, fromr −4 tor −1 . We then renormalise the original integrand by adding four counter terms. We integrate the resulting finite integrand numerically to obtain V 1 , which we then differentiate with respect to h to get β 1 . The numerically thus evaluated V 1 (whose only component is along x) is shown in Fig. 2b as the red solid line.

D. Number of missing particles and physical mass
Vortex lines are usually associated with a core region of depleted particle number density due to large kinetic energy densities near the vortex filament. This depletion is quantified by the number of missing particles which usually takes negative values. The missing particle number N s is closely related to the physical mass defined by Eq. (3). In fact, for a zero-velocity solitonic vortex [43][44][45] M ph = mN s , where the zero-velocity condition is equivalent to setting h = D/2. In the following we outline the procedure for calculating N s . Since the individual terms n i and any finite approximation to the series expansion of the density diverge near the vortex position, the integral (31) has to be renormalized. As a consequence, the approximate density may become negative in a small region near the vortex position. A renormalisation procedure can then be defined by excluding the area B from the integral where the approximate density has negative values, defined by n(r) = 0 at r ∈ ∂B. With n = n 0 + n 1 and ξ/D → 0 this condition yields a a disk B R ρ with a radius R = ξ. Assuming convergence of the infinite series expansion for the density n, the excision radius R will shrink to zero, as the density in the hydrodynamic description of a compressible superfluid vanishes only at the position of the vortex.
We start by formally expanding N s in powers of : From Eq. (11), considering only n 1 , we get Considering n 2 , we can write but this integral is more difficult to evaluate in closed form, and essentially would give S 1 as well.
The procedure of excising a disk of radius R, which we then identify with the healing length ξ, introduces complications with terms of order 2 and higher. In par-ticular, at h = D/2, we find both analytically and numerically the following dependence on R: The leading order terms do not depend on h. After the identification R = ξ, this means that we need all the terms in the perturbation series just to get the leadingorder O(ξ 2 ) term for N s . Fortunately, it is possible to evaluate this contribution from all the terms of order 2 and higher using the general expression (12). Summing over all the higher-order contributions of integrals of the form Ω\B R ρ |∇S 0 | n d 2 r, we obtain where p F q is the generalised hypergeometric function.
The quantum pressure corrections contribute to n 2 , and so have an effect on the missing particle number at order 2 and higher. As is the case without the quantum pressure, the excision procedure introduces problems when evaluating N s . For example, we can analytically evaluate the first quantum pressure correction which, as before, means that we must evaluate higherorder corrections to obtain the full contribution to N ( ) s . While this is possible in principle, we have not performed a fully detailed calculation to obtain this contribution.

IV. EQUATION OF MOTION UNDER HARMONIC TRAPPING
In the previous Secs. II and III we have considered a homogeneous slab extending infinitely in the x-direction, which led to the vortex moving at constant velocity parallel to the x-axis. In this section we consider the effects of a weak harmonic trapping potential V trap (x) = 1 2 m p ω 2 x x 2 , where m p is the mass of the elementary particles confined by the trapping potential [46]. An important property of the solitonic vortex is that it is localised along the x-axis. Indeed it can be seen from Eq. (5) that the incompressible phase field S 0 (r) approaches a constant value (a vacuum) exponentially with a characteristic length scale D/π on either side of the vortex core and the exponential localisation remains valid at all orders of the perturbation expansions (11) and (14) for the density and phase of the compressible superfluid solutions. For this reason, in the presence of a trapping potential, the dynamics of the solitonic vortex can be treated in the framework of Landau quasiparticle dynamics, as outlined in Sec. I and previously considered in Ref. [16].

A. Period of small amplitude oscillations
In the harmonic trapping potential the buoyancy-like restoring force F takes the form and the solutions of Newton's equation (1) are oscillations, as shown in Fig. 4. In the regime of smallamplitude oscillations, the inertial and physical masses M * and M ph can be replaced by their respective values in the center of the channel. Equation (1) now becomes that of a harmonic oscillator with an oscillation period T 0 given by where T trap = 2π/ω x . The coefficients a 2 and a 4 are at most logarithmic functions of ξ 2 /D 2 with explicit expressions as follows: where F ≡ p F q 1, 1, 2 − γ; 2, 3; . The oscillation period is evaluated with the leading-order (incompressible) vortex velocity, and is shown relative to the period of small-amplitude oscillations T0, which is given by the mass ratio at h = D/2. The period Ts depends on the oscillation amplitude, but for amplitudes below 0.1RTF the motion is approximately simple harmonic. The dashed lines in (b) show the analytic expan- (γ = 1) F = 0 and quantum pressure terms are required to obtain the full coefficients. In principle, the coefficients can be evaluated to any order.

B. Trajectories as energy contours
Beyond the limit of small-amplitude oscillations, the oscillation period of the solitonic vortex depends on the amplitude. The trajectories of the solitonic vortex and its oscillation period can be calculated from the free energy E s (µ, V ) in local density approximation (LDA) where µ is replaced by the local chemical potential at the position of the vortex µ(X) = µ 0 − V trap (X). Evaluating the energy expression in the homogeneous slab to leading order yields where we have regularised the integral by excising a disk with radius R as in App. B. In the presence of the harmonic trapping potential by virtue of the LDA the background density n 0 is replaced by where we have made use of the polytropic equation of state, and R TF = 2µ 0 /(mω 2 x ) is the Thomas-Fermi radius. Applying the LDA to the energy expression of Eq. (44) and identifying the cutoff radius R with the healing length ξ, we obtain Instead of the constant f 0 , we can use the turning points to have x max , the amplitude of the oscillations, as our control parameter. Eliminating f 0 using Eq. (46), we obtain an implicit function for y and x in terms of x max that gives the same trajectories:

C. Oscillation period
The oscillation period can be calculated by integration over the trajectory as Figure 4(b) shows the resulting dependence of the oscillation period on the amplitude max obtained from numerical integration, where the velocity V = V 0 + V 1 includes numerically obtained compressible corrections up to O(ξ 2 /D 2 ). An analytic approximation for the oscillation period can be obtained by approximating the velocity by V = V 0 as given in Eq. (17). Integrating over the energy contour (47) and taking the limit x max → 0 reproduces T 0 of Eq. (41) to leading order in ξ 2 /D 2 . Expanding around x max = 0, we obtain The second order expansion of the oscillation period is shown in Fig. 4(b) as dashed lines.

V. DISCUSSION
The mass ratio (41) can become large for a solitonic vortex compared to a dark soliton, which is consistent with the interpretation of recent experiments [16,17] performed in cylindrical, all-harmonic trapping potentials. Future experiments in a rectangular slab geometry with hard-wall potentials could provide quantitative comparison of the oscillation period with the results of this work -potentially at very high accuracy. This could lead to identifying effects that are not yet included in the theory including the Kopnin mass, a mass contribution to the vortex expected to arise from fermionic bound states in the vortex core of a fermionic superfluid [12]. Also the role of quantum pressure contributions to the hydrodynamic equations of Fermi gases in the BEC to BCS crossover is currently unclear and experiments or Monte-Carlo simulations could clarify the situation by comparison to the theory presented here.
Our results for the inertial vortex mass may be compared with the recent discussion of vortex mass contributions by Sonin [12]. Based on physical arguments but without recourse to a specific confining geometry, Sonin introduces the 'core' mass and the 'compressibility' mass. While no mass contribution that scales as the leading contribution to the inertial mass (29) appears in Ref. [12], the 'core' and 'compressibility' masses resemble the next to leading order contributions ∝ mq 2 n 0 ξ 2 to the inertial mass of Eq. (29). The 'core' mass is related to the fluid displacement by the vortex core and is consistent with the non-logarithmic contribution to the ξ 2 term in Eq. (29). It is not possible to compare the terms quantitatively as any constant can be absorbed in the logarithmic term that depends on the geometry. Sonin's 'compressibility' mass has the same scaling as but a different sign than the ξ 2 ln(2D/πξ) term in Eq. (29), which does depend on the confining geometry. We find that the overall sign of this term becomes negative when ξ/D 0.35 whereas both contributions discussed in Ref. [12] as being relevant for bosonic superfluids are positive. It should be noted that the contributions at this level are geometry-dependent and thus it is essential to define the geometry as done with the slab geometry in this paper in order to make quantitative predictions.
The series expansion for the physical and inertial vortex masses established here provides the basis to further improving the theoretical understanding of stronglycorrelated quantum liquids like the superfluid Fermi gas beyond the currently known hydrodynamic model (7).
By measuring the oscillation period T 0 , future highprecision experiments in the slab geometry can determine the expansion coefficients order-by-order and inform the modelling of quantum pressure, transverse forces [4][5][6], and vortex core filling [47]. Moreover, our results can be further developed to make experimental measurements of certain microscopic features such as the level spacing of the Andreev bound states of the vortex core [47].
Acknowledgements: The authors thank Sandy Fetter for inspiring discussions. Both of the authors contributed equally to the conception of the work and to the writing of the manuscript. L.A.T. carried out the detailed calculations and J.B. designed the incompressible model.

Appendix A: Derivation of the Euler equation
Here we outline the derivation of Eq. (7) in the main text. At zero temperature, changes in the chemical potential µ for a bulk system are related to changes in the pressure p by the Gibbs-Duhem relation dp = ndµ/m, where n is the number density. The third law of thermodynamics states s = 0 for a perfect crystal at T = 0. On the other hand, for an isentropic process (ds = 0, where s is the entropy), we have dp = ndw [48], wherew is the enthalpy. Let us therefore identify µ = mw at zero temperature for isentropic processes. The Nernst-Simon formulation of the third law of thermodynamics states that ds → 0 for any reversible isothermal (dT = 0) process as T → 0, i.e. any reversible isothermal process at T = 0 is isentropic.
The Euler equation in Eq. (2.9) of Ref. [48] is given as: Applying a Galilean boost where we have usedv b = 0 and that its gradients also vanish. Importantly, in the boosted frame that travels alongside the vortex, the velocity does not change, i.e. ∂v ∂t = 0. In the main text we calculate v . We obtain (dropping the primes)  (7) of the main text now follows after multiplying with m by replacing mw with the chemical potential at the local density and denoting the constant with µ 0 .
Appendix B: First-order compressible phase correction S1: regularisation and renormalisation 1. Poisson equation for S1 The first-order velocity correction can be obtained from Eq. (14), which amounts to solving Poisson's equation in the domain Ω with zero Neumann boundary conditions (where Ω represents the channel as per Fig. 1a of the main text): where S 1 is the first order phase field, and n 1 = − 2 2m n0 µ0 γ|∇S 0 | 2 is the first-order density correction. Note that the right hand side would be zero for a free vortex since there the velocity field lines ∇S 0 are orthogonal to the density gradient ∇n 1 . For the solitonic vortex this is not the case and thus the right hand side is non-zero.

(B2)
The solution is then found by integrating the Green's function with the source term (B3) where we have regularised the otherwise divergent integral by excising a disk B R ρ with radius R. This integral can be solved numerically, and a renormalisation procedure that allows R to be taken to zero is outlined below. Further analytical progress is made noting that ∇ 2 S 0 = 0 in the integration domain and writing G∇n 1 ·∇S 0 = G∇·(n 1 ∇S 0 ) = ∇·[Gn 1 ∇S 0 ]−n 1 ∇S 0 ·∇G: (B4) The first term on the right hand side can be analytically evaluated as a boundary term on ∂B R ρ , which is the circle where |r − ρ| = R. Setting ρ = (0, h), the position of the vortex, we have (x ) 2 + (y − h) 2 = R 2 , so that x = R cos(θ) and y = h + R sin(θ), where θ is measured at the vortex. We have to integrate over all θ, and obtain as an approximation for S 1 that ignores the second term on the right hand side of Eq. (B4) where r = (x, y). The (R/D) 2 series corresponds to a 'multipole-like' expansion (Fig. 5). We take the limit R → 0 so that in effect we are excising a point.

Renormalisation
At first sight, the integral (B3) for S 1 is divergent. However, as we outline here, our model is renormalisable, and the physical properties of the solitonic vortex can be meaningfully assigned.
First, we regularise the divergence by excising the disk B R ρ . The analytically solvable term S approx 1 equates to a multipole-like expansion in the cut-off R (Fig. 5a,b,c), while the term that was ignored in Eq. (B4) is divergent as R → 0. More specifically, it integrates to a divergent dipole term (Fig. 5d), and an R-independent quadrupole term (Fig. 5f). The dipole terms from both the analytically solvable term and the numerically evaluated term vanish at h = D/2, and S 1 takes the form of a quadrupole field there (Fig. 5f). Since the quadrupole field has no R-dependence in the full S 1 (it has R 2 dependence only in the analytically solvable term S approx 1 ), the integral for S 1 is already renormalised and indeed gives an identical result as the full renormalisation procedure.
To renormalise the divergences, we expand the inte- /(q 3 ξ 2 /D 2 ) to leading order i.e. no R-dependence, (b) next-to-leading order R 2 /D 2 , (c) and order R 4 /D 4 . (d) Numerical evaluation of S1/(q 3 ξ 2 /D 2 ) with an excision (R/D = 0.01) at the vortex core. Note the (divergent) scale from having a small R. (e) Renormalised, cut-off independent S1/(q 3 ξ 2 /D 2 ). In (a)-(e), h/D = 0.7. (f) Renormalised, cut-off independent S1/(q 3 ξ 2 /D 2 ) at h = D/2. grand in Eq. (B4) around the vortex position, and add counter-terms to cancel the terms in the expansion with negative powers ofr, wherer is the distance from the vortex core. We have − 1 n 0 ∇ · (n 1 (r )∇S 0 (r ))G(r, r ) = a 3 r 3 + a 2 r 2 + a 1 r + . . . , (B6) and so we need three counter terms. Here a 1−3 are an-alytic expressions that depend on r and r . Then, we take R → 0 and obtain finite values for the remaining integrals. By doing this we find that S 1 is in general a combination of a dipole term and a quadrupole term (Fig. 5e). As discussed, only the finite quadrupole term survives at h = D/2. Although the integral (B3) for S 1 formally diverges when h = D/2, the underlying physical model is thus renormalisable with only a finite number of counter terms required.