Two-dimensional simulation of vortex points and tracer particles in counterflowing He-II

The article presents results obtained based on numerical simulation of two-dimensional vortex points and inertial tracer particles in two-component counterflowing superfluid He-II. In the low temperature limit (no normal fluid, no friction) our model would reduce to Onsager’s famous vortex gas. The flow of the normal component of the He-II is assumed uniform, while the superfluid velocity field is induced by vortex points which model three-dimensional quantized vortex lines. Probability density functions of velocity and acceleration of tracer particles and superfluid velocity field are obtained. We find that tails of probability distributions follow power-laws with various exponents, except in the case of sufficiently coarse-grained superfluid velocity field, where Gaussian shape is observed. The decay of the number of vortices is also studied, yielding results in agreement with Vinen’s phenomenological model of quantum turbulence.


Introduction
Liquid 4 He cooled below certain temperature (so-called λpoint, approx.2.17 K at SVP) becomes superfluid and is called He-II.Properties of this liquid phase are strongly influenced by quantum-mechanical effects and are conveniently, with sufficient accuracy, described by a two-fluid model.This model identifies two components within He-II -the superfluid component with zero viscosity that carries no entropy, and the viscous normal fluid component that carries all the entropy of the liquid.Due to quantummechanical constraints, the vorticity of the superfluid component cannot be arbitrary, as in classical viscous fluids, but is concentrated to thin topological defects of the size of the order of inter-atomic distance, usually imagined as lines, around which circulation is quantized in units κ ≈ 10 −7 m 2 called quanta of circulation.Turbulence is possible in both components.In the normal component it is assumed that the turbulence is essentially the same as in classical viscous fluids.In the superfluid component, however, the turbulence takes the form of a complex tangle of vortex lines.
Turbulence in He-II has been extensively studied both experimentally and numerically for several decades.The numerical approach usually builds upon Schwarz's model of quantized vortex lines as three-dimensional (3D) space curves [1].Because of the 3D discretization required and the nonlocal character of the vortex-vortex interaction, such calculations are computationally expensive.The two dimensional (2D) model here considered [2] is significantly simpler and idealized, as the space curves representing vortex lines are reduced to vortex points; at zero temperature a e-mail: varga.emil@gmail.com(T = 0) the model would reduce to the famous vortex gas of Onsager.
From the physical point of view, the type of turbulence considered here is peculiar to superfluid 4 He and is known as thermal counterflow.In order to generate thermal counterflow experimentally, a heater is placed at the closed end of a channel.When power is supplied to this heater, the normal component of He-II starts to flow away from the heater, carrying with it all the entropy deposited by it.The superfluid component flows in the opposite direction towards the heater to maintain constant density and zero mass flux (that is, ρ n + ρ s = constant and ρ n v n + ρ s v s = 0).The superfluid and normal fluid velocities along the channel are respectively [4] where Q is the heater's power per unit area (heat flux), ρ s , ρ n and ρ are the superfluid, normal fluid and total densities, respectively, S the specific entropy per unit mass, and v s , v n the superfluid and normal fluid velocities (positive velocities are in the direction away from the heater, negative towards it.)Counterflow turbulence is usually characterized by the vortex line density, L, (defined as the quantized vortex length per unit volume), which is related to the counterflow velocity w = |v n − v s | by Vinen's phenomenological equation [3] where B is the mutual friction coefficient and c 1 , c 2 are parameters of order unity.2) are found by setting w = 0; we obtain: where L 0 is the initial condition at t = 0.The article is structured as follows.The next section introduces the necessary equations used to run the numerical simulations.The third section describes the simulations and the chosen parameters.The fourth section is dedicated to results, and the fifth concludes the paper.

Equations of motion
Let r denote position on the xy plane.Let r j be the location of the vortex point j = 1, • • • N of (positive or negative) circulation κ j .The superfluid velocity at r j is the sum of any applied superflow v ext S and the superflow induced by The vortex points j = 1, • • • N move such that Magnus and friction forces are in balance [1]: where v ext n is the imposed normal fluid velocity and α, α are known temperature-dependent mutual friction coefficients.
In the absence of gravity, the equation of motion of a tracer particle located at r p is [5] where u p = dr p /dt is particle's velocity, τ the viscous relaxation time, ρ p the particle's density, and ρ 0 = ρ p + ρ/2.The assumptions behind equation ( 6) are discussed by Poole et al. [5]; in particular the particle is smaller than any other relevant length scale and does not affect the flow (one-way coupling), the normal fluid's flow around it is laminar (hence Stokes's drag formula applies), and it does not interact with other particles.The substantial derivatives are defined as In the present work the normal fluid velocity is assumed to be uniform (constant in time and space), v ext n = v n and particles are chosen to be neutrally buoyant, ρ p = ρ.With these assumptions equation ( 6) reduces to Fig. 1.A snapshot of vortices and particles in the computational domain.Red "+" symbols represent vortices with positive circulation and blue "x" symbols those of negative circulation.Green circles represent tracer particles.
Assuming that tracer particles are spheres of radius a, the relaxation time τ [5] is τ = 2a 2 ρ 0 /(9μ n ) where μ n is the viscosity of the normal fluid.
Equations ( 5) and ( 8) do not model all the behaviours than needs to be simulated.Firstly, it is known that when two vortex lines get sufficiently close to each other they reconnect [9,10].Since our model does not contain details at this microscopic level (which would require the use of the Gross-Pitaevskii equation), this effect needs to be addressed algorithmically.The reconnection process cannot be ignored because the mutual friction force can push two vortices of opposing signs very close together.In 2D the reconnection algorithm is simple: when two vortices of opposing signs approach each other (closer than an arbitrary reconnection length, see table 1), they annihilate, and, unless we study the decay of turbulence, they are re-inserted back into the computational domain at random positions to maintain the statistical steady state.
Secondly, tracer particles can become trapped on vortices.To take particle trapping into account we implement the following algorithm: when a particle becomes closer to a vortex than an arbitrary trapping length (see table 1), we tag it as trapped and let it move as the vortex hereafter.We allow a trapped particle to become free (detrapping) if the vortex to which it was trapped undergoes a reconnection with a vortex of the opposite sign.Reducing the reconnection distance by a factor of 100 did not alter the conclusion.

Numerical calculations
The computational domain is a square box of size D with periodic boundary conditions, implemented by repeating the contents of the computational domain at each boundary once, i.e., north, south, east and west.Equations ( 5) and ( 8) are then solved in the dimensionless form.The choice of units is such that (i) the side of the square computational domain is 1, and (ii) Once the spatial scaling is chosen to satisfy (i), condition (ii) defines the time scale uniquely.The size of the square box is set to fix the mean inter-vortex distance = D/ √ N (where N is the number of vortices) to experimentally known values taken from reference [6].
All calculations used N = 300 vortices (150 positive and 150 negative) and N p = 200 tracer particles.The vortices were initially (t = 0) set at random locations.The particles's initial position was also random, with velocity matching v n .The system was evolved for one (dimensionless) time unit, sampling positions, velocities and accelerations of tracer particles periodically.Data sampled this way were treated as a single ensemble for calculating statistics, so histograms of velocity or accelerations which we present represent time averages.In order to sample the superfluid velocity field we set up a regular grid of 1000 × 1000 points where we evaluated velocity and Lagrangian acceleration, that is, the substantial derivative of the superfluid velocity field defined by equation (7).Data were sampled 10 times during the calculation, and histograms were calculated from the entire ensemble.All helium fluid parameters (densities, specific entropy and mutual friction coefficients) were taken from Reference [11].The parameters which we used are listed in table 1.An adaptive 4th order Runge-Kutta scheme was used for integrating the equations of motion in time.
A illustrative snapshot of the positions of vortices and particles in the computational domain is shown in figure 1.

Results
The most important statistical properties are the mean values of the magnitudes of velocity and acceleration over the time evolution; these quantities are shown in figures 2 and 3 as functions of the imposed heat flux.
It is apparent from figures 2 that the free particles' velocity closely follows that of the counterflowing normal fluid, while the velocity of trapped particles systematically exceeds that of the counterflowing superfluid.By construction, the latter is also the velocity of the vortices.This difference is probably due to the first mutual friction term in equation (5), which imparts to the vortices (and hence to the trapped particles) a motion in the direction transverse to the direction of the imposed counterflow.The result that the mean velocity of all particles is close to the mean velocity of trapped particles is probably due to the fact that, as a detrapping mechanism, 2D vortex reconnections is not effective enough -at any instant after an initial transient most particles in the simulations are trapped.
Mean acceleration magnitudes are shown in figure 3, averaged over all particles or separated as for velocities.This set of results allows direct comparison with experiments, for example, this quantity was measured by La Mantia et al [6].Experimental values from that work are included in the figure.The agreement with the experiment is rather poor, within an order of magnitude.A possible explanation for this discrepancy (besides the fact that our model is only 2D) is that experimental values are inferred from observed trajectories, whereas in our simulations, the instantaneous acceleration is calculated from the known positions and velocities of the vortices.Such quantity will be subject to strong fluctuations which, however, do not cancel out in the mean when calculating acceleration magnitudes.
Our data also allow for calculation of normalized histograms, or probability density functions (PDFs).PDFs of Fig. 2. Mean velocity magnitudes of tracer particles calculated at T = 1.75 K.Here the particles are separated into two groups: free and trapped particles.Mean values are taken over the two groups and over the entire ensemble as well.Black squares denote magnitudes of velocity, averaged over all particles, green triangles over trapped particles only, and red circles over free particles only.The upper black line is the imposed normal fluid velocity, the lower red line is the imposed superfluid velocity; both lines are calculated from equation ( 1). .Fig. 3. Mean magnitudes of accelerations of tracer particles.The particles are separated into groups as in figure 2. Green upwardspointing triangles correspond to free particles, blue downwards pointing triangles to trapped particles, and red circles to average over all particles.Black squares represent experimental values from [6].Temperature is 1.75 K.
transverse (to the mean direction of counterflow) component of velocity and acceleration for tracer and fluid particles are in figures 4 -7.In these plots, the x axes are normalized by the standard deviation: numerical values are in table 2. For tracer particles' PDFs, no distinction between trapped or free particles was made.In all four cases, power-law tails are observed, as clearly visible in the insets of these figures.
For velocities of fluid particles, the obtained exponent of the power-law tail is consistent with the expected value of −3 [7].It is worth noting that the inertial tracer particles are not found to faithfully follow the statistics of the fluid particles.However, it is possible that small size statistics for the tracer particles could be the reason.
The tracer particles' PDFs do not overlap as closely as the PDFs of the fluid particles.However, the apparent order appears to be random and is not correlated with the power to the heater.

02124-p.3
Table 1.Parameters used in numerical simulations.The meaning of columns is as follows: at given temperature T , D denotes the size of the computational box, T t and U = D/T t are the units of length, time and velocity, respectively; Q denotes the heat flux supplied to the counterflow heater; v n and v s stand for the counterflow-imposed velocities of the normal fluid and the superfluid; τ is the viscous relaxation time.In all cases, particle radius was set to be 5μm, reconnection length D/500 and trapping length D/500 + (particle radius).Fig. 6.Probability distribution of the superfluid velocity component perpendicular to the mean counterflow, normalized by the standard deviation (see table 2).The figure contains 4 curves with colors and styles as in figure 4. The inset shows comparison with the predicted power law with the −3 exponent, in log-log plot.
Fig. 7. Probability distribution of superfluid Lagrangian particle acceleration component perpendicular to the mean counterflow, normalized by the standard deviation (see table 2).The figure contains 4 curves as figure 6.The inset shows the same data on a log-log scale with a power law of fitted exponent −1.43.
The observed sharp cutoff in the fluid particles' accelerations, at around 9 standard deviations, is probably due to the fact that points where the acceleration is sampled do not get arbitrarily close to vortices.
The superfluid velocity field (velocities of the fluid particles) are sampled on a regular square grid.This allows for coarse-graining of the velocity field, thus modelling the finite-size of a measurement region.Figure 8, with numerical values for standard deviations given in table 3, shows the progression of change of the PDF as the grid gets more and more coarse.The power-law tails disappear Fig. 8. Probability distribution of coarse-grained superfluid velocity.The legend in the figure corresponds to the curves in order from top to bottom.The coarse-graining is accomplished by imposing a coarser mesh on underlying mesh of 1000×1000 points where velocity is calculated and averaging the velocity in the cells of the coarser mesh.The last curve, 17 × 17, approximately corresponds to coarse-graining on the scale of inter-vortex distance.The numbers in the legend show how many cells are in the coarse mesh.The temperature is 1.75 K and heater power is 608 W. and the distributions assumes an approximately Gaussian form.Similar behaviour is observed in [12].
It is also possible to model the decay of turbulence by not reinserting back the vortices that annihilate.Figure 9 shows the decay starting from three different numbers of vortices at 1.65 K. Figure is kept in dimensionless time, because to fix the time scale vortex line density must be known.This can be inferred from the intended heater power, but since during the decay the heater is switched off, the choice would be arbitrary.All curves, however, have the same time scale.The curves are found to be in good agreement with Vinen's decay (3).

Conclusions
We performed computer simulations of vortex points and tracer particles in a simple two-dimensional model of superfluid helium.We found that the mean particle velocity closely follows the velocity of either component of counterflowing He-II, depending on whether particles are trapped into vortices or not.The tracer particle accelerations are much higher than those observed in the experiments; this effect is probably caused by the 2D trapping and detrapping mechanisms which we used, and will require further attention.
The probability distribution functions of velocity and acceleration components of tracer particles and fluid parti-   2).Line color and style have the same meanings as in figure 4. Inset shows the same data in log-log coordinates, highlighting the power-law tails with fitted exponent −1.78.
T (K) Q (W/m 2 ) std(v x ) (mm/s) std(a x ) (mm/s 2 ) grid, std(v x ) (mm/s) grid, std(a x ) (mm/s  cles are found to have power-law tails.The power-law exponent of the fluid particle velocity is consistent with the expected value of −3.Fluid particles and tracer particles have different exponents.Finally, we found that the decay of the number of vortex points obey Vinen's prediction, indicating that the movement of vortex points is random.

Fig. 4 .
Fig. 4. Probability density function of tracer particle velocity component perpendicular to the mean counterflow, normalized by the standard deviation (see table 2). Figure contains 4 curves, corresponding to 4 heat flux values (or mean counterflow velocities) as in figures 2 to 3, in W/m 2 : 414 (solid black curve), 590 (dashed red), 608 (dotted green) and 691 (dot-dashed blue).The inset shows the same data in log-log scale, highlighting the power-law tails with fitted exponent -4.3.The exponent of the power law tail is −4.3 ± 0.3.Temperature is 1.75 K.

Fig. 5 .
Fig.5.Probability distribution of tracer particle acceleration perpendicular to the mean counterflow normalized by the standard deviation (see table2).Line color and style have the same meanings as in figure4.Inset shows the same data in log-log coordinates, highlighting the power-law tails with fitted exponent −1.78.

Fig. 9 .
Fig. 9. Decay of the number of vortex points, N, starting from three initial numbers of vortices N 0 : 8192, 4096 and 2048.The curves corresponding to N 0 of 4096 and 2048 are the result of ensemble averaging of 10 individual calculations.The scaling of the dimensionless time with respect to the real time is the same for all curves and time or axes of the curves are not manipulated in any way: curves naturally collapse onto each other.

Table 3 .
Standard deviations for the curves in figure8.