Modelling turbulent flow of superfluid $^4$He past a rough solid wall in the $T = 0$ limit

We present a numerical study, using the vortex filament model, of vortex tangles in a flow of pure superfluid $^4$He in the $T = 0$ limit through a channel of width $D = 1$ mm for various applied velocities $V$. The flat channel walls are assumed to be microscopically rough such that vortices terminating at the walls are permanently pinned; vortices are liberated from their pinned ends exclusively through self-reconnection with their images. Sustained tangles were observed, for a period of 80 s, above the critical velocity $V_c \sim 0.20$ cm s$^{-1} = 20 \kappa/D$. The coarse-grained velocity profile was akin to a classical parabolic profile of the laminar Poiseuille flow, albeit with a non-zero slip velocity $\sim$ 0.20 cm s$^{-1}$ at the walls. The friction force was found to be proportional to the applied velocity. The effective kinematic viscosity was $\sim 0.1\kappa$, and effective Reynolds numbers within $\mathrm{Re'}<15$. The fraction of the polarized vortex length varied between zero in the middle of the channel and $\sim$ 60% within the shear flow regions $\sim D/4$ from the walls. Therefore, we studied a state of polarized ultraquantum (Vinen) turbulence fuelled at short lengthscales by vortex reconnections, including those with vortex images due to the relative motion between the vortex tangle and the pinning rough surface.


Introduction
Flow of superfluid 4 He through a channel is only non-dissipative when its velocity is below the channel-specific critical velocity v c .Above it, the dissipative regime sets in, in which the chaotic motion of quantum vortices (called Quantum Turbulence) results in the transfer of the flow momentum to the channel's walls, i. e. a friction force.
Experimental evidence has shown that this friction force becomes greatly reduced for T < 0.7 K [1] and that vortex pinning becomes much weaker with lowering temperature below T < 0.4 K [2].Approaching T = 0, the density of viscous normal component vanishes and the mutual friction which couples it to the vortex tangle of turbulent superfluid is negligible, thus the direct interaction between vortex lines and irregularities of the walls of the channel (vortex pinning) must be considered.
In these computer simulations, we model a flow of superfluid 4 He at T = 0 between two parallel solid walls.We assume the limit of extremely rough walls, where the areal density of sharp protuberances is greater than the density of vortex lines in the vortex tangle.Then the processes of un-pinning (due to a self-reconnection with an image vortex) of each vortex line from the pinning protuberance and re-pinning at the nearest protuberance occur on the lengthscale smaller than the characteristic scale of the vortex tangle (of order mean inter-vortex distance ℓ), i. e. independently of other vortices.This is different from the the generation of dense vortex tangles by larger-scale irregularities of the profile of solid wall simulated by Stagg et al. [17].
In numerical simulations vortex lines are represented by chains of discrete points with the set inter-point scale δ, the smallest resolved lengthscale is ∼ δ.Within the mechanism of pinning-unpinning mentioned above, the elementary distance between the unpinned-pinned ends of a vortex line becomes of order ∼ δ.Hence, in our framework, δ is a proxy of the scale of roughness of the solid wall.

Numerical Methods
The VFM [18] has been used to great effect to simulate and visualise the dynamics of vortices in superfluid helium.In the limit of zero temperature, using the VFM, the local self-induced vortex velocity ṡ evaluated at the point r along a vortex line may be described entirely by the Biot-Savart integral where κ = h/m 4 is the quantum of circulation in superfluid 4 He.The line integral over L represents inclusion of the complete vortex configuration, which is discretized into points s i for (i = 1, ..., N ).The discretized integral, after removing the singularity by separating the local and non-local contributions, becomes where l + and l − represent the arc lengths to adjacent vortex points s i+1 , and s i−1 , s ′ i and s ′′ i are the local tangent and curvature respectively, and a ∼ 1 Å is the size of the vortex core in He-II.Here the first term is akin to the local induction approximation (LIA) and the second term describes the non-local contributions, including the effect of image-vortices due to any solid boundaries present [9].To reduce the number of computations required to simulate vortex tangles in a timely manner, the integral has been reduced to only include contributions from points within a cut-off radius r c of the target point s i , defining a sphere of contributing non-local vortex points as shown in Figure 1 (right).The cut radius was chosen to be half the container size r c = D/2.
The simulation volume was a cubic cell D × D × D with periodic boundary conditions in the x, y directions.The boundaries at z = 0 and z = D were modelled as a rough solid plane surface with strong pinning.To satisfy the solid boundaries the method of images allows the duplication and reflection of the vortex configuration across each solid boundary.This combined with a periodic wrapping of the system at the x, y boundaries gives 26 copies of the original volume: 2 reflected at solid boundaries, 8 simply periodic and 16 periodic-reflected copies.An illustration of the boundary conditions in the x and z directions is given in Fig. 1 (left).All 27 cubes were considered in evaluating Equation 2. The spatial resolution δ of the simulation is preserved by removing points on the same vortex line that move within δ/2 of each other.Similarly, new points are introduced, according to the local curvature of the line, when adjacent points are more than δ apart.Additionally, if two non-adjacent points on the same vortex line or otherwise, become closer than the critical distance δ/2, a line-line reconnection occurs, altering the vortex configuration and causing dissipation from the removal of vortex line length.The reconnection procedure used was of Type-III [19].
To produce channel flow in the T = 0 limit a positive finite superfluid velocity V is applied in the x direction such that instantaneous velocity of the discrete vortex point is The time-evolution followed a third order Adams-Bashforth scheme: where ∆t is the size of a time step and n = t/∆t is the current time in integer steps [20].
The presence of solid wall boundaries alters the vortex motion since no flow is allowed through the walls.The solid boundaries are assumed to be microscopically rough such that the ends of vortex lines, terminated at the walls, are permanently pinned and fixed in position with u(z = ±D/2) = 0. Liberation of the vortex line is only possible through its self-reconnection with its image following the same reconnection criterion that triggers a line-line event.Thus, a vortex point, which comes within a distance δ⁄4 of the wall, reconnects with its image, thus liberating the line before it becomes pinned again ∼ 0.7δ away, as demonstrated in Figure 2. Here, the reconnection was triggered by the blue line segment, which is then removed from the simulation alongside the adjacent red line segment.The new end of the vortex is then instantaneously pinned by artificially adjusting the z coordinate to coincide with the surface, which generates a kink in the filament just above the wall.The model mimics a vortex line "walking" along a flat rough surface, with the vortex end jumping between sharp protuberances spaced on the scale of the spatial resolution δ [21].
Effectively, this mechanism assumes the existence of a critical angle between the vortex and the wall for unpinning.Such an assumption was recently used in simulations of vortex lines attached to a microelectromechanical (MEMS) oscillator, with a chosen critical angle of θ c = π/6 [16].
The potential benefits of our approach are in that one could investigate the role of Kelvin waves emitted on vortex lines by discrete steps; although to assure an adequate resolution of those waves one might need to introduce a separate lengthscale for the kink generated by an image-reconnection such that the wavelength of injected distortions is larger than the resolution.
To examine the rate of momentum transfer to the walls, two equivalent methods of determining the friction force were used.Both have been validated by comparison with the analytical solution for a vortex semi-ring travelling along a flat rough surface.The rate of change of the total impulse of all vortex lines in the simulation volume, termed the integral method: where ρ s is the superfluid density and ξ is the arc length along a vortex line.The second method assumes that the vortex line tension f t (energy per unit line length) is constant.The angle θ between the surface and the pinned segment can be used to obtain the components of the line tension for a single terminated line.Thus, summing over all the pinned ends gives the total friction force in the streamwise direction where b was a cut-off length scale chosen as δ/2 which is the effective size of pinning sites in the vortex walking and a = 1 Å.We term this the tension method.
To examine the velocity profile, the instantaneous mean velocity as a function of the wall-normal direction z must be calculated.This was achieved by slicing the simulation volume into N slabs with thickness D/(N − 1), the mean velocity of all line segments within a single slab was then time-averaged during the steady state to obtain the coarse-grained velocity profile u(z) .The component in the direction of flow, u x then gives the cross-channel velocity profile.Fixed segments of the vortex line were omitted from such calculations to avoid skewing the average velocity, with u pinned = 0, in the top and bottom slabs.The variation in the coarse-grained vortex line density across the height of the channel L(z) was obtained similarly, but without the need to exclude fixed segments.

Results
The parameters for simulations were: D = 0.1 cm, δ = 2 × 10 −3 cm, with temporal resolution dt = 8 × 10 −5 s.Initially, at t = 0, the volume was populated with 80 randomly placed and randomly oriented vortex rings with equivalent radii of 0.012 cm; with an instantly acting applied superflow of V = 0.30 cm s −1 .Snapshots of the early development of the tangle is shown by the inset panels of Fig. 3.At early times (t < 1 s) frequent reconnections between vortex rings in the centre of the channel cause an immediate drop in the total vortex line length, Λ, of the system, which is also shown in Fig. 3. Gradually vortex rings become attached to the walls and begin "walking" across the surface in the direction of flow.Within t = 4 s both surfaces are decorated with several terminating vortices.As the flow develops further vortex lines become stretched along it causing a gradual increase in the total line length as a vortex tangle forms between the walls.Both line-line and image-line reconnections occur frequently, each process removing a small portion of vortex length, which becomes approximately stable when the rate of dissipating length from reconnections is similar to the rate of vortex growth at the walls.At t = 80 s the vortex tangle is comfortably in the steady state with the total vortex length in the simulation volume Λ = 4.14 cm, and the walls are decorated with 14 pinned vortex lines.The tangle at t = 80 s is the same as the configuration shown in Fig. 1.
To simulate flows at any other applied velocity in the range 0.08 cm s −1 ≤ V ≤ 0.34 cm s −1 for further 80 s, the same initial configuration of the vortex tangle was used, namely, the one for V = 0.30 cm s −1 at t = 80 s.
The evolution of Λ for a selection of V is shown in Fig. 3 in the range 80 s ≤ t ≤ 160 s.Vortex tangles, sustained for at least 80 s as in Fig. 4, were observed above the critical applied flow velocity V c = 0.19 ± 0.01 cm s −1 , i. e. V c D/κ ∼ 20.Below V c , the tangle is not sustainable; all vortex lines eventually detach from the walls at some point and move with the flow without gaining any energy -hence the tangle decays.
Fig. 4 Snapshots of vortex tangles in the steady state at t = 160 s for all V > Vc.The flow is directed towards the right.
For each applied superfluid velocity V , the mean flow velocity was calculated as u x = D −1 D 0 u x (z)dz and is plotted against V in Fig. 5a.They are proportional to each other with relation u x ≈ 0.96V .
The stream-wise component of friction force was evaluated from equations ( 5) and ( 6), time-averaged over the steady state and normalised by the reduced vortex line tension f t = ρκ 2 4π ln b a , where b = δ/2 is the effective size of a pinning site.The observed force-velocity relation, shown by Fig. 5b, was approximately linear within the range of flow velocities selected.Friction force per number of pinned vortex ends was constant across all velocities at (0.626 ± 0.007)f t demonstrating that the total friction force is proportional to the number of pinned ends in the simulation box and the effective critical angle in the direction of flow is cos −1 (0.626) = 51 • .Both integral and tension methods of determining the friction force consistently agreed to within 1%.
The total vortex line length Λ fluctuates with the dominant period ∼ 10 s (corresponding to the time constant ∼ 2 s).This could be compared with the time D/V c ∼ 0.5 s for the mean flow to pass the cell size, the intrinsic quantum time constant D 2 /κ = 10 s and our chosen observation time of 80 s.Fig. 5c demonstrates the variation of time-averaged vortex line density, L = Λ /D 3 , which appears to be cubic with velocity.Interestingly, the values of L for unsustained tangles do not fall far from the V 3 line.
The coarse-grained stream-wise velocity profile u x (z) was steady and nearly parabolic albeit with a non-zero slip velocity, as seen in Fig 6 .We fit the profile with the form of a classical Poiseuille profile, adapted for finite slip boundary conditions at the walls, where u min represents the slip-velocity at the walls and u max is the maximum local velocity in the centre of the channel.The agreement is relatively good, but the fits consistently underestimate u max .The slip-velocity varied only weakly across the simulations with values in the region of V c as seen in Fig. 7 (left).The universality of the flow profiles is shown in Fig. 6 (right).Deviations from the universal profile were only observed for the lowest V which still produced a steady flow -possibly due to the smaller number of line segments leading to less averaging off of fluctuations.The effective kinematic viscosity ν ′ = F 8D(umax−umin)ρ was in the range 0.3κ -0.1κ -shown in Fig. 7 (right) -meaning the effective Reynolds number Re ′ = D(V − V c )/ν ′ was between 0 and ∼ 15.Clearly this was insufficient to support quasi-classical turbulence in the coarse-grained velocity field.Thus, we have polarised ultra-quantum turbulence driven by injections of short-wavelength Kelvin waves, enabled by the frequent reconnection of vortex ends due to the relative motion between the vortex tangle and the rough wall.
The coarse-grained vortex line density profile, shown in Fig. 8a, demonstrates the tangle formed maximally dense regions ∼ D/4 away from each surface, with minima in the centre of the channel and at the walls.Such a spatial variation in L is similar with simulations at finite temperature where flows were driven by a normal fluid flow described by a Poiseuille profile [22].The peaks are mirrored by the polarised line density L y (z) (Fig. 8b) and their ratio Ly(z) L(z) (Fig. 8c) shows the polarisation of the vortex tangle, whose total was between 0 and 60% depending on the distance from walls, in all simulations.

Conclusions
These simulations demonstrate a critical velocity V c = 0.19 ± 0.01 cm s −1 ≃ 20 κ D for pure superfluid flow in channels, above which the turbulent state of the vortex tangle sustained for at least 80 s.There is a finite slip velocity at the solid boundary, whose value is similar to V c , resultant from the "walking" of terminated vortex lines along the rough surface due to frequent reconnections with their images.The coarse-grained velocity profiles u x (z) were similar to classical laminar profiles (albeit with a nonzero velocity at the walls) -consistent with the effective Reynolds number Re ′ ≤ 15 too small to sustain quasi-classical turbulence that typically requires Re > 3000.We thus had ultra-quantum (Vinen's) regime of quantum turbulence.The fraction of the polarised vortex length reached ∼ 60% within the shear regions near both walls, and decreased to zero in the middle of the channel.The effective kinematic viscosity ν ′ , computed for the shear regions, tends towards ∼ 0.1κ with increasing the applied flow velocity -in quantitative agreement with experimental observations for quantum turbulence in superfluid 4 He [1].The flows exhibited an approximately proportional dependence of the friction force F x (V ) and near cubic dependence of the vortex line length Λ(V ) on the mean flow velocity V .These dependences held even below V cfor metastable vortex tangles of lifetime shorter than 80 s.

Fig. 1 (
Fig.1(Left) Visualisation of simulated boundary conditions.The periodicity of the boundaries is shown only in the x direction such that 9 out of 27 cubes are displayed.The original computational volume in the centre is shown by the black box outline.Flat grey surfaces represent the solid boundaries and are extended into the periodic boundaries in the x-direction.The vortex tangle in each cube is coloured to indicate the type of transformation applied: periodic-translation (red), solid-refelction (blue) and a combination of the two (purple).(Right) Illustration of the sphere of velocity contributions centred on a point, s i , near the corner of the volume.Red vortex segments, within the sphere, contribute to ṡi and vortex lines outside of the sphere are neglected.

Fig. 2
Fig.2Progressive snapshots of a single vortex line, which is initially pinned vertically between two flat rough surfaces, evolving under an applied velocity V = 0.05 cm s −1 .Vortex lines drawn in grey are shown at times: 50 ms, 100 ms, 150 ms, 200 ms, 255 ms and 260 ms.The region displayed is the lower 5% of the channel with the solid green line marking the solid boundary.The dashed green line shows the critical distance for vortex-wall reconnection.The solid black lines show the filament in the instant before and after the first time the vortex reconnected with the wall at t = 253.69ms.Coloured line segments (red and blue) show the length removed by the event, and the blue segment is responsible for the reconnection.The size of arrows indicates the spatial resolution, δ and the distance between new and old pinning sites.The inset illustrates the moment of reconnection, with the image-vortex drawn symmetrically behind the boundary.

Fig. 3
Fig. 3 Vortex line length as a function of time shown for 5 different simulations.Solid lines show data and dashed lines show fitted curves.The inset plots are early-time snapshots of the developing vortex tangle for V = 0.30 cm s −1 .

3 Fig. 5 (
Fig. 5 (a) Time-averaged mean flow velocity as a function of the applied superflow shown for both sustained (closed symbols) and unsustained (open symbols) turbulence.(b) The coarse-grained friction force exerted by vortex lines on the channel walls, normalised by a constant vortex line tension, shown to be proportional to the mean flow.(c) Time-averaged vortex line density as a function of the mean flow.Open circles represent unsustainable line densities.

Fig. 6 Fig. 7 (
Fig. 6 Cross-channel coarse-grained velocity profiles.(Left) Absolute values, reduced by the slipvelocity at the walls.(Right) Normalised by the range of velocity present in the flow.