1 Introduction

The interaction between topological defects and their environment is responsible for a host of fascinating phenomena in physical and biological systems. One of the most important examples thereof is the pinning of topological defects, e.g. rotating waves in cardiac muscles [1,2,3], vortices in active matter [4] and nematic defects in liquid crystals [5, 6]. In quantum fluids, such as superfluids and superconductors, the nucleation and motion of quantized vortices play crucial roles. For instance, in superfluid He-II, the presence of vortices can lead to dissipation of the superflow [7]. Additionally, vortices can be readily pinned to obstacles with a length scale comparable to the superfluid healing length [8], such as bumps in a superfluid container and defects in superconductors. Furthermore, vortex pinning results in a correction to the Berezinskii–Kosterlitz–Thouless transition in thin-film He-II [9,10,11], holding magnetic flux in type-II superconductors [12,13,14], an increase in critical counterflow velocity [5]. Vortex pinning can also be implemented in order to manipulate atomtronic devices [15,16,17].

Quantum fluids are also thought to be present in cosmological and astrophysical systems such as dark matter and neutron stars. In the ultra-light dark matter model [18,19,20,21], vortices are found to be unstable in the central region of dark matter halos [22,23,24] but are associated with the granule size with a turbulence-like characteristics in the outer regions [25]. In the interior of a neutron star, both neutrons and protons can be in the superfluid phase [26]. The neutron fluid typically contains of order \(10^{18}\) vortices, which can pin to the nuclear lattice in the star’s outer crust [27, 28] and to magnetic flux tubes in the star’s core [29, 30]. In this system, vortex pinning prevents the superfluid from spinning down at the same rate as the crust, thereby creating a rotational lag. It is believed that, when this lag reaches a critical value, the vortices depin and transfer angular momentum from the superfluid to the crust, resulting in a sudden increase in the observed rotational frequency known as a glitch. The exact process by which a glitch occurs is not fully understood but may involve an avalanche of millions of vortices that depin and thereby trigger further depinning [31]. Given the complexity of the mechanism underlying these spindown glitches, it is pertinent to first study the dynamics of a small number of vortices and thereby understand the conditions under which they depin, which is the aim of the present work.

From the perspective of a pinned vortex, the rotational lag in the neutron star crust manifests as an ambient superflow, which exerts a lateral Magnus force on the vortex [32]. Above a critical superfluid velocity, this force causes the vortex to depin, after which it moves with essentially the same velocity as the ambient flow. This critical velocity was first studied numerically by Schwarz [33] using the vortex filament method (in which a vortex is modelled as a one-dimensional line); in this model, the pinning site was a hemispherical ‘bump’ on the boundary of the superfluid. Subsequently, [8] demonstrated that multiple vortices can be trapped on the same bump. More recently, [34] investigated a related problem—the scattering of a superfluid vortex by an obstacle. For an obstacle of large width, they found that the critical velocity for pinning could be predicted using the equation of motion for a vortex under the small-displacement approximation [35].

Despite this, the depinning process of a vortex subjected to a background superflow still lacks a detailed theoretical investigation. Thus, we carefully study the depinning dynamics of an initially pinned vortex due to the presence of a background superflow in a two-dimensional superfluid, as illustrated in Fig. 1. The superfluid system is modelled by the Gross–Pitaevskii equation with a phenomenological dissipation. The phase diagram of the critical depinning velocity is explored in terms of the height and width of the pinning potential, using an energetic argument to understand the transition between the pinned and free vortex states. A complementary set of simulations is also conducted in a different numerical setup with a different boundary condition to validate our finding.

Fig. 1
figure 1

Schematic plots of a vortex (red dot) subjected to a pinning site of width w and a background superflow \({\textbf{v}}_s=-v_y\hat{{\textbf{y}}}\). Initially (a), the vortex is trapped within the pinning site but (b) depins when \(v_y\) exceeds a critical value \(v_c\). The vortex trajectory is marked by a long triangle with a colour gradient. The colour gradient indicates the vortex positions at different times, and the width of the triangle is the velocity of the vortex

Full size image

This work is structured as follows. In Sect. 2, we describe our model and discuss the effects of dissipation on the vortex motion. We then proceed to address the numerical setup of simulating the dynamics of depinning and the results of these simulations in Sect. 3. Subsequently, in Sect. 4, we propose a phase diagram for the critical depinning velocity and an energy landscape for vortex depinning alongside a validation check of the depinning simulations employing quasiperiodic boundary conditions (QPBCs). Finally, we summarize and discuss our findings in Sect. 5.

2 Theoretical Formalism

2.1 Gross–Pitaevskii Theory

We model a two-dimensional superfluid using the damped Gross–Pitaevskii equation (dGPE) [36,37,38,39,40,41],

$$\begin{aligned} i\hbar \frac{\partial }{\partial t}\psi ({\textbf{r}},t)=(1-i\gamma )\left( {\hat{H}}_\textrm{GP}-\mu \right) \psi ({\textbf{r}},t) \end{aligned}$$
(1)

where \(\psi\) is the order parameter for superfluid particles with mass m, \(\mu\) is the chemical potential and

$$\begin{aligned} {\hat{H}}_\textrm{GP}=-\frac{\hbar ^2\nabla ^2}{2m}+V_\textrm{pin}({\textbf{r}})+g\left| \psi ({\textbf{r}},t)\right| ^2. \end{aligned}$$
(2)

The superfluid number density and velocity are related to \(\psi\) by \(|\psi ({\textbf{r}},t)|^2\) and \({\textbf{v}}_s=({\hbar }/{m})\nabla \textrm{arg}[\psi ({\textbf{r}},t)]\), respectively. The parameter \(g>0\) describes the self-repulsion of the superfluid and is related to the s-wave interaction strength [42]. For the pinning potential, \(V_\textrm{pin}({\textbf{r}})\), we assume a Gaussian profile [43, 44],

$$\begin{aligned} {V_\textrm{pin}({\textbf{r}})=V_0e^{-(x^2+y^2)/2w^2},} \end{aligned}$$
(3)

with height \(V_0\) and width w. The dimensionless parameter \(\gamma\) in Eq. (1) represents dissipation arising from the interaction between superfluid and normal componentsFootnote 1. In the context of Bose–Einstein condensates, \(\gamma\) is proportional to temperature and is usually considered to be spatially constant with a numerical value of \(\lesssim 10^{-3}\) [45,46,47]. However, for reasons described in Sect. 2.2, our model includes a spatially dependent \(\gamma ({\textbf{r}})\) that suppresses waves and other dynamics far from the pinning site.

As described later, we choose the value of \(\mu\) in all of our simulations such that a constant density, \(|\psi |^2=n_0\), is maintained at large distances from the pinning site. The system then has a characteristic length scale given by the healing length, \(\xi _0=\hbar /\sqrt{mgn_0}\), and a characteristic time scale \(t_0=\hbar /gn_0\). The ratio of these determines the speed of sound in the superfluid, \(c_0=\sqrt{gn_0/m}\) [42]. Provided that \(\gamma \ll 1\), dissipation only plays a role on a much longer time scale, \(\tau _0 = t_0/\gamma\) [48, 49].

In this work, we theoretically investigate the critical depinning velocity of a vortex initially pinned by a pinning potential in a superfluid with a background flow \({\textbf{v}}_s=-v_y\hat{{\textbf{y}}}\). The width of the pinning potential, w, is chosen to be of the same order as the healing length, \(\xi _0\), which is the relevant parameter regime for the neutron star crust. Unlike several previous works that assume a wide obstacle, we cannot use either the vortex filament or the Thomas–Fermi (TF) approximation in this regime. Instead, we use the dGPE model to determine the critical velocity. For later convenience, we introduce the following notation to represent different possible states of the system:

  1. (i)

    \(\psi _{0,v_y}=\sqrt{n_0}e^{-imv_yy/\hbar }\): a homogeneous density solution with a background flow;

  2. (ii)

    \(\psi _{\textrm{v},v_y}\): a single vortex state subjected to a background flow;

  3. (iii)

    \(\psi _{\textrm{ps},v_y}\): a vortex-free state with a pinning site subjected to a background flow;

  4. (iv)

    \(\psi _{\textrm{pv},v_y}\): a pinned vortex state subjected to a background flow.

  5. (v)

    \(\psi _{\textrm{fv},v_y}\): a free (depinned) vortex state subjected to a background flow far away from the pinning site.

States (ii)–(v) cannot be fully described analytically, and instead, we study them numerically, focusing on the dynamics in Sect. 3 and the energetics in Sect. 4.2.

2.2 Numerical Setup

Fig. 2
figure 2

a A schematic plot of the initialized states with the positions of vortices and anti-vortices marked in red and blue, respectively. The targeted vortex sits at the centre while the remaining three vortices (one vortex and two anti-vortices) are generated by the phase frustration at the edges due to the periodic boundary condition. The sponge layer, Eq. (4), is represented by the round rectangular contours for \(x^{10}+y^{10}=R_s^{10}\) (grey dashed line), \(x^{10}+y^{10}=(R_s - w_\textrm{abs})^{10}\) (grey dotted line) and \(x^{10}+y^{10}=(R_s+w_\textrm{abs})^{10}\) (grey dashed-dotted line). b Density profiles corresponding to the numerical solutions in the absence of a background superflow for a free vortex (\(n = |\psi _{\mathrm {v/pv},v_y=0}|^2\), \(V_0 = w = 0\)) and a pinned vortex (\(n = |\psi _{\mathrm {v/pv},v_y=0}|^2\), \(V_0=2E_0\) and \(w=\lbrace \xi _0,\,2\xi _0,\,3\xi _0\rbrace\)). The grey dashed line is the density profile of the analytical ansatz, Eq. (17), for the free vortex while the dotted lines are the numerical solutions for \(\psi _{\textrm{ps},v_y=0}\). The hollow circles mark the Thomas–Fermi radius \(R_\textrm{TF} = w\ln (V_0/E_0)\) where the density vanishes for \(r<R_\textrm{TF}\) under the TF approximation, suggesting that the vortex density profiles are poorly described by the TF solution

Full size image

We solve Eq. (1) using a Fourier pseudospectral method and a 4th-order Runge–Kutta scheme with a time step of \(\Delta t=10^{-3}t_0\). The computational domain is a square box of width \(L_x=L_y=256\xi _0\) with \(N_x=N_y=512\) grid points, thereby giving a grid spacing of \(\xi _0/2\). The pinning site, when present, is located at the centre of the box.

The use of Fourier transforms requires periodic boundary conditions in both x and y, and thereby imposes a vanishing net circulation around the boundary of the domain [50]. Although we are only interested in the dynamics of a single vortex, it is, therefore, necessary to simulate an alternating lattice of vortices and anti-vortices as illustrated in Fig. 2a. Given the size of the domain, we expect the effect of these additional vortices to be small, and we minimize their effect further by introducing a sponge layer via a spatially dependent dissipation [40, 51, 52]. We adopt a smooth, rounded-rectangular sponge as plotted in Fig. 2a by defining

$$\begin{aligned} \gamma ({\textbf{r}})=\gamma _0+\frac{\gamma _\text {abs}}{2}\left[ 1+\tanh \left\{ \frac{\left( x^{10}+y^{10}\right) ^{1/10}-R_s}{w_\textrm{abs}}\right\} \right] \end{aligned}$$
(4)

where \(\gamma _0\) is the dissipation induced by the interaction between superfluid and normal components. The strong dissipation set by \(\gamma _\text {abs}\) at the boundaries of the domain damps any sound waves, thereby isolating the central vortex from its periodic neighbours. However, this dissipation also induces a transverse motion of the vortices near the edge of the domain, with a velocity \(-\gamma \kappa v_y\hat{{\textbf{x}}}\) [17, 35]. Fortunately, this transverse motion ends once the vortices reach the edge of the sponge layer. Thus, for a simulation domain large enough that the sponge layer is sufficiently distant from the central vortex, the effects of the edge vortices are marginal. Here, we set \(\gamma _\textrm{abs} = 1\), \(R_s=98\xi _0\) and \(w_\textrm{abs}=10\xi _0\).

Another consequence of our periodic boundary conditions is that the imposed superflow velocity, \(v_y\), must be an integer multiple of \(\Delta v = 2\pi c_0\xi _0/L_y \approx 0.025c_0\). This is because the phase of the order parameter, \(\arg \psi\), must also be periodic, and it is the gradient of the phase that determines the superflow. This restriction limits the precision with which the critical velocity can be determined but, given the size of the domain, this does not represent a significant limitation of our methodology.

2.3 Initialization

As noted earlier, we specify the value of the chemical potential \(\mu\) in order to fix the (background) density to \(n_0\). The correct choice is not only dependent on g and \(n_0\) but also on the background superflow velocity \({\textbf{v}}_s=-v_y\hat{{\textbf{y}}}\). In the dGPE simulation, the dissipative term ultimately drives the system towards a state satisfying \(({\hat{H}}_\textrm{GP}-\mu )\psi \approx 0\), and so we must choose the value of \(\mu\) such that in this steady state, the density is equal to \(n_0\) far from the pinning site. To determine the appropriate value, we consider the grand-canonical energy of the system,

$$\begin{aligned} F=\int _A d^2{\textbf{r}}\psi ^*({\textbf{r}})\left[ -\frac{\hbar ^2\nabla ^2}{2m}+V_\textrm{pin}({\textbf{r}})+\frac{g}{2}\left| \psi ({\textbf{r}})\right| ^2\right] \psi ({\textbf{r}})-\mu N_A \end{aligned}$$
(5)

where A is the domain area, g is the s-wave interaction strength and \(N_A=\int _Ad^2{\textbf{r}}\left| \psi ({\textbf{r}})\right| ^2\) is the particle number in the domain. The damped Gross–Pitaevskii Eq. (1) can be expressed as \(i\hbar \partial _t\psi =(1-i\gamma )\delta F/\delta \psi ^*\), and the steady state corresponds to a minimum of the the functional F. In the absence of any vortex or pinning site, we expect the steady state to be \(\psi =\psi _{0,v_y} = \sqrt{n_0}e^{-imv_yy/\hbar }\), and this is a minimum of Eq. (5) provided that

$$\begin{aligned} \mu =\mu _{0,v_y}=\frac{mv_y^2}{2}+E_0, \end{aligned}$$
(6)

with a background energy \(E_0=gn_0\). Therefore, \(\mu\) is set to the value given by Eq. (6) in our simulations.

To initialize the system with a vortex (with or without a pinning site), we first set \(v_y=0\) and imprint a phase winding by setting \(\psi = \sqrt{n_0}e^{i\phi }\), where \(\phi\) is the polar angle relative to the centre of the domain. The GPE is then solved in imaginary time (i.e. with \((1-i\gamma )\) in Eq. (1) replaced by -1), with \(\mu = E_0\), until the quantity

$$\begin{aligned} \frac{1}{N_A}\int _A d^2{\textbf{r}}\,\psi ^*{\hat{H}}_\text {GP}\psi \end{aligned}$$
(7)

converges to within \(10^{-7}\) of \(E_0\). This converged state is then used as the initial condition for the dGPE. As mentioned earlier, because of the periodic boundary conditions, this state actually features a lattice of vortices and anti-vortices but, for a sufficiently large domain, the central vortex is essentially independent of the others.

As examples of the solutions of Eq. (1), in Fig. 2b, we plot the density profiles of the free (\(\psi _{\textrm{v},v_y=0}({\textbf{r}})\)) and pinned (\(\psi _{\textrm{pv},v_y=0}\)) vortices (both represented by solid lines), as well as the vortex-free state subject to the pinning potential (\(\psi _{\textrm{ps},v_y=0}\), represented by dotted lines), for \(V_0 = 2E_0\) with different pinning widths w. The vortex solutions show a density depletion with a width similar to the pinning site but it is also evident that they do not agree well with the Thomas–Fermi density,

$$\begin{aligned} n_\textrm{TF}({\textbf{r}})= & {} n_0[1-V_\textrm{pin}({\textbf{r}})/E_0]\Theta (|{\textbf{r}}|-R_\textrm{TF}), \end{aligned}$$
(8)
$$\begin{aligned} R_\textrm{TF}= & {} w\sqrt{\ln V_0/E_0}, \end{aligned}$$
(9)

as the radii within which the densities are strongly depleted are not in agreement with the TF radii, \(\lbrace R_\textrm{TF}\rbrace\), which we have presented in Fig. 2b as hollow circles. This suggests that the pinning potential height and width we have specified are too small for a vortex to be in the TF regime [34].

3 Depinning Dynamics

With the initial state prepared, we impose a background flow by introducing a phase gradient

$$\begin{aligned} \psi _{\textrm{pv},v_y}({\textbf{r}},t=0)=\psi _{\textrm{pv},v_y=0}({\textbf{r}})e^{-imv_yy/\hbar }, \end{aligned}$$
(10)

where \(v_y\) must be an integer multiple of \(\Delta v=2\pi c_0\xi _0/L_y\), for reasons explained earlier. At the same time, we update the value of \(\mu\) to \(\mu = E_0+mv_y^2/2\) so that the system remains in a steady state far from the vortices and pinning site. The simulations are performed with three values of dissipation: \(\gamma _0=0\), \(5\times 10^{-4}\) and \(5\times 10^{-3}\). We find that the value of \(\gamma _0\) has very little effect on the dynamics, including the depinning process, except that after depinning the vortex drifts with respect to the ambient superflow by an amount proportional to \(\gamma _0\). This transverse drift is expected for reasons mentioned in Sect. 2.2. Although the total energy and momentum are not conserved quantities due to the presence of dissipation, in the simulations, we find that they only vary by a few per cent.

Fig. 3
figure 3

a Examples of vortex trajectories for \(V_0=2E_0\), \(w=\sqrt{2}\xi _0\) and \(\gamma _0=5\times 10^{-3}\), with (i) \(v_y\approx 0.172c_0\), (ii) \(v_y\approx 0.196c_0\) and (iii) \(v_y\approx 0.221c_0\) within the time span \(t\in (0,200t_0)\). b–g The density fluctuations, measured by \(\delta n=[n({\textbf{r}},t)-n_\textrm{ps}({\textbf{r}})]/n_0\), are plotted when b \(t=0\), c \(t=6t_0\), d \(t=14t_0\), e \(t=40t_0\), f \(t=60t_0\) and g \(t=120t_0\), where it is evident that the radiated sound waves spiral out from the vortex. The vortex moves to the left of the pinning site due to the Magnus force but, once it escapes from the pinning site, it drifts downward with the superflow at a speed \(v_y\). The width of the pinning potential is indicated by the grey dashed circle. No significant difference between the dynamics for \(\gamma _0=0\), \(5\times 10^{-4}\) and \(5\times 10^3\) can be noted

Full size image

Three types of dynamics are found in our simulation: (1) If \(v_y\) is sufficiently small, then the vortex is displaced a small distance from the centre of the pinning site, but remains pinned; (2) if \(v_y\) exceeds a critical value, \(v_c\), the vortex depins and is carried away by the ambient superflow and (3) if \(v_y\) is a significant fraction of the sound speed, \(c_0\), additional vortices nucleate at the pinning site. However, the present work is not concerned with regime (3), which has been studied in detail by others, e.g. Refs. [43, 44]. Figure 3a illustrates typical vortex trajectories in regimes (1) and (2) for \(V_0=2E_0\) and \(w=\sqrt{2}\xi _0\). The vortex position in the figure is determined by locating the phase defect and density minimum in the wavefunction, after interpolating to the sub-grid level, and the trajectory is subsequently tracked by linking timeframes using the Hungarian algorithm. In Fig. 3, the orange–black colour transition of the trajectories tracks the evolution of the vortex position from \(t=0\) to \(200t_0\). It is evident that, in regime (1), the vortex follows a spiral trajectory and eventually reaches an equilibrium position to the left of the pinning site centre. In this new equilibrium, the pinning force \(\propto \hat{{\textbf{x}}}\) balances the Magnus force \(\propto \hat{{\textbf{z}}}\times (-v_y\hat{{\textbf{y}}})\) [32, 35]. The maximum spiral radius increases with \(v_y\), as shown in Fig. 3a (i) and (ii) for \(v_y\approx 0.172c_0\) and \(0.196c_0\), respectively. Conversely, if \(v_y > v_c\approx 0.221c_0\) the system is in regime (2); the vortex initially follows a similar trajectory but moves far enough from the pinning site centre that it escapes and ultimately drifts with the ambient superflow.

In both regimes (1) and (2), the motion of the vortex excites sound waves that carry energy away from the vortex. To illustrate this, in Fig. 3b–g, we plot the density fluctuation \(\delta n=[n({\textbf{r}},t)-n_\textrm{ps}({\textbf{r}})]/n_0\), where \(n_\textrm{ps}({\textbf{r}})=|\psi _{\textrm{ps},v_y=0}|^2\) is the density in the absence of a vortex or flow, together with the vortex trajectory for \(t/t_0=0\), 6, 14, 40, 60 and 120. The plots show the case with \(\gamma _0=5\times 10^{-3}\), but we observe very similar results for the other values of \(\gamma _0\) considered. In these plots, the vortex location is marked by a red circle, and the orange–black line is the trajectory up to time t. Similarly to the case of a vortex in a stirred condensate [53], while the vortex remains close to the pinning site, the sound waves form a dipolar pattern that spirals out from the vortex. For comparison, we have also performed a simulation with the same parameters but without a vortex; in that case, sound waves are still produced by the sudden imposition of a superflow past the pinning site at \(t=0\), but these sound waves have a roughly circular pattern, and the density perturbations are smaller in magnitude by a factor of about 2.

The emission of sound waves becomes negligible in the later time dynamics whether the vortex remains pinned or depins, and those emitted at early times rapidly dissipate within the sponge layer around the domain boundary. We note that, as shown in Fig. 3g, density perturbations persist around the pinning site even after the vortex has depinned. This reflects the effect of the ambient flow on the pinning site. Furthermore, as mentioned earlier, at sufficiently high-flow velocities, vortices are nucleated in the pinning potential [43, 44]. If the pinning potential is too high (\(V_0\ge 3.5E_0\)) or too wide (\(\sqrt{2}w\ge 3.5\xi _0\)), then we find that this regime is reached before the superflow velocity reaches the critical depinning velocity \(v_c\). In this regime, the depinning process is complicated by the involvement of multiple vortex interactions [34], and therefore, we only present data for pinning sites with smaller values of \(V_0\) and w.

4 Critical Depinning Velocity

4.1 Phase Diagram of the Depinning Velocity

As illustrated in Sect. 3, if the imposed flow, \(v_y\), is not sufficient to depin the vortex, then it generally settles to a new equilibrium position within a duration of \(200t_0\). This new equilibrium is generally within a distance w of the pinning site centre. In our results, we, therefore, consider the vortex to be depinned if it is displaced by more than 1.5w from the pinning site within a duration of \(200t_0\), and we define \(v_c\) to be the smallest value of \(v_y\) for which depinning occurs. (We recall that our periodic boundary conditions only allow us to increase \(v_y\) in steps of \(\Delta v\approx 0.025c_0\). Therefore, our results for \(v_c\) may overestimate its true value by up to \(\Delta v\).) In Fig. 4a, we show the values of \(v_c\) obtained for \(V_0/E_0\in [0.5,3.5]\) and \(w/\xi _0\in [1,3.5]\). This figure shows the results in the cases \(\gamma _0 = 0\) and \(\gamma _0 = 5\times 10^{-4}\), which are identical; the results for \(\gamma _0 = 5\times 10^{-3}\) are nearly identical, except for the case \(V_0=0.5E_0\) and \(w=\xi _0\), where we find \(v_c\) to be larger by \(\Delta v\). This indicates that dissipation plays a negligible role in determining the critical velocity in our simulations.

Fig. 4
figure 4

a Phase diagram of the critical velocity \(v_c\) for \(\gamma _0=0\) and \(5\times 10^{-4}\). b Comparison between the critical velocity and the empirical formula, Eq. (13). The lower errorbar indicates the possible range for \(v_c\) given the velocity resolution, \(\Delta v\approx 0.025c_0\)

Full size image

From dimensional analysis, the value of \(v_c/c_0\) ought to depend only on \(V_0/E_0\), \(w/\xi _0\) and \(\gamma _0\). Thus, for sufficiently small \(\gamma _0\), we ought, therefore, to be able to fit our results to a function of the form \(v_c/c_0 = f(V_0/E_0,w/\xi _0)\). Before performing this fitting, we first briefly review earlier work concerning the form of the function f. Stockdale et al. [34] used the equation of motion of a vortex under the small-displacement approximation [35], in the absence of dissipation, to investigate the pinning of a vortex in 2D with a circular top-hat barrier and a background flow. Making a Thomas–Fermi approximation for the density, they found that the vortex remains pinned if its equation of motion has a fixed point within the pinning site. In terms of our notation, the critical depinning velocity is then given by

$$\begin{aligned} v_c/c_0=\frac{\xi _0}{4\pi }\left( \frac{2w+\delta }{2w\delta }\right) \frac{V_0}{E_0} \end{aligned}$$
(11)

where w is the radius of the obstacle, and \(\delta\) is a phenomenological screening parameter of order \(\xi _0\). This result is only valid in the regime where \(w > \delta\) and \(V_0 < E_0\), however, which is not the regime of interest of the present work.

Earlier, [33] used a vortex filament model to simulate a vortex pinned to a hemispherical bump on the boundary of the domain. Here, the critical depinning velocity was found to be

$$\begin{aligned} v_c/c_0\approx \frac{\xi _0}{4\pi D}\ln \frac{w}{\xi _0} \end{aligned}$$
(12)

where w is the radius of the obstacle, and D is the channel width of the superfluid container. Again, this result assumes \(w \gg \xi _0\) and is intrinsically three-dimensional, and thus, it is not directly applicable to our results.

Motivated by the above literature, and by the results in Fig. 4a, we seek a fit of the critical velocity to the form

$$\begin{aligned} v_c/c_0=f\left( \frac{V_0}{E_0},\frac{w}{\xi _0}\right) =a\ln \left[ 1+ \left( b\frac{w}{\xi _0}+c\right) \frac{V_0}{E_0}\right] \end{aligned}$$
(13)

where the parameters a, b and c are assumed to be constant. Applying a least-squares fit to our data, we find

$$\begin{aligned} a = 0.1039\pm 0.0175, \; {b = 1.0721 \pm 0.5736} \; \text{ and } \; c = 1.188\pm 0.4880, \end{aligned}$$
(14)

and for these coefficients, the empirical formula predicts \(v_c\) to within \(15\%\) across our data set. In Fig. 4b, we plot \(v_c\) as a function of w for three values of \(V_0\), alongside our empirical formula. The result shows that in our considered range of parameters, the logarithmic trend is not strong, and it would be interesting to investigate a wider range of parameters to better constrain this fit. In particular, for larger values of w, we expect ghost vortices to be generated within the pinning site [34], which would provide a more stringent test of the formula.

4.2 Vortex Energy

In our simulations, energy is injected into the system by a sudden increase in the superflow from 0 to \(v_y\). Part of this energy is converted into sound, and part may be used to depin the vortex. In general, we would not expect the vortex to depin unless doing so reduces the overall energy in the system, and this fact can potentially be used to determine a lower bound on the critical velocity \(v_c\). Specifically, we seek to determine the value of \(v_y\) for which the energy of a pinned vortex state exceeds that of a depinned vortex state. This energy can be defined as the energy cost \(\alpha _{\psi }\) to create the targeted state from a relevant reference state, which is conventionally chosen to be the ground state of the free energy function given by Eq. (5). For a free vortex state, \(\psi _{\mathrm {v,v_y}}\), this is equivalent to the energy of a vortex state, subtracted by the vortex-free state energy, for a fixed number of \(N_A\) superfluid atoms in a domain of area A [54]. This complexity can be reduced by carrying out the calculation of the grand-canonical energy, Eq. (5), at a fixed chemical potential instead [42, 55]. In this case,

$$\begin{aligned} \alpha _{\text {a},v_y}=F[\psi _\text {a}]-F[\psi _\text {gs}] \end{aligned}$$
(15)

where \(\psi _\text {a}\) represents the target state, and \(\psi _\text {gs}\) is a suitable reference groundstate. According to the aforementioned energetic argument, a lower bound on \(v_c\) can be determined by finding the value of \(v_y\) above which the excess energy of a free vortex state far away from the pinning site, \(\alpha _{\textrm{fv},v_y} = \alpha _{\textrm{v},v_y} + \alpha _{{\textrm{ps},v_y}}\), is less than that of a pinned vortex, i.e.

$$\begin{aligned} \alpha _{\textrm{fv},v_y}\le \alpha _{{\textrm{pv},v_y}} \end{aligned}$$
(16)

for \(v_y\ge v_c\). Here, we specifically consider A to be the area of a circular disc with radius R, giving \(A=\pi R^2\), and a vortex or pinning potential is positioned at the centre of the disc.

The creation energy of a free vortex can be estimated by the required energy added into a homogeneous superfluid, namely \(\alpha _{\textrm{v},v_y}=F[\psi _{\textrm{v},v_y}]-F_0\) with \(F_0=(mv_y^2-E_0)N_A/2\). The wavefunction of a free \(\kappa\)-charged vortex can be well approximated by an analytical profile [42, 54, 56],

$$\begin{aligned} \psi _{\textrm{v},v_y}=\sqrt{n_0\frac{r^2/\xi ^2}{1+r^2/\xi ^2}}e^{i\kappa \theta }e^{-imv_yy/\hbar } \end{aligned}$$
(17)

where \(r=\sqrt{x^2+y^2}\) and \(\xi\) is the radius of the vortex core. In the limit \(R\gg \xi\), the leading order terms give

$$\begin{aligned} \alpha _{\textrm{v},v_y}\approx \frac{\pi \hbar ^2\kappa ^2n_0}{m}\ln \frac{R}{\xi }+\frac{\pi \hbar ^2 n_0}{4m}+\frac{\pi gn_0^2\xi ^2}{2}, \end{aligned}$$
(18)

and the minimization of \(\alpha _{{\textrm{v},v_y}}\) with respect to \(\xi\) gives \(\xi =\left| \kappa \right| \xi _0\) ,Footnote 2. This solution, presented in Fig. 2b as a grey dashed line, agrees well with the corresponding numerical solution represented by the black unbroken line in the same figure. The free vortex energy \(\alpha _{\textrm{v},v_y}\) can then be simplified as,Footnote 3

$$\begin{aligned} \alpha _{{\textrm{v},v_y}}\approx \displaystyle \frac{\pi \hbar ^2n_0}{m}\ln 2.12\frac{R}{\xi _0} \end{aligned}$$
(19)

For a singly-charged vortex, \(\kappa =\pm 1\), the numerical result of \(\psi _{\textrm{v},v_y=0}\) agrees well with this prediction as demonstrated in Fig. 6. We note, however, that the \(v_y\)-dependent terms in this analytical calculation cancel each other out and cannot provide a further prediction of \(v_c\). Furthermore, we reiterate that due to the healing length-sized pinning potential, the Thomas–Fermi density profile is not a particularly good approximation for the vortex state as shown in Fig. 2b.

Fig. 5
figure 5

Condensate density profiles of \(\psi _{\textrm{pv},v_y}\) for \(V_0=2E_0\) and \(w=\sqrt{2}\xi _0\) with a \(v_y=0\) and b \(v_y=0.221c_0\) and c that of \(\psi _{\textrm{v},v_y}\) for \(v_y=0.221c_0\). Panel (i) presents the two-dimensional profile while panel (ii) shows density slices along the x and y axes with the green dashed lines in (i) and (ii) representing the pinning potential. These clearly depict an axial asymmetry of the density profiles due to the background flow. The width of the pinning potential is indicated by the green dashed circle and vertical dashed lines in (i) and (ii), respectively

Full size image

To evaluate the vortex energy with a background flow, we numerically obtain \(\psi _{\textrm{pv},v_y}({\textbf{r}})\) and \(\psi _{\textrm{v},v_y}({\textbf{r}})\) by propagating an initial trial vortex state, Eq. (17), in imaginary time with \(\kappa =1\). When the flow amplitude is small, we found that \(|\mu _\textrm{num}-\mu _{0,v_y}|\) can reduce to less than \(10^{-7}\) quickly. In contrast for large \(v_y\), the convergence is much slower, and the central and boundary-induced vortices would annihilate as \(\tau \rightarrow \infty\). In addition, when \(v_y\) is large enough, the vortex is detached from the pinning potential. Hence, we reduce the numerical tolerance so that imaginary time propagation is run until \(|\mu _\textrm{num}-\mu _{0,v_y}|<10^{-6}\). We note, due to this relaxation of numerical tolerance, this state might be merely metastable. Indeed, even for \(v_y>v_c\), it is possible to find pinned vortex states, and some of the vortices may remain pinned after we move to real-time evolution.

In Fig. 5, we illustrate the density profiles of \(\psi _{\textrm{pv},v_y}\) for \(v_y=0\) and \(0.221c_0\) with \(V_0=2E_0\) and \(w=2\xi _0\) as well as that of \(\psi _{\textrm{v},v_y}\) with \(v_y\approx 0.221c_0\) for the sake of comparison. When the superflow velocity is nonzero, the density becomes locally axially asymmetric about the vortex and pinning site for both \(\psi _{\textrm{pv},v_y}\) and \(\psi _{\textrm{ps},v_y}\). The latter case can be observed in the plot of \(\delta n\) in Fig. 3g where a quadrupolar density fluctuation profile is evident.Footnote 4 The combination of the background and vortex flows, and the interaction with the pinning potential locally affects the density and has a maximum influence along the direction orthogonal to the flow. Therefore, this effect is minimized parallel to the flow in \(\psi _{\textrm{pv},v_y}\), as shown in Fig. 5b and c, where it is apparent that n(0, y) remains almost symmetric along the y-axis. We note that during the propagation of Eq. (1) in real time over timescales comparable to the characteristic relaxation timescale, \(\hbar /\gamma gn_0=200t_0\), the stability of these vortex states is ambiguous even when \(v_y>v_c\). This indicates that a depinning process requires additional energy or instabilities compared to those considered here.

Fig. 6
figure 6

The vortex energies for a pinned vortex \(\alpha _{\textrm{pv},v_y}\) and free vortex \(\alpha _{\textrm{fv},v_y}=\alpha _{\textrm{v},v_y}+\alpha _{\textrm{ps},v_y}\) as a function of \(v_y\) for \(V_0=2E_0\) and \(w=\sqrt{2}\xi _0\). Here \(\alpha _{\textrm{v},v_y}\) and \(\alpha _{\textrm{ps},v_y}\) are the vortex energies in a homogeneous condensate without a pinning potential and the excess energy to create a flow in a condensate with a present of a barrier, respectively. The green dotted line plots the analytical result of a free vortex, Eq. (18). The vertical dotted line marks the energetically preferred \(v_c\), and the grey dashed-dotted lines are the values of \(v_c\) obtained from the dGPE simulation with a background superflow, while the dashed lines are the values of \(v_c\) obtained from the complementary, quasiadiabatic advective simulations discussed in Sect. 4.3

Full size image

We have also computed \(\alpha _\psi\) from the numerical simulations for a domain given by a disc, centred on the pinning site, with an area \(A=\pi R^2\) where \(R=50\xi _0\). In Fig. 6, we show the numerical results of \(\alpha _{{\textrm{fv},v_y}}\) and \(\alpha _{\psi _{\textrm{pv},v_y}}\) together with \(\alpha _{{\textrm{v},v_y}}\) and \(\alpha _{{\textrm{ps},v_y}}\) for \(V_0=2E_0\) and \(w=2\xi _0\) as an example. Here, one can see that the pinned vortex energy \(\alpha _{{\textrm{pv},v_y}}\) is approximately constant while the free vortex energy shows a stronger dependence on \(v_y\) and decreases monotonically as a function of \(v_y\). We also note that the value of \(\alpha _{\textrm{v},v_y=0}\) agrees with the analytical calculation presented earlier. Given that \(\alpha _{\textrm{v},v_y=0} > \alpha _{\textrm{pv},v_y=0}\) and \(d\alpha _{\textrm{v},v_y}/dv_y > d\alpha _{\textrm{pv},v_y}/dv_y\), there exists a value of \(v_y\) where \(\alpha _{\textrm{fv},v_y} = \alpha _{\textrm{pv},v_y}\), and Eq. (16) is satisfied. Above this velocity, we expect that it is energetically favourable for the vortex to depin, with the black vertical dotted line in Fig. 6 showing the estimated critical velocity for the set of parameters. However, what we observe is that the value of \(v_c\) based on energetic considerations alone is found to be less than \(v_c\) in the dGPE simulations of Sect. 3 by roughly \(20\%\) (note the grey vertical dot-dashed line in the figure). Such caveats notwithstanding, we argue that this estimate of \(v_c\) based on purely energetic arguments offers a robust order of magnitude estimate. The small difference between the two values points to additional physics playing a role in the process of depinning, which could be the focus of future studies.

4.3 Quasiadiabatic Advective Acceleration

For reasons explained in Sect. 2.2, the use of periodic boundary conditions has consequences for our numerical simulations. For instance, as demonstrated in Fig. 2a, the requirement that \(\psi\) be circulation-free inside the computational cell results in a vortex lattice of alternating circulations when a vortex is imprinted during the imaginary time propagation of Eq. (1). Crucially, the periodicity of \(\psi\) also restricts the allowed values for the imposed flow velocity \(v_y\). To test the robustness of our results, we have, therefore, also taken a complementary approach, using a different numerical scheme to slowly accelerate the imposed flow from zero to a desired value, after which its value is held fixed. To do so, we first perform a Galilean transformation of the Gross–Pitaevskii equation into the frame moving with velocity \(v_y\hat{{\textbf{y}}}\), hereafter referred to as the advective Gross–Pitaevskii equation (aGPE):

$$\begin{aligned} i\hbar \left( \frac{\partial }{\partial t} - v_y\frac{\partial }{\partial y}\right) \psi ({\textbf{r}},t)=\left[ -\frac{\hbar ^2\nabla ^2}{2m}+V_\textrm{pin}({\textbf{r}})+g\left| \psi ({\textbf{r}},t)\right| ^2-\mu \right] \psi ({\textbf{r}},t). \end{aligned}$$
(20)

We solve this equation together with quasiperiodic boundary conditions (QPBCs),

$$\begin{aligned} \psi (x + L_x, y)= & {} \exp \left[ \frac{i\pi }{L_y}\left( y + \frac{L_y}{2}\right) \right] \psi (x, y), \end{aligned}$$
(21)
$$\begin{aligned} \psi (x, y + L_y)= & {} \psi (x, y), \end{aligned}$$
(22)

that account for the phase winding arising from a single vortex within the domain. Initially introduced in the context of superfluid vortices in Ref. [57] and subsequently applied in Refs. [58,59,60], QPBCs have proven to be advantageous for the study of vortex configurations of nonzero net circulation in a homogeneously uniform background fluid. Because the velocity \(v_y\) is slowly increased in this new setup, we anticipate that little sound will be produced, and therefore, we do not include any dissipation in the model. It is for this reason that we refer to this setup as “quasiadiabatic,” in contrast with the instantaneous acceleration of the flow considered in Sect. 3.

Fig. 7
figure 7

Vortex trajectories of the aGPE simulation with a quasiadiabatic acceleration of \(v_y(t)\) from zero till a final value a \(v_{y,f}= 0.174c_0\), b \(v_{y,f}=0.196c_0\), c \(v_{y,f}=0.221 c_0\) and d \(v_{y,f}=0.245c_0\), followed by a duration of evolution at constant \(v_y\), for \(V_0 = 2E_0\) and \(w = \sqrt{2}\xi _0\). The hollow squares mark the vortex position at the end of the advective acceleration, and the red circle shows the position of vortex at \(t=400t_0\). The width of the pinning potential is indicated by the grey dashed circle

Full size image

Initially, we obtain the stationary state with a vortex at the origin by imprinting a superfluid phase Footnote 5 [61, 62]

$$\begin{aligned} \arg \left\{ \vartheta _3\left[ i\pi \left( \frac{x}{L_y} + \frac{L_x}{L_y}\right) - \pi \left( \frac{y}{L_y} + \frac{1}{2}\right) ; e^{-\pi \left( \frac{L_y}{L_x}\right) }\right] \right\} \end{aligned}$$
(23)

and solve the aGPE in imaginary time to find the ground state. Subsequently, the dynamics of this vortex are simulated by propagating Eq. (20) forward in time with a slow acceleration of the advection as given by

$$\begin{aligned} v_y(t)=\left\{ \begin{array}{ll} \displaystyle \frac{\Delta v}{25t_0}t &{} t\le 25n_f, \\ \\ \displaystyle v_{y,f}=\Delta vn_f, &{} t > 25n_f, \end{array}\right. \end{aligned}$$
(24)

where the final advection index \(n_f\in {\mathbb {Z}}\) sets up the advective flow after the acceleration, and the same spatial grids and timestep used are the same as those in Sect. 4.1. The advective flow velocity takes \(25n_ft_0\) to reach the final value \(v_y\) according to Eq. (24), a duration we believe to be sufficient for the vortex to adjust adiabatically to a finite advection. Subsequently, the aGPE is propagated up to \(t = 400t_0\) to allow time for the vortex to depin, à la Fig. 3. This numerical setup is designed to inject as little energy as possible into the system during the acceleration of the flow, and therefore, we expect the critical depinning velocity, \(v_c\), that we obtain to be an upper bound for \(v_c\) in any other numerical setup.

Throughout the quasiadiabatic advective simulations, we use a plaquette-based subgrid interpolation method [63] to locate the vortex as precisely as possible and track the evolution of its position. In Fig. 7, the vortex trajectories of the aGPE simulation for \(V_0 = 2E_0\) and \(w = 2\xi _0\) with \(n_f=7,\;8,\;9,\;10\) are shown. From this, we identify the lowest \(v_{y,f}\) that dislodges the initially pinned vortex as the critical velocity. Similar to the simulations presented in Sect. 3, we see that the vortex initially always drifts to the left of the pinning site centre as a consequence of the Magnus effect. Once \(v_y\) reaches its final value, \(v_{y,f}\), if this value is sufficiently small, then the vortex quickly inspirals to a new equilibrium position; this final equilibrium position is comparable to that in Sect. 3. However, if \(v_{y,f}\) exceeds a critical velocity, \(v_c\), then the vortex is able to escape the pinning potential, as illustrated in Fig. 7d. With this quasiadiabatic setup, we still observe the emission of a spiral pulse of sound waves, but their amplitude is about an order of magnitude smaller than in the instantaneously imposed flow simulations. The critical velocity obtained via this approach is comparable to that found in Sect. 3, but somewhat larger, as illustrated by the grey dashed vertical line in Fig. 6.

5 Conclusion

In these proceedings, we have numerically simulated how an imposed superflow can detach a vortex that is initially pinned to an obstacle of size comparable to the superfluid healing length. We find that, below a critical velocity, the Magnus force causes the vortex to drift laterally with respect to the imposed flow before ultimately spiralling into a new equilibrium position. The motion of the vortex is responsible for the emission of a spiral of sound which contributes to the loss of energy and stabilization of the vortex. If the imposed superflow exceeds a critical value, then the vortex escapes the pinning region and is carried away by this superflow. The critical depinning velocities obtained from dGPE simulations are presented in a phase diagram as a function of the height and width of the Gaussian potential which defines the pinning site. Inspired by theoretical predictions existing in the literature [33, 34], we propose a new empirical formula for the critical depinning velocity which is in quantitative agreement with our numerical findings. When the pinning potential is sufficiently high and wide, vortices can be nucleated in the pinning potential by the imposed flow.

We then provide an energetic argument for a vortex transitioning from a pinned state to a free (depinned) state by evaluating the vortex nucleation energy under different conditions. The free vortex energy can be analytically evaluated by an axially symmetric trial solution, but the flow dependence in energy is not characterized in this Ansatz. Therefore, we utilize the imaginary time propagation method to search for pinned vortex states which are nonexistent when the flow is too rapid. The pinned states are asymmetric, not only in the off-axis location of the vortex in the pinning potential but also a flow-induced asymmetry of the density which is nonetheless small in magnitude. We find that these solutions are not always stationary over long periods of imaginary time propagation but still provide sensible estimations of the pinned states, suggesting that it might be prudent to consider improved numerical methods in future investigations [55, 64]. Proceeding to examine the energy landscape of free and pinned vortices, we observe that the energy of a vortex decreases monotonically as a function of superflow velocity. However, while the energy of a pinned vortex is lower than that of a free vortex for low velocities, its rate of decrease as a function of the superflow velocity is also correspondingly lower. This results in a crossover superfluid velocity above which free vortices are energetically preferable, a velocity that is comparable to but generally lower than the critical velocities found by direct dGPE simulations with an imposed superflow.

Lastly, we verify our finding with a set of complementary simulations where quasiperiodic boundary conditions are employed to describe a vortex subject to a pinning site being accelerated almost adiabatically to the desired final value of the superflow velocity. These studies, conducted in the reference frame advected with the superflow, provide slightly higher estimates for the critical velocity. Together with our analysis of the energy landscape of pinned and free vortices, this suggests that vortex depinning requires additional energy or instabilities aside from effects arising at the boundaries.

In conclusion, this work develops new insights into the depinning dynamics in a pure 2D superfluid system from a narrow obstacle that is amenable to experimental investigations of superfluids and applicable to theoretical studies of neutron stars. This opens the door to combining studies of the influence of a background superflow with that of ambient sound waves in depinning a vortex from a pinning site. For instance, it has been proposed that the sound waves emitted by a moving vortex around the pinning site can induce the depinning of vortices in its vicinity [65]. Thus, we believe that a natural direction for research in this field would be the study of the vortex–vortex and vortex–sound scattering processes in a system of multiple, initially pinned, vortices in the presence of superflow.