Critical velocity for quantized vortex formation in a super�uid with a plate-shaped obstacle

We theoretically investigate the dependence of the critical velocity for quantum vortex formulation on the size L of a plate-shaped obstacle moving in a uniform Bose-Einstein condenste in quasi-two-dimensions. Numerical simulations of the Gross-Pitaevskii equation reveal that the critical velocity v c is proportional to L − 1 / 2 when L is much larger than the healing length. The power law behavior is quantitatively explained by the potential ﬂow theory of classical hydrodynamics for the large L limit.


Introduction
The analysis of wakes, formed behind an obstacle moving through a fluid, is fundamental to understanding flow dynamics in classical hydrodynamics [1].Wakes are relevant to important aspects of flow phenomena, including the transition between laminar and turbulent flows, often characterized by the Reynolds number.The wakes exhibit various vortex dynamics associated with the generation of twin vortex and von Kármán vortex streets.
Quantum hydrodynamic instability plays a crucial role in the study of superfluid turbulence and related superfluid phenomena [2].A wake in the superfluid has been studied in some systems, especially in a cold atomic-gas Bose-Einstein condensate (BEC) [3][4][5][6][7][8].There have been various theoretical studies of superfluid wakes using the Gross-Pitaevskii (GP) model [9][10][11][12][13][14][15][16].The characteristics of superfluid wake are (i) the critical velocity above which the vortex nucleation occurs, and (ii) generated vortices are quantized.Motivated by these features, studies of the wake in classical fluid systems have been extended to superfluid systems, such as quantum von Kármán vortex street [7,17] and superfluid Reynolds number [17].In weakly interacting Bose gases, Landau's argument on superfluidity predicts that the critical velocity in a uniform system is given by the sound velocity c s .Additionally, according to the incompressible flow in classical fluids around a cylinder, the critical velocity is expected to be 0.5c s because the local velocity on the lateral sides of the cylinder is as twice as the bulk velocity.The calculations based on the GP model show that the critical velocity is 0.37c s when the size of the cylinder is sufficiently large relative to the healing length [18,19].In cold atomic Bose-condensed systems, localized repulsive potentials made by laser beams have been used in superfluid wake experiments.Such obstacle potentials are often modeled by the Gaussian form in the theoretical or numerical analysis.However, the critical velocity for vortex formation depends on the quantitative details of the obstacle potential [5,13].It makes us difficult to predict analytically universal understanding of the superfluid wake and its dynamics.
In this study, we discuss the size dependence of the critical velocity in a uniform two-dimensional BEC by using a plate-shaped potential as an fundamental case toward a systematic understanding of superfluid wake.The size dependence can be studied more clearly by the plate-shaped potential because its relevant scale is given only by the size of the plate if the thickness of the plate is much smaller than the healing length.By varying the plate width L over a wide range, the relationship v c ∝ L −1/2 between the critical velocity v c and the plate size L is obtained when the width is much larger than the healing length.When the plate moves with a velocity smaller than v c , the velocity field around it can be described by the complex potential flow in fluid mechanics, which suports the asymptotic power-law behavior of the critical velocity.

Formulation
In this section, we introduce a formulation to consider the problem of the wake in a uniform superfluid system using the GP equation.We consider a wake caused by a plate-shaped obstacle moving with the constant velocity u = (0, v) in the y direction in a BEC, which is uniform in the xy-plane and is strongly trapped in the z-direction.We consider a reference frame comoving with the obstacle.The dynamics of the wave function ψ(r, t) in the reference frame can be described by Here, µ is the chemical potential, g in the nonlinear term is the coupling constant in the 2D system, and V p is the obstacle potential given by Below, length, time, and energy are scaled with ξ = ℏ √ 2mµ , ω −1 = ℏ/µ and µ, respectively.

Numerical Methods
In this work, we discuss the critical velocity for vortex nucleation by numerically calculating Eq.( 1) using the imaginary time propagation method by replacing t by iτ .The generation of vortices through the imaginary time evolution means that the energy of the system is decreased by this dynamics induced by the energetic instability.Since this method is associated with the steepest decent of the GP energy functional, it can determine the velocity at which the energy barrier between the vortex-free and non-vortex-free states.In the real-time dynamics, the energy difference between the two states can be consumed by exciting collective fluctuations, i.e., emitting sound waves.If fluctuations exist in advance and are large enough to make the transition over the energy barrier, the vortex nucleation can occur below the critical velocity that obtained in the imaginary time propagation.As such fluctuations do not exist in our numerical simulations, we will obtain a theoretical upper limit of the critical velocity.
The sizes of the space and time grids are ∆x = ∆y = 0.5ξ and ∆t = 0.005ω −1 , and periodic boundary conditions are imposed for both the x and y directions.The potential height of the plate-shaped obstacle potential is V 0 = 100µ.Here, we determined the critical velocity by observing a sudden sharp decrease in energy caused by vortex generation in the following steps.First, we calculate the steady-state solution with v = v 0 < v c .Next, we increase the velocity as v i = v 0 + ∆v from the steady-state solution and check whether the system becomes steady or not.If the steady-state solution is stationary, we increase the velocity again as v i+1 = v i + ∆v and check the stability of the solution.This process is repeated until the steady state becomes unstable (at i = n), and consequently the critical velocity is evaluated as v c = v n − 0.5∆v.Here, we take ∆v = 0.002(ξω), regarded as the error of v c in our numerical analysis.Figure 1 shows the imaginary time development of the density |ψ| 2 and the phase θ = arg(ψ) when the velocity is higher than the critical velocity.We can see that quantum vortices are emitted in the lengthwise direction of the plate.These vortices disappear at the numerical boundary and subsequently the vortex emission repeats.These processes are relevant to phase slippage that leads to the decay of the superflow.In particular, when the plate size is large enough, the relationship v c /c s ∝ L −1/2 is confirmed.

Comparison with the potential flow theory
To explain the numerical results, we consider the potential flow theory in fluid mechanics around the plate-shaped obstacle.We expect that superfluid flow, with a velocity slower than the critical velocity, can mimic a flow of the perfect fluid, an incompressible fluid without viscosity in two dimensions.For our case, a flow around a plate-shaped obstacle of width L can be described by the Joukowski transformation of that around a cylinder of radius L/4 [1].The Joukowski transform has been successfully introduced to describe the elliptical superfluid velocity field around a quantum vortex in the polar phase of spin-1 BECs [20], where the vortex is considered the Joukowski transform of a conventional vortex with the circular velocity field.
First, we consider the velocity potential around a plate under a background velocity (0, v).The complex velocity potential W in the complex plane z = x + iy is given by ( The double sign takes a plus for y ≤ 0 and a minus for y < 0. By differentiating the complex velocity potential W in the coordinates of the complex plane z, we obtain 0.1 1 10 100 Fig. 2 The critical velocity as a function of the plate size L for V 0 = 100µ.The dots represent the numerical result obtained by the method in Sec.3.The critical velocity is scaled by the sound velocity cs = √ 2ξω.The solid blue line is the critical velocity obtained from the potential flow described in Sec.4.Here, the parameter l in Eq.( 7) is chosen as l = 0.531ξ.The two results agree well with each other as the obstacle size L increases.The dash red line is the asymptotic form [Eq.( 8)] of Eq.( 7) for a large plate size.

the complex velocity
Here, v x and v y are the velocities in the x and y directions.Since this complex velocity is based on the flow around a static obstacle, it should be converted to the reference frame moving with the plate as done in the simulation in Sec.4.This is done by adding −iv to the complex velocity w as In Fig. 3, we show the comparison between the velocity field given by Eq.( 5) and that by our numerical simulation, where we set the velocity just below v c .From Figs. 3(a where the velocity potential for superfluids is given by the phase of ψ.The lines in Fig. 3(a) indicate that the complex velocity potential reproduces the streamline around the plate obstacle.To see the consistency more quantitatively, we plot in Fig. 3(c) the magnitude of the velocity along the lengthwise direction of the plate from its edge.In Fig. 3(d) shows the ratio of the velocity obtained by the numerical calculation to the one obtained by the potential flow.Here we compare the local velocities obtained from the potential flow and the GP model for L = 15ξ (Green),50ξ(Black),100ξ (Red).The velocity distribution are close agreement with each other at a sufficient distance from the edges, but is misaligned near the plate, probably due to the fact the quantum pressure and/or compressibility is not included in the potential flow theory.

Critical velocity for macroscopic plates
We analytically estimate the critical velocity by the local velocity around the edge of the plate.It is supposed that the vortex nucleation could occur when the local velocity at the position with a distance l(x = −L/2) in the lengthwise direction of the plate exceeds the sound velocity.To get reasonable results the distance l would have to be on the order of ξ, the thickness of a vortex.Substituting z = L/2 + l in Eq.( 4), the local velocity v loc at a distance l away from the edge of the plate is given by The value of w(z = L/2 + l) is purely imaginary, since the flow only has a component in the direction perpendicular to the plate.Thus, the critical velocity v c is determined by the relation v loc = c s , being written as When the plate size L is large enough, Eq.( 7) is asymptotic to The blue solid line in Fig. 2 is the critical velocity obtained by Eq.( 7) with a fitting parameter l = 0.531ξ, where it is chosen so as to reproduce the behavior at a large L. The asymptotic power law of Eq.( 8) is plotted together by the red dash line for comparison.Our fitting satisfies the condition l/ξ ≲ O(1), as we expected before.We see that the potential theory deviates more from the numerical result for smaller values of L/ξ; the potential flow theory is applicable to a flow around a plate with a macroscopic size L ≫ ξ while a different property unique to the superfluid occurs on the microscopic scale ξ.

Conclusion
In this study, we numerically explored the critical velocity for a moving plate in a uniform superfluid and discussed its dependence on the plate size.When the size L of the plate is large, the critical velocity is proportional to ∝ L −1/2 , which is consistent with the prediction based on the potential flow theory.We thus consider the superfluid wake to be a classical-like wake for large plate-shaped potential.For future research, we intend to clarify the critical velocity for a small plate, which is expected to be different from the classical regime.

0 1 𝑣Fig. 1
Fig. 1 Snapshots of the density profie (top) and the phase profile (bottom) of the condensate wave function by the imaginary time development for L = 40ξ and v = 0.32(ξω).The spatial region of the plot is −150ξ < x, y < 150ξ.The vortices are emitted from the edge of the plate in the direction of the plate, and these vortices disappear at the boundary.

Fig. 3
Fig. 3 Comparison between the result from the potential flow theory and that from the GP model.The profile of the velocity potential Re[W ] obtained by the potential flow theory is shown in the panel (a), while the profile of the phase θ = arg(ψ) of the condensate wave function is shown in the panel (b) for L = 50ξ and v/cs = 0.142.In the panel (a), the several lines overlapped with the potential profile are streamlines Im[W ] of the potential flow.The two velocity potentials have similar distributions, which suggests that the complex velocity potential can describe the profile of the superfluid.Panel (c) shows the distributions of the velocity field from the edge of the plate along the x-direction for v/cs = vc/cs − 0.056.Solid and dashed curves represent the magnitudes of the velocity fields obtained by the numerical calculation and the analytical potential flow for L = 15ξ (Green), 50ξ (Black), 100ξ (Red).The two velocity fields are misaligned in the vicinity of the plate due to the variation of the wave function at the extent of the healing length but they are coincident far from the plate.Panel (d) shows the variance rate of the velocity obtained by numerical calculation to the one obtained by potential flow.We can see that the velocity field in the superfluid as the plate size increases approaches the velocity field obtained by the potential flow.