Critical Velocity and Arrest of a Superfluid in a Point-Like Disordered Potential

Superfluid flow past a potential barrier is a well studied problem in ultracold Bose gases, however, fewer studies have considered the case of flow through a disordered potential. Here we consider the case of a superfluid flowing through a channel containing multiple point-like barriers, randomly placed to form a disordered potential. We begin by identifying the relationship between the relative position of two point-like barriers and the critical velocity of such an arrangement. We then show that there is a mapping between the critical velocity of a system with two obstacles, and a system with a large number of obstacles. By establishing an initial superflow through a point-like disordered potential, moving faster than the critical velocity, we study how the superflow is arrested through the nucleation of vortices and the breakdown of superfluidity, a problem with interesting connections to quantum turbulence and coarsening. We calculate the vortex decay rate as the width of the barriers is increased, and show that vortex pinning becomes a more important effect for these larger barriers.


I. INTRODUCTION
A prototypical study of turbulence in fluids is that of the wake behind a cylinder in a flow [1].In classical fluids, the degree of turbulence in the flow can be encoded by the dimensionless Reynolds number Re = vD/η, where v is the velocity of the uniform flow, D is the size of an obstacle in the flow, and η is the kinematic viscosity of the fluid.Dynamical similarity allows us to map flows with different v, D and η to the same flow pattern, so long as the combination vD/η is the same.In a superfluid flow, although η → 0, it has been shown that quantum fluids exhibit dynamic similarities in the same way classical fluids do [2].
A superfluid is characterized by frictionless flow in the absence of viscous effects.For a sufficiently small velocity, the flow around an obstacle is steady laminar flow and no vortices are nucleated [3].Above a critical velocity, the flow around an obstacle creates a drag force which is responsible for the nucleation of quantized vortices [4,5].These vortices signal the breakdown of superfluidity in the system at zero temperature [4,[6][7][8][9][10].Immediately above the critical velocity, pairs of oppositely charged vortices are shed periodically from opposite sides of the obstacle [11].As the velocity of the flow around the obstacle increases, there is a transition from the regular shedding of vortex dipole pairs to an irregular shedding of larger clusters of same-sign vortices, indicating that the system has become turbulent [3,12].The transition to turbulence in superfluid flow past a potential obstacle has been the focus of recent theoretical [2,3,[13][14][15][16] and experimental [11,12,17,18] work.These works have investigated the effect of obstacle shape [14,19,20] and finite temperature effects [15,16] on the critical velocity for vortex nucleation past a single obstacle.
As no real system is truly free of imperfections, disorder is an important consideration in interacting Bose systems, with the interplay between disorder and particle-particle interactions providing a rich test bed for many-body quantum physics.Studies into disorder in BECs have employed impurities [21], rough boundaries [22,23] and optical speckle patterns [24][25][26][27][28][29][30][31][32][33][34]; the latter playing a role in the prediction of a lowered superfluid transition temperature in 2D [28] and 3D [25,26,28], the realisation of Anderson localization [29], and the transition to an exotic Bose glass [30,33].While disorder is an important consideration, few studies have considered the case of a superfluid flow in the presence of a pointlike disordered potential.Such a disorder potential is now experimentally realisable, as new optical techniques employing technologies such as digital micromirror devices (DMDs) allow experiments to have an unprecedented level of control in creating arbitrary shaped potentials [35][36][37].
Forcing a quasi-two-dimensional (2D) superfluid through a disordered potential faster than the critical velocity is a process which injects vortices into the system.These then decay by a process of vortex-antivortex annihilation, which is similar to the coarsening process which takes place after a thermal quench.Such coarsening is a current topic in 2D Bose gases, with investigations into the phase-ordering kinetics of this process being performed in conservative [38][39][40][41][42][43][44] and dissipative [40][41][42][43][44][45] situations, as well as in systems of binary BECs [46], spinor BECs [47,48], and excitonpolariton condensates [49][50][51][52].Previous works on singlecomponent Bose gases have conducted quenches by starting from non-equilibrium initial conditions which tend to rapidly seed an approximately isotropic distribution of vortex dipoles [15,16,44,53,54], while other studies have imprinted a random distribution of vortices with unit charge [43,55] or multiple charges [42].Here, we also observe a system which transitions from a non-equilibrium state containing many vortices towards an eventual equilibrium state, and use an energyand number-conserving description similar to the conservative studies mentioned above.Our system, however, has several key differences.Firstly, the vortex injection in our system is different; unlike the initial conditions discussed above, the vortices which are created by a series of barriers have an anisotropic initial position which depends on the details of the barriers and the flow velocity of the superfluid.Secondly, the vortex injection is not instantaneous; rather vortices are shed over time from the barrier as the barrier moves through the superfluid above the critical velocity.Despite these differences, the system we describe provides a relatively simple way to generate non-equilibrium conditions which can be used to study related coarsening behaviour in a BEC.
In this paper we investigate the dynamics of dense 2D superfluid flow through a point-like disorder potential: a scenario which combines disorder, turbulence and coarsening in a 2D Bose gas.We impose the point-like disorder through an external trapping potential which is taken to be zero everywhere, apart from at a series of points where a localised repulsive barrier is placed.These repulsive barriers, which are Gaussian in shape, may be thought of as a set of blue-detuned laser beams whose intensity can be controlled at any point in space [35][36][37].Unlike the disorder which is imposed by an optical speckle pattern, a key feature of this work is that the barriers which comprise the disorder potential are sufficiently separated (i.e.several healing lengths apart) so that the fluid is homogeneous away from the centre of the barrier.This ensures that it is possible to have a global superfluid phase, since localization of the condensate does not play a role [29], and we can treat quantities such as the speed of sound and the healing length as effectively constant across the system.
The rest of the paper is structured as follows.In Sec.II, we describe the system and its equations of motion.In Sec.III, we calculate the critical velocity for vortex nucleation for different point-like potentials.We begin by placing two identical point-like barriers in a superfluid flow, and study the interplay between relative separation and the incident angle of the barriers on the critical velocity.We then look at a system with many point like barriers, and investigate the link between the density of these point-like barriers and the critical velocity of the system.In Sec.IV, we study the long term behaviour of an initially non-equilibrium superfluid flowing through a disordered potential at varying initial velocities.We measure the condensate fraction, the superfluid fraction, and the superfluid velocity during this process.This illustrates how at short times the superflow breaks down, accompanied by vortex generation and depletion of the condensate fraction.At intermediate times, the momentum of the Bose gas continues to be arrested by interaction with the barriers, as vortex-antivortex annihilation begins.Over longer times vortices continue to annihilate and thermalization takes place; the gas recondenses and superfluidity is restored.In Sec.V, we investigate the effect of varying the effective barrier width on the vortex decay rate.For small point-like barriers (radius on the order of the healing length), the vortex decay rate follows the expectation for a thermal quench.We show that for sufficiently large barriers this changes, and at the same time vortex pinning becomes an important effect in the dynamics of the system.Sec.VI contains our conclusions.

II. SYSTEM AND NUMERICAL IMPLEMENTATION
We consider an obstacle which is moving at a steady velocity v through a superfluid which is otherwise uniform in the xy plane, and trapped strongly enough in the z direction that all excitations are suppressed in this direction.Such a 2D system, when comprised of a weakly interacting atomic Bose gas at finite temperature, can be described by a wavefunction Ψ which obeys the projected Gross-Pitaevskii equation, PGPE, Here, µ 2D is the chemical potential and the strength of the atomic interactions is parameterized by g 2D = √ 8π 2 a s /ml z , where m is the atomic mass, a s is the s-wave scattering length, and l z = /mω z is the harmonic oscillator length in the z direction.We impose a uniform flow with velocity v in the x direction by multiplying the initial wavefunction by a phase gradient (see, for example, Ref. [19]).The crucial feature of the PGPE, beyond the ordinary non-projected Gross-Pitaevskii equation, is the projection operator P which implements an energy cutoff in the basis of non-interacting single particle modes.When working at finite temperature, this allows one to set the cutoff so that modes below the cutoff are highly occupied.In this regime quantum fluctuations are relatively small and the classical field description is accurate [56].
Alternatively, we can consider the system in which the obstacles are dragged through the fluid at some velocity v.In this system, the coordinate of the obstacle reference frame is r = r L + vt, and the lab-frame wavefunction Ψ(r, t) = Ψ L (r L , t).The PGPE governing the lab-frame wavefunction is given by where the Gallilean shift to the obstacle frame (from the labframe) is given by the v • p term, with p = −i ∇ the usual quantum momentum operator [57,58].
To simulate an obstacle which is a collection of point-like barriers, we use the sum of N B repulsive Gaussian potentials, which have their centers at (x k , y k ).These barriers each have an effective cylinder width which can be estimated from the zero density region of the Thomas-Fermi approximation, 2a ln (V 0 /µ 2D ).In contrast to previous works which use hard-walled barriers [3,14], we use soft-walled barriers (with V 0 = eµ 2D ) and where the critical velocity is lower [19].Unless otherwise stated, we take barriers to have a narrow waist, a = ξ, thus providing a point like potential with an effective cylinder width 2ξ.
In what follows, we take v = −v obst x.The single-particle modes for this system are plane waves satisfying |k| < k cut , for some wave-number cutoff k cut .To implement the cutoff, we ignore V obj as it will not affect the potential on the scale of 1/k cut , and hence does not affect our choice of basis functions.The PGPE is evolved numerically, with doubly periodic boundary conditions, using an adaptive Runge-Kutta method (implemented using XMDS2 [59]) on a L x × L y grid with N x × N y grid points.We take the energy cutoff to be k cut = πN x / (2L x ) − π/L x .In the rest of the paper, we typically express quantities with reference to energy µ 2D , healing length, ξ = / √ mµ 2D , density, ρ = µ 2D /g 2D , and the speed of sound, c = µ 2D /m.Consequently, times are expressed in units of τ = /µ 2D .The number of grid points in our simulations is chosen such that there are two computational grid points per healing length.

III. CRITICAL VELOCITY OF POINT-LIKE DISORDERED POTENTIALS
A. Method In order to find the critical velocity, we first find the ground state of the condensate in the presence of the point-like potentials.To do this, we evolve the damped PGPE, found by multiplying the right hand side of Eqn. ( 1) by (1 − iγ), where γ is a phenomenological damping parameter [60], with stationary barriers v obst = 0, and for γ = 1, up to t = 5000τ.This converges to a wavefunction which is approximately the ground state of the system, and which will be the initial condition for all of the following simulations.We then set γ = 0 and evolve Eqn.(2) , whilst smoothly ramping up the velocity [14] according to Smoothly increasing the velocity in this way prevents the generation of sound which would be caused by instantaneously setting v obst = v f .This simulation is run for 1000τ.The value of v f is increased discretely in small increments until vortices are observed to be shed from the potential.For reference, the critical velocity of a single point-like barrier is v crit /c = 0.5625 ± 0.0025.

B. A Pair of Point-like Barriers
We begin by finding the critical velocity of two point-like barriers, as we vary the relative distance and angle between these barriers.Without loss of generality, we place one barrier at the origin, and one barrier at (−R cos α, −R sin α).The results of this are plotted in Fig. 1.
When α, the angle between the barriers in the direction of the flow, is small, the system has an increased critical velocity as the barriers are behind each other in the direction of the flow, becoming streamlined.As α increases, the critical velocity decreases since the barriers become a more like an effective elliptical obstacle, causing a denser wake [14].An important observation that we make is that in the case where R = 4ξ the two barriers act as one larger (essentially elliptical) barrier, and for v v crit will shed only one dipole pair of vortices.In the cases where R ≥ 8ξ, the barriers act independently and both of the point-like potentials will emit a dipole pair, for flow speeds just above v crit .The flattening of the curves indicates that, as we would expect, v crit tends towards the single barrier result as R → ∞.

C. Multiple Barriers
Having found the critical velocity for a pair of point-like barriers, we now find the critical velocity for N B barriers which are placed at random in the cell, subject to a minimum separation of 4ξ.A vortex detection algorithm similar to Ref. [61] is used to automate the search.
In Fig. 2 we plot the critical velocity of a disordered system, as a function of the angle between the nearest neighbour pair of point-like potentials in the particular disorder realiza- tion, n.n.α.We choose this measure because we anticipate that the critical velocity of a particular potential will be most sensitive to the configuration of the pair of barriers with the smallest separation, as shown by the range of values in the blue curve of Fig. 1.The panels of Fig. 2 correspond to the binning of all realizations according to the nearest neighbour distance between the closest two point-like barriers in each re-alization, while the type of marker represents the total number of barriers in the system, N B .The gray shaded area indicates the region which contains the critical velocity of a system of 2 point-like barriers whose separation distance corresponds to the separation distance of the panel.For larger n.n.α, the nearest neighbour interactions of the closest pair of point-like barriers dominate the critical velocity, as can be seen by the points lying within the gray shaded region.Where the closest nearest neighbour barriers form a streamlined barrier, given by smaller n.n.α, the critical velocity is smaller than the two barrier case; this is due to two factors.Firstly as N B increases, so does the probability that other (non-closest) pairs of nearestneighbour barriers are separated by a similar distance but have a large angle against the flow, creating an efficient vortex emitter.Secondly, given that there are multiple barriers in the system, the critical velocity is limited by the single barrier case -any barrier which is sufficiently separated ( 20ξ) from the other barriers will act independently, and cause vortices to be present in the system as soon as the flow velocity is greater than the critical velocity for a single point-like barrier.Indeed, we observe that the critical velocity of a point-like disordered potential is bounded above by the lowest of: (a) the critical velocity of a single barrier; (b) the highest critical velocity of the two barrier test case for equivalent nearest-neighbour separation of the closest two barriers.

IV. ARREST OF A SUPERFLOW: VELOCITY DEPENDENCE A. Overview
Driving a superfluid through a disordered potential faster than the critical velocity injects vortices into the system.The resulting non-equilibrium dynamics are a key object of study in two dimensional quantum turbulence, and have been employed as the initial conditions of studies into quenches both in the highly turbulent clustered case [42], and the dipole dominated case [15,16,44,53,54].In this section, we consider a superfluid which is initially flowing through a disordered point-like potential, with an imposed velocity which is greater than the critical velocity of the potential.We observe that the reaction of the fluid is to be arrested by the barriers, with effectively viscous effects entering the system, before the system equilibrates.The manner in which this disordered system reaches an equilibrium state has connections with quantum turbulence and coarsening in 2D Bose gases.
In this section, we consider one disordered potential with N B = 25 barriers in a system with dimensions L x = 256ξ by L y = 64ξ.As the system consists of a superfluid initially moving through the point-like barriers above the critical velocity, the formation of elementary excitations causes the system to fall out of equilibrium [63,64].To investigate the turbulence in such a system, we measure the condensate and non-condensate fractions, the velocity of the condensate and non-condensate fractions, the superfluid and normal fluid fractions, and the number of vortices which are nucleated by the obstacle.A similar setup has also been considered in a 4v crit , and column (iv) v obst = 1.6v crit .Row (a) displays the condensate (blue circles) and non-condensate (red pluses) fractions.Row (b) is the velocity of the condensate mode (blue circles), the velocity of the the non-condensate mode (red pluses), and the approximate velocity of the normal fluid (black curve) given by Eqn.(11); the grey dashed line indicates v obst , while the grey dotted line indicates zero velocity.Row (c) shows the superfluid fraction computed using the current-current correlations (blue circles), the approximated superfluid fraction described in the text (blue dotted line), and the normal fluid fraction (red pluses).Row (d) plots the vortex number; insets show the vortex number on a log-log scale.The markers are added to help distinguish between curves, rather than indicating individual data points.In the Supplemental Material [62] we provide example movies of these simulations.
1D Bose gas subject to a series of randomly positioned delta scatterers [65], and more recently in a 2D Bose gas flowing through a blue-detuned speckle potential [32,34], however our random potential consists of point-like barriers between which there is a dense superfluid flow.
In order to perform ensemble averaging we add a small amount of complex white noise to the groundstate of the wavefunction (with amplitude approximately equal to 1% of the background density).This small amount of initial noise ensures that the system dynamics, and in particular vortex motion, differs in each realization, such that the statistics are not dominated by particular vortex trajectories.Averaging over this ensemble allows us to reliably calculate condensate fractions, condensate velocities, and superfluid fractions, as will be described in the following subsections.
In order to perform an analysis of this system in the long time limit, we evolve the PGPE prescribed in Eqn.(1).This is expressed in the frame where the barriers are at rest and the wavefunction is given an instantaneous initial boost, where Ψ (g) (r) is the wavefunction in the ground-state of the system, is the amount of noise to be added, v int = v obst /∆v ∆v, where • is the ceiling function, and ∆v = 2πcξ/L x is the smallest velocity representable on the grid in the x direction.The random variables are (r) ∼ U[0, 1] and ϕ(r) ∼ U[0, 2π), and we renormalize such that the initial condition has the same normalisation as the ground state.We choose to evolve Eqn.(1) as it keeps the late-time, closeto-equilibrium momentum distribution of the system close to symmetric about k = 0. Since the k-space cut-off imposed by the projector is symmetric about k = 0, this choice ensures that the system is evolving towards a well-defined PGPE equilibrium, and hence that the calculation of the momentummomentum correlations required to compute the superfluid fraction (Sec.IV D) may be performed without the need to perform a gauge transformation.

B. The Condensate and Non-Condensate Fractions
Where obstacles are dragged through a system at a speed sufficiently above v crit , a large number of vortex-antivortex pairs are nucleated, forming a complicated phase field [12].Annihilation events between the vortex-antivortex pairs lead to the generation of sound in the system, which causes a depletion to the condensate fraction; this marks the onset of a dissipative regime.
Using the criterion of Penrose and Onsager [66], within the c-field formalism [56], we calculate the condensate and noncondensate fractions from the one body density matrix where • T indicates short time averaging.This fraction is calculated for each of the trajectories before averaging over all trajectories.The condensate number can be identified as the largest magnitude eigenvalue of the one-body density matrix, while the corresponding eigenvector, ψ 0 , is the condensate mode.Under this formalism, we deconstruct the wavefunction into contributions from the condensate mode and a non-condensate mode, where n 0 is the condensate fraction and n nc is the noncondensate fraction, with n 0 + n nc = 1.The non-condensate mode, ψ nc , is the sum of the eigenvectors of G 1B excluding ψ 0 .Within the c-field formalism, the condensate and noncondensate modes are orthogonal.The average condensate fractions for obstacles which are dragged through the fluid with velocity v obst ≥ v crit are plotted in Fig. 3, row (a).In systems where v obst ≥ v crit , there is an initial depletion of the condensate fractions as the barrier sheds vortices which are subsequently annihilated, ultimately heating the system [9].We observe that the size of the depletion of the condensate fraction (and therefore the spike in the non-condensate fraction) monotonically increases with the velocity of the obstacle which is consistent with the finding of Refs.[4,19].Eventually, the energy that is injected into the system by annihilation events is distributed among phonons [4,64], and the system relaxes to a uniform flow.This is shown by the increase and then plateauing of the condensate fraction, indicating that the system has reached equilibrium and no further vortices are shed.We see that the long term behaviour of the condensate fraction depends on the speed of the obstacles; in the system where v obst = v crit the final condensate fraction is n 0 ≈ 0.84, where as for v obst = 1.6v crit the final condensate fraction is n 0 ≈ 0.63.This is to be expected, as the faster initial boost injects more energy into the system creating a hotter final state.

C. The Velocity of the Condensate and Non-Condensate Modes
As the shedding of vortices causes the depletion of the condensate fraction, we would expect the presence of thermal effects to lower the critical velocity [15] which in turn would lead to the nucleation of more vortices, until the condensate is depleted.In fact, since the long term behaviour of the condensate fraction is to equilibriate, we deduce that the system stops shedding vortices.This indicates that the system dynamically reacts to the obstacle velocity.
The velocity of the condensate mode ψ 0 (equivalently, the non-condensate mode ψ nc ), is [67] where the index k ∈ {0, nc} and ∇ = ξ∇.We calculate the average velocities of the condensate, v 0 (t x, and plot them in the barrier reference frame in Fig. 3, row (b).
In the presence of the barriers, the fluid nucleates vortices where v obst ≥ v crit .As these vortices nucleate, a phase winding is imparted on the wavefunction, which accelerates the fluid in an attempt to match the speed of the obstacle.In the long time limit, we observe that the velocity of the condensate and non-condensate modes is arrested by the barrier, suppressing further vortex nucleation.The drag force which is exerted on the obstacle potential by the fluid can also be measured [3], and we find that this vanishes as the system evolves.
We note that we also expect to see a variation in the velocity in the y direction, as different configurations of the barriers act as airfoils, [20], causing a lift effect.In our simulations, since v = v obst x, this variation depends on the configuration of the barriers.In any realization it is small in comparison to the velocity change in the x direction.

D. The Superfluid and Normal Fluid Fractions
In order to understand the mechanism by which the velocity of the condensate and non-condensate modes are arrested by the barrier, we calculate the superfluid and normal fluid fractions.We calculate the superfluid fraction in two ways; further details of each approach are given in Appendix A.
Firstly, we assume that the current, J, of the wavefunction can be decomposed into contributions from a superfluid component (which flows without energy loss) and a normal fluid component (which is subject to viscous effects).We expect that the normal fluid will move with the barriers, and so in the frame of reference where the barriers are stationary, the velocity of the normal fluid will vanish in equilibrium.Since the superfluid velocity is locked to the condensate velocity [68,69], assuming zero normal fluid velocity leads to an estimate of the superfluid fraction f s using J = ρ f s v 0 , where v 0 is the condensate velocity introduced in the previous section.
Secondly, we compute the superfluid fraction by noting that the (α, β) element of the current-current correlations of the system in thermal equilibrium can be written as in the limit of vanishing momentum [70], where α, β ∈ {x, y}.
Here f s and f n are the superfluid and normal fluid fractions, and F (J) is the Fourier transform of the current of the wavefunction.The angled brackets • T R indicate that the correlations are found by short-time averaging and by averaging over the ensemble of initial conditions.It is then possible to extract the superfluid and normal fluid fractions by fitting the current-current correlations of the wavefunction to the right hand side of Eqn. ( 9) [61].Formally this method is only valid at equilibrium.Here we employ it with ensemble-and shorttime-averaging to give a dynamic measure.We expect this measure to be quantitatively accurate at late times as equilibrium is approached, since we observe the current-current correlations are well fitted by the expected functional form of Eqn. ( 9) at later times.At earlier times, further from equilibrium, fits to the expected functional form of the correlations fail, and the measure only provides a qualitative indication of the lack of superfluidity.The superfluid fraction of the system is plotted in Fig. 3, row (c).Where the superfluid fraction (computed using the current-current correlations) is negative, or the normal fluid fraction is greater than 1, it is clear that the condition of vanishing momentum is not met.This condition is better fulfilled at later times, where the velocity of the fluid has been arrested by the barrier [see Fig. 3 row (b)].The fluid must respond to the boost which is initially imposed, and so the velocity of the normal fluid at t = 0 is not necessarily zero.This explains why, at very early times, the superfluid fraction computed by decomposing the momentum of the wavefunction is greater than 1.
At earlier times, there is a jump in the normal fluid fraction which equates to the absence of superfluidity.It is this mechanism which causes the fluid to be arrested by the barriers: the appearance of many vortices is associated with a rise in the normal fluid component which is subject to viscous effects causing the fluid to be decelerated by the barrier.At later times, the superfluid fraction grows and equilibrates with the fluid velocity now approximately zero.By the end of the simulation, both measures of the superfluid fraction are close to each other.As the velocity of the barriers is increased, the final superfluid fraction decreases.In the system where v obst = v crit the superfluid fraction, averaged over the last 20% of the simulation is fs = 0.94 (with standard deviation 0.0009) using current-current correlations, and fs = 0.93 (with standard deviation 0.0015) by decomposing the current of the wavefunction; in the system where v obst = 1.6v crit these values are fs = 0.85 (with standard deviation 0.0014) and fs = 0.83 (with standard deviation 0.0047) respectively.This is an analogous result to the depletion of the final condensate fraction as v obst increases, as discussed in Section IV B.
While a slow but non-zero final velocity of the superfluid (i.e., a slow remnant superflow) is not physically unexpected, it is interesting that we do not observe the non-condensate velocity v nc reaching zero over the timescale of our simulations, as can be seen in Fig. 3 row (b).However, as noted above and as can be seen in Fig. 3 row (c), while the superfluid fractions found by current-current correlations and found by assuming that the momentum of the fluid is due entirely to the superflow are close, the two quantities are not equal.There is also a substantial difference between the superfluid and condensate fractions at late times.These observations suggest that part of the non-condensate fraction contributes to the superflow.This prompts a useful consistency check on our results: the expected total momentum of the fluid can be written in terms of either the superfluid and normal or the condensate and noncondensate components.Therefore, for the average velocities in the x-direction we should have where v s and v n are the superfluid and normal fluid velocities.Assuming that the superfluid velocity is locked to the condensate velocity, v s = v 0 , but relaxing the assumption that v n = 0, we can use our estimates of n 0 , n nc , v 0 and v nc from the Penrose-Onsager analysis and our estimates of f s and f n from the current-current correlation analysis to extract the normal fluid velocity from Eqn. (10) as The value of normal fluid velocity v n obtained from Eqn. (11) is shown in Fig. 3 row (b).One can see at late times that the zero value v n = 0 is within the range of the time-variation in this quantity, and the average v n is much closer to zero than the average v nc .We expect the slight remaining offset of the average v n from zero results from the combination of the various statistical uncertainties in our simulations and the Penrose-Onsager and current-current correlation analysis that feed into Eqn.(11).Overall our simulations and analysis show a consistent picture that, over the timescale of our simulations, interaction with the barriers has resulted in a remnant superflow at well below the cricital velocity coexisting with a normal fluid component that has been slowed to very close to zero velocity with respect to the barriers.

E. The Vortex Number
Since the reaction of the fluid is to accelerate to catch up with the barriers, vortex-antivortex pairs are shed from the barrier only at the beginning of the simulation.This leads to a peak in the vortex number as seen in row (d) of Fig. 3.It is evident that the amplitude of the peak in N v increases as v obst increases; this is because the vortex shedding frequency increases with the velocity of the obstacle [19].
At the end of the simulation it is possible that a small number of vortices remain in the system.The average number of such vortices at late times increases as the late-time condensate and superfluid fractions decrease.Typically, for the barriers considered in this section, this small number of vortices are not pinned to barriers but are free to move, and hence consistent with thermal vortices in the fluid.We will discuss the role of free vortices and vortices which become pinned to the barriers in more detail in the next section.
It should be emphasized that, while the results presented in this section are measurements of one disordered potential averaged over an ensemble of ten initial conditions, these results are applicable to other disordered potentials.We have checked that the results presented in Fig. 3 are qualitatively the same for other N B , so long as v crit is the same (within error bars).The effect of simulating a system which has a higher (lower) v crit is simply to steepen (flatten) the curves seen in Fig. 3, while the long-term behaviour is unchanged.

V. ARREST OF A SUPERFLOW: SCALING AND TURBULENCE
A. Overview Until now we have only considered disordered potentials which consist of a number of point-like barriers, with an effective radius of 1ξ, randomly placed in a periodic cell.We now extend our parameter space to consider disordered potentials consisting of barriers with a greater effective radius, and focus on analyzing vortex decay processes.In this section we consider a square domain with dimensions L x = L y = 256ξ.
The numerical simulations which are carried out in this section can be related to practical experiments.Periodic boundary conditions, such as those imposed in our simulations, can be realised in one direction in experiments using ring traps [71] It is possible to impose a persistent superflow current in such a geometry by stirring [72] or optical methods [71], creating a superflow in the periodic direction.Technology such as DMDs could be used to paint the stationary disordered potential in part or all of the ring trap [37].For a large, annular (i.e., tightly confined in the z-direction) ring trap, the main difference from our simulations here would be the lack of periodic boundary conditions perpendicular to the flow.We do not expect that difference to play a crucial role in the dynamics as long as the difference between inner and outer radii of the annulus is a large number of healing lengths.Interestingly, in addition to the studies performed here, in such a system one could switch off the disorder potential after the initial burst of vortex injections; this could be used as a controllable way to inject a vortex distribution and study the resulting coarsening dynamics without the point-like disorder.

B. Vortex Decay Rate
The rate at which a distribution of vortex dipoles in a quasi-2D Bose gas decays has been the subject of much discussion over the last decade [17,39,54,73,74].The vortex decay rate is expected to be connected to the growth of the correlation length of a system, L c .As the system relaxes after a quench, L c should become the only relevant length scale, and it is predicted that L c grows as L c (t) ∼ t 1/z , where z is the dynamical critical exponent [75].It is also predicted that, for randomly distributed vortices in a homogeneous system, the vortex number and the correlation length are linked as N v ∼ L −2 c .Based on experimental observations, the suggested phenomenological rate equation for N v is [17] Single vortex annihilations are prohibited as vortices are topologically protected, meaning that Γ 1 N v describes the drifting of vortices out of the condensate at boundaries (a onevortex mechanism), while Γ 2 N 2 v represents the rate of vortexantivortex annihilations (a two-vortex mechanism, in this model).However, the decay rate given by Eqn.(12) does not match with the results of zero-temperature GPE simulations [39,54,73,74].This has led to the proposal of a corrected idealized decay rate [74] where it is argued that the drift and annihilation processes have a N 3/2 v and N 4 v dependence respectively.It has since been shown [54,73,74] that for a homogeneous system at zero temperature N v ∼ t −1/3 which is indicative of a four-vortex process, while the addition of dissipation (finite temperature effects) or trapping potentials removes the need for a fourth vortex [55] (the N 4 v scaling which describes a four-vortex annihilation process was also observed numerically in Ref. [39]).
Due to the large proportion of our simulations which occur after the peak in vortex number, it is possible to study the long-time behaviour of vortex decay in our disordered potential systems in a similar fashion.As discussed earlier, we use N v → N v − 1 → + 1   the plaquette technique [61] to enable vortex detection.Unlike before, where we focused on barriers with effective radii 1ξ, for barriers which have an effective radii 2ξ there is a significant zero density region where the phase of the condensate is ill-defined.Naively applying the plaquette technique here leads to the detection of spurious vortices.However, it is also possible for a net number of quanta of circulation to genuinely be present at this low density region: we define this number of quanta as the winding number of the barrier W k (for the kth barrier).The winding number can also be interpreted as a number of pinned vortices.Hence, when computing the vortex number we detect both the number of mobile vortices N v , using the plaquette technique and excluding the density-depleted regions, and the total number of pinned vortices which is computed using a loop integral technique described in the next section.The evolution of the vortex numbers for a system with N B = 25 barriers of varying effective radii is shown in Fig. 4 (a)-(d).In Fig. 4(a) and (c) we plot only the number of mobile vortices N v .In Fig. 4(b) and (d) we plot the total number of vortices (mobile and pinned), N v + W. For the narrowest barriers we consider, the vortex decay rate appears to follow a t −1.1 power law for effective barrier radii of ξ, and a t −1.2 power law for effective barrier radii of 3ξ/2, as can be seen in panel (a).In a system where the vortex number only decays via vortex-antivortex annihilations, Eqn.(12) predicts that N v ∝ t −1 .The fact that the observed power laws are relatively close to t −1 for the narrowest barriers is indicative of the fact that vortex decay is a two-vortex process in this system.For barriers which are larger than the typical size of a vortex core (i.e., have an effective width which is greater than a few healing lengths), the vortex number appears to decay exponentially, as can be seen in panel (c).This is consistent with a solution to Eqn. (12) where a one-vortex mechanism is dominant, i.e., N v ∝ exp (−Γ 1 t).This suggests that for wider barriers, at late times in the simulation, vortices are colliding with a barrier more often than they are colliding and annihilating with a vortex of the opposite sign.We discuss the effects of vortices colliding with barriers in the following section.

C. Pinning to Barriers
As well as measuring the rate at which the number of vortices decay, we have also measured the number of vortices which become pinned to the barriers.The pinning and unpinning of superfluid vortices is an important physical process for understanding the mechanism of neutron star glitches [76][77][78][79][80], and is also of interest in systems with macroscopic container defects [81][82][83], as well as spin-down experiments with helium [84,85], and laboratory BECs [86].The microscopic process by which a vortex becomes pinned to a density depleted region has recently been studied by Ref. [87].For systems where impurities exist, it is energetically favourable for a vortex to be contained within the zero-density region, as there is no cost in energy to create a vortex core [88].
As described above, we define pinned vortices in terms of the net quanta of circulation around a barrier, which is well defined as the branch cut representing a discontinuity in the phase extends into the non-zero density region of the condensate (i.e., it is not a spurious vortex caused by the phase not being well defined in the zero density region at the centre of the barrier).For each barrier in a given potential, we can measure the winding number W k by integrating around a loop containing the barrier (see Appendix B for details of the numerical method).Examples of the phase of a barrier with no pinned vortices, one pinned vortex, and two pinned vortices are shown in Fig. 5 panels (i), (j) and (k).Also shown is the approximate location of the radius of the circular "exclusion zone" which we choose when counting the number of mobile vortices.A slightly larger circular loop is used to measure the winding number.It should be noted that in any one trajectory the time-dependent values of the numbers of mobile and pinned vortices may display fluctuations in time that depend on the precise choice of radii for these circles, especially when two or more barriers are close together.While we were unable to find choices that eliminate these fluctuations in any one trajectory, we find the averaged results are relatively insensitive to the choice of radii.
At early times, the system is in a highly non-equilibrium state, and many vortices are periodically shed by the barriers.However, by t In the Supplemental Material [62] we provide example movies of these simulations.
each barrier is steady.This can be seen in Fig. 4 panel (e).
For larger barriers, the number of mobile vortices in the system decays as N v ∝ exp (−Γ 1 t), suggesting that the vortices are annihilating with the barriers.From our observations of the simulations, we suggest that there are 3 processes taking place here.Process I: a vortex collides with a barrier which has a number of vortices with the same sign pinned to it.Here the number of mobile vortices decays, N v → N v − 1, while the number of pinned vortices increases, W → W + 1. Process II: a vortex collides with a barrier which has a number of vortices with the opposite sign pinned to it.Here the number of mobile vortices decays, N v → N v − 1, but the number of pinned vortices also decays since the mobile vortex annihilates with one of the pinned vortices, W → W − 1. Process III: a dipole pair collides with a barrier which has a number of vortices pinned to it.Here, the number of mobile vortices decreases by two, N v → N v − 2, however the number of pinned vortices remains the same, since one of the dipole pair will annihilate with the vortices of opposite sign in the barrier, while the other vortex in the dipole pair will remain and will become pinned to the barrier, W → W. As each barrier sheds an equal number of vortices and anti-vortices, Processes I and II take place with approximately the same frequency, conserving the pinning number W. Process III, which also conserves the winding number, happens far less frequently.However, this process may perhaps explain the slight modifications to the exponential decay which we see in Fig. 4. We assume that collisions between three or more vortices and a barrier are so rare as to be negligible.The probability of observing a given winding number can be seen in the histograms in Fig. 5, where the data is taken from 10 4 τ ≤ t ≤ 2 × 10 4 τ.As we can see, for narrow barriers vortex pinning is not an important feature.However, for barriers which are significantly larger than a vortex core, a significant number of the barriers do have a vortex or anti-vortex pinned to them (W k = ±1), and the largest barriers which we consider support the pinning of multiple vortices (|W k | > 1).Examples of this behaviour can be seen in Fig. 6.
It can be seen in Fig. 4 that the rate at which the number of mobile vortices decays becomes quicker as the effective radius is increased past 2ξ, and is at its fastest for barriers which have an effective radius of ≈ 5ξ.This may be attributed to the fact that for barriers with an effective radius greater than 2ξ we have observed that is is more likely for a barrier to support the pinning of vortices; this provides a mechanism to lose mobile vortices via Process I, above.For barriers which have a larger effective radius than 5ξ, we have observed that it is possible to have multiple vortices pinned to a barrier.Multiple pinning creates a stronger velocity field around the barrier than single pinning does; this could explain why the rate at which the number of mobile vortices decays slows slightly as the effective barrier radius increases above 5ξ.

VI. CONCLUSION
In this paper we have studied the effect of dragging a disordered point-like potential through a superfluid which is initially in the ground state.We have seen how the critical velocity of two point like barriers depends on the relative distance and angle between the barriers.We have then determined the critical velocity for a system which has up to 50 point like barriers at randomized locations, and shown that the critical velocity of such a system can be mapped on to the two-barrier case by considering the separation and angle with respect to the flow of the closest nearest-neighbour pair of barriers in the disorder potential.
Using PGPE simulations, we investigated the evolution of a system in which an initial superflow, moving at or above the critical velocity, is disturbed by a stationary point-like disorder potential.This strongly non-equilibrium initial condition causes the nucleation of vortices and depletion of the condensate and superfluid fractions.We observe that the reaction of the fluid is to accelerate to a final velocity closer to the obstacle velocity.This suppresses the nucleation of further vortices, and the fluid re-condenses and some superfluidity is restored.
We extended our parameter space to consider the effect of larger barriers in the system, and investigated the way in which this affects the decay of the number of vortices in the system.It is clear that the presence of randomly placed barri-ers that have an effective width which is larger than the characteristic size of a vortex core, modifies the form of the vortex number decay from the behaviour identified in previous theoretical works without a disordered potential.Within the limits of our numerical analysis, it appears as though the vortex decay rate no longer follows a t −1 power-law scaling which is indicative of vortex-antivortex annihilations, but rather the vortices collide with the barriers which make up the potential, causing an exponential decay.This one-vortex decay process is confirmed with our observations of the simulations.Finally, we observe that for these larger barriers vortex pinning becomes a relevant phenomenon, with the largest barriers which we consider supporting the pinning of multiple vortices.
With an appropriate trapping geometry, it may be possible to experimentally study a system equivalent to the one studied here in which the disordered potential arrests a superflow.the unwrapped phase difference between neighbouring points, ∆θ j, j+1 = θ α j − θ α j+1 .(B2) It is necessary to unwrap the phase in this way to ensure that the phase is continuous between neighbouring points [61], however working on a discrete grid this continuity is poorly defined as there may be jumps in the phase of 2π; to correct for this we add multiples of 2π so that |∆θ j, j+1 | < π .The winding number is then computed as

FIG. 1 .
FIG.1.The critical velocity of two point-like barriers, with separation distance R, and angle α incident to the direction of the flow.Blue circles represent barriers with separation R = 4ξ, red squares represent barriers with separation R = 8ξ, green crosses represent barriers with separation R = 12ξ, and orange pluses represent barriers with separation R = 16ξ.The olive region is the critical velocity of a single point-like barrier, plotted as a guide to the eye (the width indicates numerical uncertainty).The error from the systematic uncertainty of increasing v x in discrete steps is smaller than the symbols used.The inset shows a schematic of the experimental set up of the two point-like barriers.

FIG. 2 .
FIG. 2. The critical velocity of a disordered potential with N B point-like barriers.Organised by nearest neighbour distance (n.n.R) between the barriers, panel (a) has 4ξ ≤ n.n.R < 5ξ, (b) has 5ξ ≤ n.n.R < 6ξ, (c) has 6ξ ≤ n.n.R < 7ξ, (d) has 7ξ ≤ n.n.R < 8ξ, (e) has 8ξ ≤ n.n.R < 12ξ, (f) has 12ξ ≤ n.n.R. Different markers represent varying barrier density.The gray shaded area indicates the region containing the critical velocity of an isolated pair of point-like barriers when their separation distance lies within the range of nearest neighbour distances for the panel.The olive shaded area is the critical velocity of one point-like barrier (the width indicates numerical uncertainty).

FIG. 4 .
FIG.4.The decay of the number of vortices in a system as the barrier width varies.Panels (a) and (c), the decay of the number of mobile vortices, N v .Panels (b) and (d), the decay of the total number of vortices, (N v + W).Panel (e), the decay of the number of pinned vortices, W. Panels (a) and (b) are plotted on a log-log scale, while panels (c)-(e) are plotted on a semi-log scale.The power law N v ∝ t −1 , black dashed line, and the exponential decay N v ∝ exp (−Γ 1 t), black dotted line, are added as guides to the eye.The markers are added to help distinguish between curves, rather than indicating individual data points.

FIG. 5 .
FIG. 5. Normalized histogram of the winding number, W k , for barriers with effective radius (a) 1ξ, (b) 3ξ/2, (c) 2ξ, (d) 3ξ, (e) 4ξ, (f) 5ξ, (g) 6ξ and (h) 7ξ.Right: examples the phase of the wavefunction, Arg (Ψ), around a barrier with (i) no pinning, W k = 0, (j) one pinned vortex, W k = +1, and (k) two pinned vortices, W k = +2; the black circle is approximately the boundary of the zero-density region of the barrier.A slightly larger circular loop is used to compute the integral required to calculate W k .

FIG. 6 .
FIG. 6.A still at time t = 14460τ from the simulations of barriers with an effective radius of 3ξ, (a) and (c), and simulations of barriers with an effective radius of 7ξ, (b) and (d).In panels (a) and (b), the density of the wavefunction is shown, while different markers indicate the winding number W k of a barrier, and the position of a vortex or anti-vortex.The phase of the wavefunction is shown in panels (c) and (d).In the Supplemental Material[62] we provide example movies of these simulations.