3.3.1. Superfluid density

The number density, $n \equiv |\psi |^2$ , when combined with (3.2) yields

(3.8) \begin{equation} n |_{y \to 0^+}=0,\qquad \left . \frac {\partial n }{\partial y} \right |_{y\to 0^+} = 0, \end{equation}

which do not arise in a classical fluid. The first of (3.8) shows that no superfluid exists at the solid wall, $y = 0$ , i.e. it constitutes a singular line bounding the flow domain.

Next, because the superfluid is unbounded in the $y$ -direction, and we are considering perturbations local to the solid wall at $y=0$ , it follows that

(3.9) \begin{equation} n |_{y\to \infty } = 1, \end{equation}

which is the bulk density.

Equations (3.8), and (3.9) are the required boundary conditions for the superfluid density, $n$ . Three boundary conditions are needed due to the existence of third-order derivatives in (3.5b ), which differs from a classical Euler fluid where only first-order derivatives appear.

3.3.2. Superfluid velocity

The phase, $S$ , in (3.4) can be expressed as $\tan S = \psi _{\textit{im}}/\psi _{\textit{re}}$ , where the subscripts ‘ $re$ ’ and ‘ $im$ ’ refer to the real and imaginary parts, respectively. It then follows that the velocity field, $\boldsymbol{u} = u \hat {\boldsymbol{x}} + v \hat {\boldsymbol{y}}$ , where $\hat {\boldsymbol{x}}$ and $\hat {\boldsymbol{y}}$ are the Cartesian basis vectors in the $x$ and $y$ -directions, respectively, is given by

(3.10a) \begin{align} u &= \frac {1}{n}\left (\psi _{\textit{re}}\frac {\partial \psi _{\textit{im}}}{\partial x} - \psi _{\textit{im}}\frac {\partial \psi _{\textit{re}}}{\partial x} \right ) = \text{Im} \left (\frac {\frac {\partial \psi }{\partial x}}{\psi } \right )\! , \end{align}

(3.10b) \begin{align} v &= \frac {1}{n}\left (\psi _{\textit{re}}\frac {\partial \psi _{\textit{im}}}{\partial y} - \psi _{\textit{im}}\frac {\partial \psi _{\textit{re}}}{\partial y} \right ) = \text{Im} \left (\frac {\frac {\partial \psi }{\partial y}}{\psi } \right )\!, \end{align}

where $n = \psi _{\textit{re}}^2 + \psi _{\textit{im}}^2$ and $\text{Im} (\ldots )$ denotes the imaginary part.

The extension of the wavefunction, $\psi$ , into $y \lt 0$ satisfies odd parity with the number density, $n \equiv | \psi |^2$ , evidently being of even parity, see Appendix A. Because the GPE, (3.1), is an odd function of $\psi$ , it follows from (3.1) and (3.2) that the general asymptotic form for $\psi$ in the neighbourhood of the solid wall is

(3.11) \begin{equation} \psi (x, y, t) = \sum _{m=0}^\infty y^{2 m +1} f_{2 m +1}(x,t) , \end{equation}

where $f_{2 m +1}(x,t)$ are functions that can be determined (though not required here). Substituting (3.11) into (3.10b ) gives the required boundary condition for the velocity component of the fluid that is normal to the solid wall,

(3.12) \begin{equation} v |_{y \to 0^+} = 0. \end{equation}

No such constraint exists for the velocity component, $u$ , that is parallel to the solid wall. The final boundary condition for the velocity field is that it matches its imposed value far from the solid wall, i.e. as $y \to \infty$ . This together with (3.12) mimics the boundary conditions for a perfect Euler fluid. We note that the derivation of (3.12) from (3.10b ) involves a removable singularity, which arises from the first equation in (3.8).

4.3.1. Boundary conditions for the linearised equations

Equation (4.8) requires four boundary conditions for its solution. Three boundary conditions are reported in § 3.3.1 while the fourth boundary condition arises from the even parity of $n$ (see Appendix A), i.e.

(4.11) \begin{equation} n_1 |_{y \to 0^+}=0,\qquad \left . \frac {\partial n_1 }{\partial y} \right |_{y\to 0^+} = 0, \qquad \left . \frac {\partial ^3 n_1 }{\partial y^3} \right |_{y\to 0^+} = 0, \qquad n_1 |_{y \to \infty }=0. \end{equation}

Equations (4.8) and (4.11) constitute a singular eigenvalue problem for $n_1$ , whose solution is reported in the following section; the corresponding solution to $S_1$ then follows from (4.7b ).

4.4.1. Convergence study

A detailed convergence study is performed for a nominal wavenumber of $k=1$ , to access the sensitivity of the numerical method to its discretisation parameters. Table 1 gives values of the minimum values of $c^2 \equiv c_{\textit{min}}^2$ , calculated according to (4.9), i.e.

(4.14) \begin{equation} c_{\textit{min}}^2 \equiv \frac {\lambda _{\textit{min}}}{4 k^2}, \end{equation}

for $y_M= 25$ and $50$ , where five-figure convergence is evident. It can also be seen that results for the same grid size, $\Delta y$ , agree to six figures, e.g. see $(y_{\textit{max}}, M) = (25,800)$ and $(y_{\textit{max}}, M) = (50,1600)$ . The eigenfunction, $n_1$ , corresponding to the eigenvalue in table 1, is shown in figure 2(a), $n_1$ is real valued and normalised to unit maximum. Results are given for all parameter sets of table 1 and cannot be distinguished to plotting accuracy.

Table 1.

Convergence of the rescaled eigenvalue, $c_{\textit{min}}^2 = \lambda _{\textit{min}}/(4 k^2)$ , with the algebraically minimum real part for $k=1$ , computed as a function of $y_{\textit{max}}$ and $M$ .

Figure 2.

Eigenfunction for $k=1$ , corresponding to the eigenvalue reported in table 1. (a) All computed eigenfunctions for all discretisation parameters in table 1; there are no discernible differences. Log–linear coordinates are used. (b) Comparison of the 10th-order polynomial approximation (solid line), constructed using a special stencil near $y=0$ , with the FD numerical solution (dots). Results here for $y_{\textit{max}}=50$ and $M=3200$ , and log–log coordinates are used.

Figure 3.

As for figure 2(b), but for (a) $n^{\prime}_1$ , (b) $n^{\prime\prime}_1$ , (c) $n^{\prime\prime\prime}_1$ , (d) $n^{\prime\prime\prime\prime}_1$ , where the dashes refer to derivatives with respect to $y$ .

Figure 2(b) gives a comparison of the polynomial approximation with the eigenfunction, $n_1$ , obtained from the first five points (see point (iii) of the numerical method above), with the FD numerical solution; for $y_{\textit{max}} =50$ and $M=3200$ . Corresponding results for the first, second, third and fourth derivatives of this eigenfunction are shown in figure 3. Excellent agreement is again observed for smaller values of $y$ , with divergence occurring for larger values as may be expected.

A further diagnostic of interest is the comparison of the FD numerical solution with the large- $y$ asymptotic solution reported in (B3), which is given in figure 4 for $y_{\textit{max}}=50$ and $M=3200$ . The constant $A \approx 1.3406$ is determined by fitting to the FD numerical solution at $y \approx 12.03125$ . The discrepancy observed for $y \gtrsim 40$ is due to the choice of $y_{\textit{max}}=50$ where the boundary condition, $n_1(y_{\textit{max}}) =0$ , is enforced in the numerical method.

Figure 4.

Comparison of the asymptotic solution in (B3) (dashed line) with the FD numerical solution (solid line). The constant, $A$ , in (B3) is obtained by equating (B3) with the FD numerical solution at $y \approx 12.03125$ . Results for $k=1$ , $y_{\textit{max}}=50$ and $M=3200$ .

4.4.2. Dispersion relation for $0.2 \leqslant k \leqslant 30$

The $c_{\textit{min}}^2$ – $k$ dispersion relation arising from the eigenvalues, $\lambda _{\textit{min}}$ , computed using the present numerical method for $0.2 \leqslant k \leqslant 30$ , is given in figure 5(a). This shows that $c_{\textit{min}}^2\gt 1$ , i.e. $\text{Im} [c_{\textit{min}}] = 0$ , establishing that (a) the wave speed exceeds the sound speed, (b) the waves are marginally stable, as per the statement following (4.6) and (c) move with supersonic speed with respect to the stationary fluid far from the solid wall, i.e. $y\to \infty$ .

It is established in Appendix B, that the required (far field) fourth boundary condition in (4.11) is satisfied when

(4.15) \begin{equation} -\frac {1}{k^2} \leqslant c_{\textit{min}}^2 \lt 1+ \frac {k^2}{4}. \end{equation}

While this fourth boundary condition is enforced in the numerical method, it remains to be shown that (4.15) applies, i.e. that the numerical solution is valid. The lower limit in (4.15) obviously holds because $c_{\textit{min}}^2\gt 1$ , but validity of the upper limit is not immediately evident. To assess its validity, figure 5(b) gives computed data for the function

(4.16) \begin{equation} \nu \equiv 1+\frac {k^2}{4}-c_{\textit{min}}^2, \end{equation}

as a function of $k$ . This shows that $\nu \gt 0$ for $0.2 \leqslant k \leqslant 30$ , and hence that (4.15) holds. The combined data in figure 5(a,b) also show that $\nu \ll c^2$ , i.e.

(4.17) \begin{equation} c_{\textit{min}}^2 \approx 1+ \frac {k^2}{4}, \end{equation}

an approximation that improves as $k \to 0$ ; it appears to be exact as $k \to 0$ . This latter finding together with the large- $y$ asymptotic solution in (B3) establishes that $n_1$ decays at a slower rate with respect to $y$ , as $k$ decreases. That is, $y_{\textit{max}}$ must be increased as $k$ is reduced, which limits our ability to probe the very small- $k$ regime, i.e. $k \lt 0.2$ , due to available computational resources.

Figure 5.

Computed $c_{\textit{min}}^2$ – $k$ dispersion relation arising from the linear stability analysis. (a) Dispersion relation $c_{\textit{min}}^2$ versus $k$ . (b) Difference function, $\nu$ , in (4.16), showing the degree to which the upper limit of (4.15) is satisfied; $\nu \gt 0$ is needed.

4.4.3. Corresponding eigenfunctions, $n_1$ and $S_1$

The eigenfunctions, $n_1$ , corresponding to $c_{\textit{min}}^2$ for several wavenumbers, $k$ , are reported in figure 6(a). These exhibit a unimodal dependence on $y$ , with maxima occurring at $y\approx 1.5$ . The decay rates of these eigenfunctions with respect to $y$ , decreases as $k$ is reduced, as per figure 5(b). There is little variation in the eigenfunctions for $k \gtrsim 5$ . The eigenfunctions, $S_1$ , arising from (4.7b ) are purely imaginary, and

(4.18) \begin{equation} \tilde {S}_1 \equiv \text{Im} [S_1], \end{equation}

is given in figure 6(b). It is seen that $S_1$ is non-zero $y=0$ , establishing finite velocity slip tangent to the solid wall, i.e. $u |_{y \to 0^+} \neq 0$ , as per (3.6) and (4.5b ).

4.4.4. Summary of the linearised solution

The linearised solution to wavefunction, i.e. correct to $O(\epsilon )$ , is

(4.19) \begin{align} \psi = \Big ( \tanh (y) &+ \epsilon \Big [n_1(y)\coth (y)\cos (k[x-c t]) \nonumber \\ & - 2 i \tilde {S}_1(y)\tanh (y)\sin (k[x-c t]) \Big ] \Big ) \exp (-it), \end{align}

from which the real and imaginary parts can be read. The corresponding real-valued solutions for the number density and phase are

(4.20a) \begin{align} n &= \tanh ^2(y) + 2 \epsilon n_1(y) \cos (k[x-c t]), \end{align}

(4.20b) \begin{align} S &= -t - 2\epsilon \tilde {S}_1(y) \sin (k[x-c t]). \end{align}

Equations (3.6) and (4.20b ) then give the velocity components,

(4.21a) \begin{align} u &= -2 \epsilon \tilde {S}_1(y) k \cos (k[x-c t]), \end{align}

(4.21b) \begin{align} v &= -2\epsilon \frac {{\textrm d}\tilde {S}_1}{{\textrm d}y}\sin (k[x-c t]). \end{align}

These linear solutions represent supersonic waves that travel without change of form, at the wave speed, $c_{\textit{min}}(k)$ , with respect to the stationary fluid at $y\to \infty$ . The corresponding solutions for a uniform base flow of velocity, $U \hat {\boldsymbol{x}}$ , follow directly from the transform in (4.2) and the text that follows that equation.

In summary, the linear stability analysis establishes that the base flow is marginally stable. That is, small/infinitesimal perturbations to the base flow neither grow nor decay in time. If vortices are to be nucleated, they must arise from finite-amplitude perturbations to this base flow, i.e. due to a nonlinear mechanism. To proceed, we therefore investigate the stability of finite-amplitude perturbations to the base flow using a DNS of the full GPE, (3.1), which includes all nonlinear effects.

Figure 6.

Eigenfunctions for wavenumbers, $k=0.5, 0.75,1.0,1.5, 5,10$ (top to bottom), computed using the numerical method. (a) Density function, $n_1$ , (b) the imaginary part of the phase function, $\tilde {S}_1$ , defined in (4.18).

5.1.1. Operator splitting

Operator splitting is applied to (3.1), with respect to the Cartesian coordinate system, requiring the sequential solution of two equivalent differential equations,

(5.1) \begin{equation} \frac {\partial \psi }{\partial t} = \frac {i}{2} \frac {\partial ^2 \psi }{\partial x^2} - i |\psi |^2 \psi , \qquad \frac {\partial \psi }{\partial t} = \frac {i}{2} \frac {\partial ^2 \psi }{\partial y^2}, \end{equation}

in the $x$ - and $y$ -directions, respectively; the external potential, $V_{\textit{ext}}$ , in (3.1) is set to zero. The corresponding boundary conditions are given in (3.2) and (3.3). The mapping,

(5.2) \begin{equation} \psi = \left (\phi + \frac {y}{y_{\textit{max}}}\right ) \exp (-i t), \end{equation}

is applied to the second differential equation in (5.1) and its boundary conditions, which become

(5.3) \begin{equation} \frac {\partial \phi }{\partial t}= i \left ( \frac {1}{2} \frac {\partial ^2 \phi }{\partial y^2} + \phi + \frac {y}{y_{\textit{max}}} \right )\!, \end{equation}

with homogenous boundary conditions,

(5.4) \begin{equation} \phi |_{y \to 0^+}=\phi |_{y \to \infty } =0. \end{equation}

This simplifies numerical solution of the second differential equation in (5.1), where the far field boundary condition in (3.3) is enforced at $y = y_{\textit{max}}$ .

5.1.2. Numerical algorithm for each differential equation

A standard split-time Fourier method is employed to solve the first equation in (5.1) over a time step, $\Delta t$ , yielding

(5.5) \begin{align} \psi (x,y, t+\Delta t) = \mathcal{F}_x^{-1} \Bigg \{\! \exp \left (\! -\frac {i}{2} K_x^2 \Delta t \!\right) \mathcal{F}_x \!\left \{ \exp \left (\! -i \big | \psi (\boldsymbol{r}, t) \big |^2 \Delta t \!\right) \psi (\boldsymbol{r}, t)\! \right \} \!\!\Bigg \} + O \big ( \Delta t^3 \big ), \end{align}

where $\mathcal{F}_x$ denotes the (spatial) discrete Fourier transform with respect to $x$ , $\mathcal{F}_x^{-1}$ is its corresponding inverse transform and $K_x = 2 \pi / J$ is the discrete wavenumber in the $x$ -direction (which is not to be confused with the wavenumber, $k$ , used in the linear stability analysis). Equation (5.5) is applied to each value of $y_m$ , for all $m=1,2,\ldots ,M-1$ , which then provides the initial condition at non-zero time, $t$ , for the second differential equation in (5.1) and its mapped version in (5.3).

Equations (5.3) and (5.4) are then solved using the discrete sine transform in $y$ together with a Crank–Nicolson integration in $t$ . This is performed for each value of $x_j$ , for all $j=1,2,\ldots ,J$ . This achieves spectral accuracy in $y$ and second-order accuracy in $t$ , identical to that of the $x$ integration rule in (5.5).

5.1.3. Initial condition at time, $t=0$

The required initial condition for the DNS at $t = 0$ , is chosen to be the solution arising from the linear stability analysis in (4.19). Equation (4.19) is discretised on the grid used for the DNS. This is performed for each fixed value of the wavenumber, $k$ , and the perturbation amplitude, $\epsilon$ . Use of this discrete initial condition enables the marginal stability prediction of the linear stability analysis to be assessed, as a function of the amplitude, $\epsilon$ . We note that completeness of the basis set of eigenfunctions in (4.19) is yet to be established. This limitation of the present approach is an avenue for future work.

Figure 7.

Linear-wave initial condition given by (4.19), $k =2$ . Solution at $t =5.658$ showing profiles of $n^{1/2} =|\psi |$ at $y = 0.7031$ . Solid, DNS with resolution $192\times 256$ ; dot–dashed, DNS with resolution $128\times 128$ ; dashed, linear solution, (4.19). Here (a) $0\leqslant x \leqslant 2\pi$ . (b) Close up of a part of solution. Differences between the two DNS solutions are of order $10^{-4}$ .

Figure 8.

Linear-wave initial condition given by (4.19), $k =2$ , $\epsilon = 0.02$ . Solution at $t =5.658$ showing profiles of $n^{1/2} =|\psi |$ at $y = 0.3125$ . Solid, DNS with resolution $192\times 256$ ; dot–dashed, DNS with resolution $128\times 128$ ; dashed, linear solution, (4.19). Here (a) $0\leqslant x \leqslant 2\pi$ . (b) Close up of a part of solution. Differences between the two DNS solutions are of order $10^{-4}$ .

5.3.1. Vortex nucleation at the wall $y=0$ with $\epsilon = 0.45$

Here DNS with $k=2, 3, 5, 10$ and $\epsilon = 0.45$ show different behaviour. These are done with different resolutions and domain heights as listed in table 2. Figure 11 shows three-dimensional surface plots of $|\psi (x,y)|$ with $k=2$ , $y_{\textit{max}} = 10$ , over the simulation domain for times $t=1.0, 1.2, 2.8, 4.0$ with resolution $(J\times M)=(192\times 256)$ . Figure 12 shows similar plots for $k=2, 3, 5, 10$ with $y_{\textit{max}} =20,\epsilon = 0.45$ , $ t= 4.0$ , all with resolution $(J\times M,y_{\textit{max}})=(128\times 256,20)$ .

Table 2.

Table of simulation parameters for $\epsilon = 0.45$ .

Figure 11.

Surface plots of $|\psi (x,y)|$ : (a) $t=0$ , (b) $t=1.2$ , (c) $t=2.8$ , (d) $t=4.0$ . $(J\times M)=(192\times 256)$ , $y_{\textit{max}}= 10$ .

Figure 11. Long description

Panel A: A 3D surface plot showing the variation of a parameter labeled psi as a function of x and y. The x-axis ranges from 0 to 6, the y-axis ranges from 0 to 10, and the z-axis represents psi ranging from 0 to 1. The surface exhibits a smooth gradient with a noticeable curvature. Panel B: Another 3D surface plot similar to Panel A, with slight variations in the surface curvature and gradient. Panel C: A 3D surface plot with more pronounced undulations and variations in the surface compared to Panels A and B. Panel D: A 3D surface plot with distinct surface features, including multiple peaks and valleys, indicating complex variations in the parameter psi across the x and y axes.

A striking feature of figure 11 is that for $t=2.8$ , $t=4.0$ , isolated zeros $|\psi (x,y)|=0$ in the $(x,y)$ plane have appeared. A zero in $|\psi (x,y)|$ requires both $\text{Re} [\psi ]=0, \text{Im} [\psi ]=0$ at the same $(x,y)$ . It will be shown later that, in the superfluid interpretation, the appearance of zeroes in $\psi$ corresponds to vortex nucleation. At both $t=2.8, 4.0$ , two zeros appear in the flow domain separated by $\pi$ in the $x$ -coordinate. This is because, over a $2 \pi$ , $x$ -periodic domain with $k=2$ , the spatial initial condition is $\pi$ periodic in $x$ . For visualisation purposes we plot $|\psi |$ rather than the density $n=|\psi |^2$ . Although $n$ is smooth ( $C^\infty$ ) across the computational domain, it vanishes quadratically in the vicinity of nodal points of $\psi$ , which reduces visual contrast in surface representations. By contrast, $|\psi |$ vanishes linearly with distance from a simple zero, producing a more pronounced local minimum and thereby providing a clearer diagnostic of vortex cores (to be discussed).

Figure 12.

Surface plots of $|\psi (x,y)|$ from DNS with $(J\times M)=(128\times 256)$ , $\epsilon = 0.45$ , $y_{\textit{max}} =20.$ , $t = 4.0.$ Here (a) $k=2$ , (b) $k=3$ , (c) $k=5$ , (d) $k=10.$ With $J=128$ fixed for all simulations, the number of $x$ -points per wavelength reduces with increasing $k$ as $J_w \approx 128/k$ . This can be seen in the decreasingly resolved region surrounding each $|\psi (x,y)|=0$ location as $k$ increases.

Comparison of figure 11(d) with figure 12(a) indicates that the effects of grid resolution and domain height on the near-wall flow structure are small for $k=2$ . This is also illustrated in figure 13(a,b), which shows profiles of $(\text{Re} [\psi ],\text{Im} [\psi ])$ versus $x$ on three values of $y$ at $t = 2.8$ and with $k=2$ , $\epsilon = 0.45$ . Each curve in figure 13 superposes the results of four simulations with $k=2$ (see table 2) displayed to show the effects of both resolution and domain height, which are small. Simulations with $(J\times M,y_{\textit{max}})=(128\times 256,20)$ and $(J\times M,y_{\textit{max}})=(128\times 128,10)$ differ only in the height $y_{\textit{max}}$ of the $y$ -domain. For the three profiles shown, the mean absolute difference in the calculated $|\psi |$ across the $x$ -domain is $O(10^{-4})$ , showing that, at least up to the nucleation time, the effect of the domain height is small. Figure 14(a) shows the power spectrum $|\psi |^2_{k_x}$ for a domain slice at $y = 0.89535$ for DNS with $(J\times M,y_{\textit{max}})=(192\times 256,10)$ and $(J\times M,y_{\textit{max}})=(128\times 128,10)$ , $k =2,\,\epsilon = 0.45$ and $t = 2.36$ while figure 14(b) shows corresponding $|\psi (x,y=0.89535,t=2.36)|^2$ profiles. Spectral convergence is clear.

For integer $k$ over a $2\pi$ $x$ -wise domain with an initial condition given by (4.19), and with a sufficiently large $\epsilon$ , we would therefore expect the simultaneous nucleation of $k$ vortices. This is illustrated in the four surface plots of figure 12 where, with $k=2, 3, 5, 10$ , $k$ zeroes can be seen in each panel. With a $2\pi$ $x$ -wise domain and with $J=128$ , we have $128/k$ points per wavelength, which accounts for the decreasing resolution of the large gradients surrounding each zero, as $k$ is increased. We found no indication of symmetry breaking in the $k$ -periodicity, in the present simulations.

Figure 13.

(a) Here $\text{Re} [\psi |$ versus $x$ , $t=2.8$ , $k=2$ , $\epsilon = 0.45$ for three values of $y$ . At $x=0$ , top to bottom; $y= 0.3125, 0.7813, 1.7187$ . (b) Here $\text{Im} [\psi |$ versus $x$ , $t=2.8$ , $k=2$ , $\epsilon = 0.45$ for three values of $y$ . At $x=0$ , bottom to top; $y= 0.3125, 0.7813, 1.7187$ . Each curve contains results from four DNS results superimposed. Solid line, $(J\times M,y_{\textit{max}})=(192\times 256, 10)$ ; dashed, $(J\times M,y_{\textit{max}})=(128\times 256, 20)$ ; dot–dashed, $(J\times M,y_{\textit{max}})=(128\times 128, 10)$ ; dotted, $(J\times M,y_{\textit{max}})=(64\times 64, 10$ ).

Figure 14.

(a) Here $x$ -Spectra $|\psi |_{k_x}^2$ at $y = 0.89535$ with $k=2$ , $\epsilon = 0.45$ , $y_{max}=10$ , $t =2.36$ . Resolutions; $(192\times 256)$ (crosses), $(128\times 128)$ (circles). The compressed range in $k_x$ is a result $\pi$ periodicity in the $2\pi$ domain with $k=2$ . (b) Here $|\psi |^2$ versus $x$ . Solid $(J\times M)=(192\times 256)$ ; dashed $(J\times M)=(128\times 128)$ .

Figure 15.

(a,c,e,g) Level sets of zero $\text{Re} [\psi ], \text{Im} [\psi ]$ in $(x,y)$ plane. Solid line; $\text{Re} [\psi ]=0$ ; dashed line; $\text{Im} [\psi ]=0$ . (b,d, f,h) Corresponding contour plots of $|\psi |$ . Here (a,b), $t =1.8$ , (c,d), $t =1.92$ , (e, f), $t =2.1$ , (g,h), $t =2.2$ . Level sets collide on $y=0$ at $t\approx 1.92$ .

Vortex nucleation is diagnosed further in figure 15. Figure 15(a,c,e,g) show zero level sets of $(\text{Re} [\psi ]=0, \text{Im} [\psi ]=0$ ) in the $(x,y)$ plane at the times indicated. This is the same plot style as for figure 10. Both $\text{Re} [\psi ]=0$ and $Im[\psi ]=0$ must form $90^\circ$ intersections with the solid wall. To show this, we can write $\psi = n^{1/2}(\cos [S(x,y,t)] + i\sin [S(x,y,t)])$ . Hence, $\text{Re} [\psi ]=0$ on the curve $x=x(y)$ at some $t$ , and with $n\ne 0$ , requires $S(x(y),y,t)= (n+1/2) \pi$ with $n \in (-\infty ,\infty )$ . Differentiating with respect to $y$ , we find that on  $\text{Re} [\psi ]=0$ ,

(5.6) \begin{equation} \left .\frac {{\textrm d}x}{{\textrm d}y} \right |_{\text{Re} [\psi ]=0} = -\frac {\frac {{\textrm d}S}{{\textrm d}y}}{\frac {{\textrm d}S}{{\textrm d}x}} = -\frac {v}{u}. \end{equation}

Since as $y\to 0$ , $v\to 0$ while $u$ is generally finite, then ${\textrm d}x/{\textrm d}y \to 0$ as $y\to 0$ at the solid wall. The same argument can be applied to $\text{Im} [\psi ]=0$ . This property is clear from figure 10. We expect this to be true for a general smooth boundary shape.

In figure 15, there is some jaggedness in $\text{Re} [\psi ]=0$ because of the interpolation used on the discrete set of grid points together with the shallow zero crossing. Nonetheless, it is clear that $(\text{Re} [\psi ]=0, \text{Im} [\psi ]=0)$ intersect the solid wall normally. Their zero level sets first collide or intersect on the singular line, $y=0$ (recall that this also corresponds to $|\psi (x,0)|=0$ ) at $t_c\approx 1.92$ , which will henceforth be referred to as the wall-zero-level-set collision time. This event marks the appearance of a zero in $|\psi |$ for the present initial-value problem.

Both prior to, at and after collision, the two level sets show equal and opposite curvature in the $(x,y)$ plane. In other words, the collision event involves mutually opposed curvature in the $(\text{Re} [\psi ]=0, \text{Im} [\psi ]=0)$ profiles. This is consistent with collision occurring at, rather than above the solid wall. For $t\gt t_c$ this geometrical property ensures that the $|\psi (x,y)|= n^{1/2}$ zero formed at collision moves off the $y$ -axis and convects with the superfluid flow. Figure 15(b,d, f,h) show corresponding contours of $|\psi (x,y)|$ . The detached zeros are clear at $t=2.1, 2.2$ . Isolated zeros in $n(x,y,t)$ correspond to the presence of discrete vortices in the superfluid, see Appendix C. Supplementary movie 1 shows an animation of the time evolution of the number density distribution, where the nucleation of zeros/vortices near the wall and their subsequent motion is evident. Note that the streamlines in the vicinity of each vortex are circular (not shown) while the level sets of the number density distribution are generally non-circular, see Appendix C. The connection between these zeros and the emergence of non-zero circulation is explored subsequently.

5.3.2. Gross–Pitaevskii invariant energy

As a further diagnostic on solution accuracy, we briefly discuss calculation of the invariant GPE energy (Gérard Reference Gérard2006), which for the present flow over one spatial period may be expressed as

(5.7) \begin{equation} E[\psi ](t) = \frac {1}{2}\int _0^{2\pi }\int _0^\infty \left [|\boldsymbol{\nabla }\psi |^2 + \big(|\psi |^2 -1 \big)^2\right ]\,{\textrm d}y\,{\textrm d}x. \end{equation}

This can be shown to satisfy ${\textrm d}E[\psi ](t)/{\textrm d}t=0$ (see Appendix D). For computation we here replace the upper $y$ integration limit by $y_{\textit{max}}$ , due to the imposed finite computational domain. To calculate the kernel in (5.7), $x$ -derivatives are calculated spectrally with fourth-order FDs for the $y$ -derivatives. The $x$ -integration uses the periodic trapezoidal rule while Boole’s rule is used for the $y$ -integration.

With $\epsilon =0$ , the base-flow solution (4.3) gives the exact analytical result, $E[\psi ]= 4\pi /3$ . A full $\epsilon =0$ , GPE simulation over $0 \leqslant t \leqslant 4$ was performed as a test case using this benchmark. The integration algorithm calculates $E[\psi ](t=0)$ accurately to machine precision. For $t\gt 0$ the energy drift $\Delta E[\psi ](t) =|(E[\psi ](t) -E[\psi ](0))/E[\psi ](0)|$ is calculated as a numerical solution diagnostic. Its time mean is ${\overline {\Delta E[\psi ]}} \approx 5\times 10^{-6}$ owing to time-integration errors. With $\epsilon = 0.45$ and $k=2$ , we find ${\overline {\Delta E[\psi ]}} \sim 10^{-5}$ independent of spatial resolution. The time variation of $\Delta E[\psi ](t)$ is shown in Appendix D. Similar results are found for other cases.

5.3.3. Superfluid velocity and domain circulation, $\varGamma (t)$

A search across all DNS time frames with $k=2$ and initial conditions using $\epsilon = 0.02, 0.1, 0.3$ show no evidence of points in the flow domain where both $(\text{Re} [\psi ]=0,{} \text{Im} [\psi ]=0)$ . That the appearance of these zeroes is the signature of vortex nucleation can be demonstrated by calculation of the flow velocity component, $u$ , on $y=0$ , and from this, the circulation in the domain. The superfluid flow velocity field, $(u(x,y,t),v(x,y,t))$ , can be calculated at each point using (3.10a )–(3.10b ). This is done spectrally for $u$ and with fourth-order FDs for $v$ . From (3.10a )–(3.10b ) it is clear that the flow velocity is singular where $|\psi |=0$ . This will be the case both at isolated zeroes of $|\psi (x,y)|= n^{1/2}$ , and also along the entire solid wall, $y=0$ . As shown in § 3.3, the latter is, however, generally a removable singularity for the velocity but one that cannot be calculated directly from the DNS. At each wall point $x_j$ , $u(x_j,0,t)$ is calculated with polynomial extrapolation to $y=0$ using the first three off-wall points, resulting in subspectral accuracy in the wall velocity. This is minor without consequences for the present study.

Figure 16.

(a) Velocity component $u(x,0,t)$ versus $x$ near wall location and time where vortex nucleation occurs. Dotted, $t=1.7$ ; dashed, $t= 1.8$ ; solid, $t=2.0$ ; dot–dashed, $t =2.2$ . (b) Circulation $\varGamma (t)/(4 \pi )$ versus $t$ . Parameters $ k=2, \epsilon = 0.45$ . A discontinuity in $\varGamma (t)$ occurs at $t = 1.992\pm 0.002$ , signalling vortex nucleation.

Figure 16(a) shows $u(x,0,t)$ calculated at $t = 1.7, 1.8, 2.0, 2.2$ over $0\leqslant x \leqslant 2$ . In the velocity plot, the first two values are for $t \lt t_c$ . Each shows a large negative peak. The values at $t = 2.0, 2.2$ are after level-set collision and each exhibits a large positive peak. Velocity distributions just before and just after $t_c=1.92$ can be calculated numerically but show narrow peaked behaviour. With $k=2$ , similar velocity profiles are found when $x$ is shifted by $\pi$ in the positive direction.

Figure 17.

(a) Instantaneous streamlines, $k=2, \epsilon = 0.45, t = 1.92$ . The beginning of vortex nucleation can be seen near $x = 0.65, y =0$ and near $x = 0.65 +\pi , y =0$ . (b) Instantaneous streamlines, $k=2, \epsilon = 0.45, t = 4.0$ . One nucleated vortex can be identified near $x = 3.14, y =1.51$ .

The domain circulation $\varGamma (t)$ can be calculated from the phase $S =\arctan (\text{Im} (\psi )/{}\text{Re} (\psi ) )$ . In superfluid variables $S$ is the velocity potential, for which changes around a closed curve $C$ gives the enclosed circulation $\varGamma (t) = \oint _C \boldsymbol{\nabla }S {\boldsymbol{\boldsymbol{\cdot }}}{\textrm d}{\boldsymbol{\ell }} =\oint _C(\partial S/\partial \sigma )\, {\textrm d}\sigma = \Delta S$ where $\sigma$ is the arclength and $\Delta S$ is the change in $S$ around $C$ . Let $C$ be the computational domain. Contributions to the integral from periodic boundaries $x=0,\,2\pi$ and from $y=y_{\textit{max}}$ sum to zero.

Two implementations are used. In the first, $S(x_j,y= y_{max}/M,t)$ is calculated at the first wall-normal point at each $x$ location. The second uses a wall-estimate on $y=0$ , where $\psi =0$ and $S$ is a removable singularity which can be determined by L’Hôpital’s rule, giving $S = \arctan (\text{Im} (\partial \psi /\partial y)/ \text{Re} (\partial \psi /\partial y) )$ with derivatives at y = 0 calculated with sixth-order, centred FDs. Both produce a raw wall array $S(x_j,y=0,t)$ . Owing to periodic boundary conditions on $x=0,\,2\pi$ , $S(x_j,y=0,t)$ is periodic in $x$ (i.e. $S(2\pi ,0,t) = S(0,0,t)$ ) but may contain discontinuities. These occur because the principal branch of $\arctan (\theta )$ is defined in $(-\pi \lt \theta \leqslant \pi )$ . On the GPE- $\psi$ numerical grid, a non-physical discontinuity in $S$ occurs across two cells when $\text{Im} (\psi )$ changes sign with $\text{Re} (\psi )\lt 0$ . Discontinuities can be eliminated by detection and branch correction, thus creating an ‘unwrapped’ array – henceforth denoted $\hat {S}(x_j,y=0,t)$ – that is discontinuity free but may not be periodic (if $\varGamma (t) \ne 0$ ) in $0 \leqslant x \lt 2\pi$ .

The domain circulation is then $\varGamma (t) = \hat {S}(2\pi ,0,t) - \hat {S}(0,0,t)$ and is obtained at discrete times, $\delta t_s = 0.004$ . This is shown in figure 16(b) versus time for $k=2, \, \epsilon = 0.45$ using the second method described in the previous paragraph. Generally, a $\varGamma (t)$ discontinuity occurs across a single $\delta t_s = 0.004$ at the time $t_s =1.92\pm 0.02$ in figure 16(b). The uncertainty is estimated from comparing results for the two methods previous paragraph with differences owing to small changes in $\text{Im} (\psi )$ zero crossings, sensitivity to the branch-correction algorithm and differing $y$ -locations. Circulation is zero for $t\lt t_c$ and $\varGamma (t)/(4\pi )=1$ for $t\gt t_c$ , corresponding to $2\pi$ for each of the two vortices nucleated with $k=2$ .

Images of instantaneous streamlines are shown in figure 17 for $k=2, \epsilon = 0.45$ . Vortex nucleation can be seen in figure 17(a) at $t = 1.92$ near $x = 0.65, y =0$ and also near $x = 0.65 +\pi , y =0$ . In figure 17(b), one vortex has migrated to $x \approx 3.14, y \approx 1.51$ . The sense of anticlockwise or positive circulation is clearly visible.

5.3.4. Nucleation times, $k=2, 5$

As discussed earlier, vortex nucleation with $k=2$ is not observed up to $t=25$ , for $\epsilon = 0.3$ . Times-to-nucleation for $\epsilon = 0.325, 0.35, 0.4, 0.45$ for both $k=2,5$ were estimated using the calculation of $\varGamma (t)$ for these simulations, and locating discontinuity times, where this is observed. With $\epsilon = 0.3,\,0.325$ , we find $\varGamma = 0$ up to $t =25$ . We thus interpret $\epsilon = 0.3$ as a lower threshold in $\epsilon$ for vortex nucleation with the present class of initial conditions, for both $k=2,5$ over the present simulation time. Values of $t_c$ for other $\epsilon$ are shown in table 3 again with uncertainty estimates as described above. These are somewhat smaller for $k=2$ , $\epsilon = 0.45$ because this simulation is done with $(J\times M) = (192\times 256)$ compared with $(J\times M) = (128\times 256)$ for other cases, and with a smaller $\Delta t_s$ . We remark that for a given $k$ we cannot rule out vortex nucleation at $\epsilon$ -values lower than those of table 3, for simulation times greater than $t = 25$ .

Table 3.

Values of time to vortex nucleation, $t_c$ , for several values of $\epsilon$ , $k$ .

The values $t_c=0$ with $\epsilon = 1/a_2$ in table 3 can be justified as follows. Consider the initial condition, (4.19). with $t=0$ . The initial velocity field in the whole domain can then be calculated using (3.10). Near $y=0$ we can write $n_1(y) = a_2 y^2 + a_4 y^4 + \ldots$ and $S_{1,i}(y) = b_0 + b_2 y^2 + \ldots$ where the constants $a_2(k), a_4(k), \ldots , b_0(k) \ldots$ are obtained numerically from the linear stability calculations described previously. Substituting these expressions into (4.19), and then the resulting equation for $\psi$ into (3.10) gives the velocity components near the solid wall for the initial condition.

The initial velocities on $y=0$ can then be determined by the limit $y\to 0$ with the result $v =0$ and

(5.8) \begin{equation} u(x,0,\epsilon ) = - \frac { 2 b_0 \epsilon k (\epsilon a_2 + \cos (kx)) }{(1 + \epsilon a_2 \cos (kx))^2+ 4\epsilon ^2 b_0^2 \sin ^2(kx) }. \end{equation}

If $a_2\gt 0$ and $\epsilon = 1/a_2$ , (5.8) contains removable singularities at $x = (2m+1)\pi /k$ . Recognising this, with $\epsilon = 1/a_2$ , the solid wall velocity component, $u_{\textit{bg}} \equiv u |_{y \to 0^+}$ , can then be evaluated as

(5.9) \begin{equation} u_{\textit{bg}}(x,k) = - \frac {2 b_0 a_2 k }{a_2^2 + 4 b_0^2 + \big(a_2^2 - 4 b_0^2\big) \cos (k x) }, \end{equation}

which we refer to as the background initial wall speed for $\epsilon = 1/a_2$ . This is analytic for all parameters. The limits, $\epsilon \to 1/a_2$ , from below and above are, however, singular. This is now discussed.

5.3.5. Singular limit

Figure 18(a) shows the variation of $a_2(k)$ and $b_0(k)$ obtained from the earlier described solution of the linear stability problem. These show little variation with $k$ for $k\gt 2$ and also $b_0\to - a_2/2$ for large $k$ . It is noted that $a_2 \approx 2$ and $b_0 \approx -1$ . For $k=2$ , $a_2 = 2.007, b_0 = -1.102$ giving $\epsilon = 0.4981$ . Figure 18(b) shows wall-parallel $x$ -velocities at the solid wall given by (5.8) with these $(a_2,b_0)$ for $\epsilon = 0.3, 0.45,1/a_2-10^{-3}$ . As $\epsilon \to 1/a_2$ , the velocity profiles become more singular and appear to indicate approach to Dirac delta function-like behaviour near $x=\pi /2, 3\pi /2$ with an almost constant background with magnitude unity away from these points. This is presently taken to indicate vortex nucleation at $t=0$ as $\epsilon$ passes through $1/a_2$ , which is confirmed by the DNS. Hence, we can take $t_c=0$ at $\epsilon = 0.4981$ for $k=2$ , and $\epsilon = 0.5046$ for $k=5$ .

Figure 18.

(a) Variation of parameters with $k$ from linear analysis. Solid line, $a_2$ ; dot–dashed, $- a_4$ ; dashed, $-b_0$ . (b) Wall $x$ -velocities at $t=0$ , $a_2 = 2.007, b_0 = -1.102$ , $k=2$ . Dot–dashed, $\epsilon = 0.3$ ; dashed, $\epsilon = 0.45$ ; solid, $\epsilon = 1/a_2-10^{-3}$ .

This is analysed further in Appendix E where it is shown that, within the subdomain $0\leqslant x \lt 2\pi /k$ , when $\epsilon \to 1/a_2$ from both below and above, then

(5.10a) \begin{align} u \!\left (x,0,k,\epsilon \to (1/a_2)^-\right ) \to u_{\textit{bg}}(x,k) -\pi \delta \! \left (x-\frac {\pi }{k}\right ), \end{align}

(5.10b) \begin{align} u \!\left (x,0,\epsilon \to (1/a_2)^+\right ) \to u_{\textit{bg}}(x,k) +\pi \delta \! \left (x-\frac {\pi }{k}\right ), \end{align}

where $\delta (\ldots )$ is the Dirac delta function, $u_{\textit{bg}}(x,k)$ is given by (5.9). Taking $a_2=2, b_0 = - 1$ then $u_{\textit{bg}}(x,k)=k/2$ . This is consistent with results shown in figure 18(b). Further, it can be shown that $\int _0^{2\pi /k}u_{\textit{bg}}(x,k) {\textrm d}x = \pi \textrm{sgn}(-b_0)$ which equals $\pi$ since we find $b_0\lt 0$ presently. It follows that $\int _0^{2\pi /k}u(x,y=0,k) {\textrm d}x = 0, \epsilon \lt 1/a_2$ and is equal to $2 \pi$ when $\epsilon \gt 1/a_2$ .

5.3.6. Point vortex embedded in initial condition, $\epsilon \gt 1/a_2$

This last result is consistent with the appearance of a point vortex contained in the initial condition, (4.19), when $\epsilon \gt 1/a_2$ . To confirm this, we write (4.19) to $O(y^2)$ near $y=0$ as

(5.11a) \begin{align} \psi _{\textit{re}}(x,y,k,\epsilon ) &= \frac {1}{3}y\big (3 -y^2 +\epsilon \big (3 a_2+ (a_2 + 3 a_4) y^2 \big ) \cos (k x) \big ), \end{align}

(5.11b) \begin{align} \psi _{\textit{im}}(x,y,k,\epsilon ) &= \frac {2}{3}y \epsilon \big (-3 b_0+ (b_0 - 3 b_2) y^2 \big ) \sin (k x). \end{align}

For $x = (2 m+1)\pi /k$ , $\psi _{\textit{im}}=0$ while $\psi _{\textit{re}}=0$ at $y=0$ and at $y$ satisfying $3 -y^2 -\epsilon (3 a_2+ (a_2 + 3 a_4) y^2 ) = 0$ . Solving this equation, putting $\epsilon = 1/a_2 + \hat {\epsilon }$ where $\hat {\epsilon }\ll 1$ and expanding to leading order in $\hat {\epsilon }\ll 1$ gives the solution

(5.12) \begin{equation} y_{\hat {\epsilon }} = \frac {\sqrt {3 \hat {\epsilon }} a_2}{ \sqrt {-(2 a_2 + 3 a_4) } }. \end{equation}

From the linear stability analysis described earlier we find $2 a_2 + 3 a_4 \lt 0$ (see figure 18 a) for all $k\gt 0$ , and so $y_{\hat {\epsilon }}\gt 0$ for $\hat {\epsilon }\gt 0$ . Hence, $\psi =0$ at $(x_v,y_v) = ((2 m+1)\pi /k,y_{\hat {\epsilon }})$ .

For the initial condition $\psi _{\textit{re}}+ i \psi _{\textit{im}}$ in the neighbourhood of $(x_v,y_v)$ to represent a point vortex it is sufficient that $\psi$ have finite gradients at $(x_v,y_v)$ , see Appendix C. Taylor expanding (5.11a ) and (5.11b ) about $(x_v,y_v)$ gives to leading order, after some algebra

(5.13) \begin{equation} \psi _{\textit{re}}(x^{\prime},y^{\prime}) + i \psi _{\textit{im}}(x^{\prime},y^{\prime}) = i \frac {2 \sqrt {3 \hat {\epsilon }} b_0 k}{ \sqrt {-(2 a_2 + 3 a_4) } } x^{\prime}+ 2 a_2 \hat {\epsilon } y^{\prime}, \end{equation}

where $x^{\prime} = x - x_v, y^{\prime}= y -y_v$ , which is a special case of (C5) of Appendix C. Hence, (5.13) represents a point vortex of strength $2 \pi$ at $(x_v,y_v)$ . We can conclude that for $\epsilon \gt 1/a_2$ , the initial condition (4.19) contains a point vortex. This is consistent with numerical simulations that show $t_c =0$ at $\epsilon = 1/a_2$ .

5.3.7. Point-vortex trajectory, $k=2$

Once nucleated, the two point vortices convect with the base flow. Their $(x(t),y(t)$ trajectories for $t\gt t_c$ can be determined numerically by locating the position of the $|\psi |=n^{1/2}$ isolated zero in the GPE solution grid at each discrete time. This is done at selected time intervals by fitting a second- or third-order interpolating polynomial to both $(\text{Re} [\psi (x,y], \text{Im} [\psi (x,y)])$ at a given time list, and then calculating the intersection of the respective zero level sets. The $y(t)$ trajectory is displayed in figure 19(a) in linear coordinates. Figure 19(b) shows $y(t)$ versus $(t-t_c)$ with $t_c=1.92$ in log–log coordinates. Also shown in both figures is the line of best power-law fit over the first $12$ points as $y(t) = 1.28 (t-1.92)^{0.46}$ . The vortex vertical speed is infinite at $t=t_c$ . The $x(t)$ trajectory is almost a straight line with speed close to the linear wave speed for $k=2$ .

Figure 19.

(a) Vortex $y$ trajectory versus $t$ . (b) Vortex $y$ trajectory versus $t$ in log–log coordinates. Dashed line in (a) and solid line in (b) $y(t) = 1.28 (t-1.92)^{0.46}$ .

The initial $y(t)$ trajectory can be understood by reference to our earlier observation that vortex nucleation is found to be initiated by the collision of level sets $\text{Re} [\psi ]=0$ , $\text{Im} [\psi ]=0$ at the solid wall, and further, that the level sets terminate normally at the solid wall and with opposite-signed curvature, see figure 15(c). Consider modelling the level-set collision by the independent movement of two ellipses with semiminor axes aligned with the $x$ -axis. Suppose that these are distinct for $\tau \lt 0$ but begin a passive intersection starting at $\tau =0$ . It is straightforward to show that the $y$ trajectory of the initial intersection point varies as $\tau ^{1/2}$ for small $\tau$ . As a simple example, consider circles of radius $r_0$ centred at $x = \pm r_0$ at $\tau =0$ and moving with respective speeds $-/+ U$ . For small $\tau$ , their intersection point is $y= (2U \tau r_0)^{1/2}\sim \tau ^{1/2}$ in agreement with the numerics. The infinite ejection speed is strictly kinematical.