# Evolution of a superfluid vortex filament tangle driven by the Gross-Pitaevskii equation

Alberto Villois, Davide Proment, Giorgio Krstulovic

► To cite this version:

Alberto Villois, Davide Proment, Giorgio Krstulovic. Evolution of a superfluid vortex filament tangle driven by the Gross-Pitaevskii equation. Physical Review E, 2016, 93 (6), pp.061103. 10.1103/Phys-RevE.93.061103. hal-03526714

HAL Id: hal-03526714

https://hal.science/hal-03526714

Submitted on 17 Jan 2022

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
PHYSICAL REVIEW E 93, 061103(R) (2016)

# Evolution of a superfluid vortex filament tangle driven by the Gross-Pitaevskii equation

Alberto Villois and Davide Proment

School of Mathematics, University of East Anglia, Norwich Research Park, Norwich NR4 7TJ, United Kingdom

Giorgio Krstulovic

Laboratoire J.L. Lagrange, UMR7293, Université de la Côte d'Azur, CNRS, Observatoire de la Côte d'Azur,
Boîte Postale 4229, 06304 Nice Cedex 4, France

(Received 2 May 2016; published 30 June 2016)

The development and decay of a turbulent vortex tangle driven by the Gross-Pitaevskii equation is studied. Using a recently developed accurate and robust tracking algorithm, all quantized vortices are extracted from the fields. The Vinen's decay law for the total vortex length with a coefficient that is in quantitative agreement with the values measured in helium II is observed. The topology of the tangle is then investigated showing that linked rings may appear during the evolution. The tracking also allows for determining the statistics of small-scale quantities of vortex lines, exhibiting large fluctuations of curvature and torsion. Finally, the temporal evolution of the Kelvin wave spectrum is obtained providing evidence of the development of a weak-wave turbulence cascade.

DOI: 10.1103/PhysRevE.93.061103

The full understanding of turbulence in a fluid is one of the oldest yet still unsolved problems in physics. A fluid is said to be turbulent when it manifests excitations occurring at several length scales. Due to the large number of degrees of freedom and the nonlinearity of the governing equations of motion, the problem is usually tackled statistically by introducing assumptions and closures in terms of correlators. This is the case in the seminal work of Kolmogorov in 1941 based on the idea of Richardson's energy cascade, where energy in classical fluids is transferred from large to small scales [1].

Superfluids form a particular class among fluids characterized essentially by two main ingredients: the lack of dissipation and the evidence that vortex circulation takes only discrete values that are multiples of the quantum of circulation [2]. Superfluid examples are superfluid liquid helium (He II) and Bose-Einstein condensates (BECs) made of dilute alkali-metal gases. Here the superfluid phase is usually modeled via a complex field describing the order parameter of the system and vortices appear as topological defects where the superfluid density vanishes.

In three spatial dimensions those defects organize themselves into closed lines (or even open lines at the boundaries if confining sides are considered) of different configurations. Any vortex line point induces a velocity field which affects the motion of any object in the system including the vortex line itself. In general, even for a single closed vortex line, the dynamics are chaotic and the problem does not have analytical solutions. Superfluid turbulence regards the study of the evolution of many vortex lines, a tangle, which induce velocity field gradients in the fluid at several length scales.

Different mathematical models have been devised to mimic the dynamics of a superfluid. An example is the vortex filament (VF) model based on the Biot-Savart law that relates vorticity and velocity [3]. This model is able to mimic the dynamics of dense vortex tangles due to a relatively fast numerical integration technique [4]. The VF model implicitly assumes that the superfluid density is constant everywhere and the vortex structure is a line with vanishing core. This assumption is generally satisfied in He II where the characteristic experimental setup sizes, and consequently the largest scales

of the motion, are order of $10^{-1}$ m and the vortex core is of the order of $1 \AA = 10^{-10}$ m. Moreover, since He II is in its liquid phase, the compressibility effects can be usually neglected. However, the VF model fails to describe vortex reconnections. These are rapid changes in the topology of the vortex configuration which occur naturally in a superfluid [5] and are one of the main mechanisms responsible for the energy transfer. Reconnections are thus introduced by some *ad hoc* mechanism.

Another superfluid model that admits quantized vortices and inherently possesses vortex reconnections is the Gross-Pitaevskii (GP) equation that describes the evolution of the superfluid order parameter $\psi$. In contrast to the VF model, the GP equation allows density fluctuations in terms of phonons and density depletion at the vortex cores. Although it has been formally derived as a mean-field theory for a dilute boson gas in the limit of zero temperature [6], it also qualitatively reproduces He II dynamics. The vortex core size here is of the order of the healing length $\xi$, the only intrinsic characteristic length scale of the model; nowadays experimental techniques are able to create BEC setups that are $10^1-10^2$ healing lengths where superfluid turbulence can develop [7]. In turbulent superfluids, vortices constantly rearrange themselves following reconnections into complex tangles with nontrivial geometrical, algebraic, and topological properties [8]. At small scales, helical excitations of vortex lines known as Kelvin waves (KWs) are believed to be the ultimate mechanism of energy dissipation via phonon emission [9]. To study such dynamics, the GP equation has the advantage that no extra modeling is needed (unlike the VF model). However, GP does not provide direct information on vortices.

In this work we apply a numerical algorithm [10] to accurately track the configuration of a turbulent vortex tangle evolving according to the GP model. First, we show that after the onset of turbulence, the vortex line density satisfies the Vinen's decay law [11] with a coefficient that is in agreement with the values measured in He II. Different algebraic and topological quantities of the tangle are then measured. The tracking allows for obtaining curvature and torsion distributions of the vortex tangle. Finally, we perform a direct measurement of
KWs during the dynamics and compute a KW spectrum that appears to be consistent with L'vov-Nazarenko's weak-wave turbulence theoretical prediction [12].

The GP model for the condensate wave function $\psi$ is

$$i\hbar \frac{\partial \psi}{\partial t} = -\frac{\hbar^2}{2m} \nabla^2 \psi + g|\psi|^2 \psi, \quad (1)$$

where $m$ is the mass of the bosons and $g = 4\pi a\hbar^2/m$, with $a$ the s-wave scattering length. Madelung's transformation $\psi(\mathbf{x},t) = \sqrt{\rho(\mathbf{x},t)/m} \exp[i\frac{m}{\hbar}\phi(\mathbf{x},t)]$ relates the wave function $\psi$ to a superfluid of density $\rho(\mathbf{x},t)$ and velocity $v = \nabla\phi$. The quantum of circulation about the $\psi = 0$ vortex lines is $\Gamma = h/m$. When Eq. (1) is linearized about a constant value $\psi = \hat{\psi}_0$, the sound velocity is given by $c = (g|\hat{\psi}_0|^2/m)^{1/2}$, with dispersive effects taking place at length scales smaller than the healing length $\xi = (\hbar^2/2m|\hat{\psi}_0|^2g)^{1/2}$.

In the simulations presented here, the mean density is fixed to the unity and the physical constants in Eq. (1) are determined by the values of $\xi$ and $c = 1$. The quantum of circulation results in $\Gamma = 4\pi c \xi / \sqrt{2}$. Numerical integration of Eq. (1) is performed using a standard pseudospectral code. We integrate an initial condition characterized by the so-called Taylor-Green flow [13], a well-studied flow in superfluid turbulence. Symmetries are not enforced during the evolution and we use resolutions of 256³ and 512³ uniformly distributed collocation points with $\xi = 2\pi/256$ and $\xi = 2\pi/512$, respectively. Mirror symmetries are broken during the evolution although traces of such symmetries will be present even at very large times. With units used in this work, the large eddy turnover time is of the order of the unity.

The Taylor-Green flow initially contains a configuration of unstable large-scale rings that develop to create a turbulent tangle. Vortices can be spotted by plotting the low-value isosurfaces of the density field as displayed in Fig. 1. Low-density regions corresponding to vortex lines are plotted in red, while density fluctuations (sound) are rendered in light blue. The initial condition is visualized in Fig. 1(a), the complex turbulent tangle at $t = 12$ in Fig. 1(c), and the final state at $t = 105$, where few vortices are present with a lot of sound in the background, in Fig. 1(e). We track the vortex lines with a recently developed algorithm [10] that allows for identifying separately each single line forming the tangle. Vortex lines are followed using the pseudo-vorticity field as in [14] and the exact vortex positions are obtained by applying a Newton-Raphson method. The algorithm is robust and accurate as it takes full advantage of the spectral resolution. The intermesh values of the field $\psi$ and its derivatives needed for the Newton-Raphson method are directly evaluated by Fourier transforms; the locations of vortices are thus found with precision given the machine-$\epsilon$ (double in the present simulations). See [10] for all technical details and a complete validation of the algorithm. Figures 1(b), 1(d), and 1(f) show the corresponding tracked vortices displayed in different colors (see Supplemental Material [15]).

We focus first on the later evolution times. During the decay, vortices radiate phonons at small scales creating a thermal bath that exchanges energy and momentum with the vortices. This process mimics mutual friction and leads eventually to the total annihilation of vortex rings [16]. In superfluids such a

<figure><img src="image_2.png" /><figcaption>FIG. 1. Left: Isosurfaces of density field at different times. Low-density regions that correspond to vortex lines are plotted in red, while density fluctuations (sound) are rendered in light blue. Right: Corresponding tracked vortices. Different colors correspond to different vortices. Snapshots taken at $t = 0$ [(a) and (b)], $t = 12$ [(c) and (d)], and $t = 105$ [(e) and (f)]. Resolution 256³.</figcaption></figure>

decay is modeled by the Vinen equation [11] for vortex line density $\mathcal{L}$:

$$\frac{d\mathcal{L}}{dt} = -\chi_2 \frac{\Gamma}{2\pi} \mathcal{L}^2, \quad (2)$$

where $\chi_2$ is a constant of the order of unity. Its solutions manifest a $\mathcal{L} \sim t^{-1}$ behavior at long times: this power-law decay has been named quantum turbulent decay and measured in He II experiments [17] and VF numerical simulations [18]. In Fig. 2(a) we show the temporal evolution of $\mathcal{L}$. It is worth noticing that it grows at the initial stages: this is caused by the instability of the initial Taylor-Green configuration and the subsequent vortex stretching due to numerous vortex reconnections. The data is compared with an estimation of $\mathcal{L}$ obtained by computing the ratio between the volume of points such that $\rho(\mathbf{x}) < 0.2$ and the corresponding surface of a perfect two-dimensional vortex profile. This latter method has become a standard technique within GP numerical simulations to compute the vortex line density [19]. Even if this technique is able to capture the qualitative behavior of $\mathcal{L}$, it fails to grasp at long times the power law predicted by Vinen's
<figure><img src="image_3.png" /><figcaption>FIG. 2. (a) Temporal evolution of the vortex line density: tracked data are plotted using blue circles, and volume estimation in a solid red line. (b) Long time decay of $\Delta\mathcal{L}$ (see text), together with Vinen's prediction $\Delta\mathcal{L}_{\text{Vinen}} = [\chi_2 \frac{\Gamma}{2\pi}(t-t_0)]^{-1}$, setting $\chi = 0.65$ (solid black line). Resolution 256³.</figcaption></figure>

equation. This is shown in Fig. 2(b) where the measured $\Delta\mathcal{L}(t) = [\mathcal{L}(t)^{-1} - \mathcal{L}(t_0)^{-1}]^{-1}$, setting $t_0 = 17$, is compared to Vinen's prediction. We can explain this discrepancy by reasoning that the vortex core size (proportional to the uniform condensate state) varies in time because more and more sound excitations are created by the superfluid decay, altering the estimation of $\mathcal{L}$ by fixing the (non-time-dependant) density threshold. The tracked data also allow for determining the numerical constant $\chi_2 = 0.65$. This value is in remarkable agreement with experimental values measured in He II in the low temperature limit [11]. Between the time of maximal vortex length ($t \approx 8.5$) and $t_0$ there is a faster decay that could be explained by the quasiclassical turbulent decay law [20], although the data (not shown here) do not allow for a precise corroboration and further studies are needed.

From Fig. 1 and the movie provided in the Supplemental Material, it is clear that the complexity of tangle first increases and then decreases. The complexity of tangle can be measured by computing the changes in some of its algebraic and topological quantities [8]. We compute the total average crossing $\bar{C} = \sum_{i \neq j} C_{i,j}$, the total linking $L_k = \sum_{i \neq j} L_{ij}$, and the writhe $W_r = \sum_i W_{ri}$, by directly performing the line integrals over the vortex ring(s) [8] as

$$ \bar{C}_{i,j} = \frac{1}{4\pi} \oint_{C_i} \oint_{C_j} \left| \frac{(\mathbf{R}_i - \mathbf{R}_j) \cdot d\mathbf{R}_i \times d\mathbf{R}_j}{|\mathbf{R}_i - \mathbf{R}_j|^3} \right|, \quad (3) $$

$$ L_{k,j} = \frac{1}{4\pi} \oint_{C_i} \oint_{C_j} \frac{(\mathbf{R}_i - \mathbf{R}_j) \cdot d\mathbf{R}_i \times d\mathbf{R}_j}{|\mathbf{R}_i - \mathbf{R}_j|^3}, \quad (4) $$

$$ W_{r_i} = \frac{1}{4\pi} \oint_{C_i} \oint_{C_i} \frac{(\mathbf{R}_i - \mathbf{R}_i') \cdot d\mathbf{R}_i \times d\mathbf{R}_i'}{|\mathbf{R}_i - \mathbf{R}_j|^3}. \quad (5) $$

<figure><img src="image_4.png" /><figcaption>FIG. 3. (a) Temporal evolution of the (normalized) total number of rings and crossing number. At $t = 0$, $N_{\text{rings}}(0) = 128$ and $\bar{C}(0) = 758$. (b) Temporal evolution of the total linking $L_k$, writhe $W_r$ and center-line helicity $H_c/\Gamma^2$. (c) Visualization of two linked rings at $t = 21$. (d) Visualization of a ring with high $W_r$ at $t = 24.5$. Resolution 256³.</figcaption></figure>

Here $\mathbf{R}_i$ corresponds to the points identifying the $i$th ring $C_i$; for the writhe number, $\mathbf{R}_i$ and $\mathbf{R}_i'$ correspond to two different points of the same ring. In Fig. 3(a) we plot the total number of rings $N_{\text{rings}}$ and $\bar{C}$ normalized by their initial values versus time. It is worth noticing that the average crossing number reaches qualitatively a maximum at the same stage of the vortex line $\mathcal{L}$ maximum, while the ring number maximum is slightly shifted forward in time. The former observation follows the idea that vortex lines simultaneously stretch, bend, and coil during reconnection events. The latter is due to the fact that longer vortex rings continue to break into pieces during the evolution until the tangle density becomes low enough and the main vortex length dissipation mechanism is given by sound interaction. We then focus on the center-line helicity $H_c/\Gamma^2 = L_k + W_r$ [21] related to the helicity in classical fluid dynamics, an important inviscid invariant. The linking number $L_k$ takes integer values and gives information about the number of linked rings present in the system, whereas the writhe takes real values and its contribution
comes from self-linked (knots), an integer contribution, and
KWs [22]. Figure 3(b) shows the temporal evolution of
these three quantities. Initially, Lk = Wr = 0, as expected for
the Taylor-Green flow. Surprisingly, during the evolution Lk
becomes nonzero, indicating the presence of linked rings, such
as the ones displayed in Fig. 3(c) [23]. This is remarkable
as in the GP model sufficiently simple vortex configurations
usually decay by reducing their complexity [24]. Once the
decay is established, no linked rings are present and only
writhe contributes to $H_c$. Note that the writhe number is not
enough to determine whether a ring is self-linked (knotted)
or not. The center-line helicity, however, fluctuates about a
zero mean, an indication of the presence of KWs. KWs are
indeed apparent in Fig. 3(d) where an unknotted ring with high
Wr is displayed. KWs have already been indirectly observed
in the Taylor-Green flow during the turbulent stage [25], in
agreement with the large values of writhe observed around
t ~ 10.

We now study statistical properties of some geometrical
quantities of the vortex filaments by exploring the time
behavior of the probability density functions (PDFs) of the
curvature κ and torsion τ of the entire set of vortices in
the system. In Fig. 4(a) we present the PDF of curvature,
normalized by its mean value, at different stages. The temporal
evolution of the mean curvature ⟨κ⟩ and its rms value κ<sub>rms</sub> are
also displayed in the inset. We can observe that ⟨κ⟩ increases
rapidly at early stages and then almost saturates, an indication
that the average vortex size (inversely proportional to the
curvature) slowly decreases at later times. The rms value of
the curvature presents the same tendency with the exception
of peaks. These are evidence of reconnection events where
high values of curvature are found in localized regions. It is
worth noticing that the PDFs, rescaled by their mean curvature,
exhibit a relatively good collapse to a self-similar form. This

<figure><img src="image_5.png" /><figcaption>FIG. 4. (a) PDFs of curvature κ normalized by their respective mean values ⟨κ⟩ at different times [same legend as (b)]. The inset displays the temporal evolution of the mean and rms values of κ. (b) PDFs of torsion τ at different times. The inset emphasizes their τ⁻³ power-law tail. Resolution 256³.</figcaption></figure>

latter observation indicates a power-law behavior ~κ¹ at
small curvature values, while an exponentially decaying tail
is present at large curvature values. A similar behavior has
also been observed within the VF model [26]. In Fig. 4(b) we
plot the torsion PDFs at the same stages. The mean torsion is
always about zero and there is no evidence of any skewness in
the PDFs. The distributions' tails show a universal power-law
behavior of τ⁻³ at all times, meaning that the second and
higher moments of the torsion diverge during the decay. The
same scaling appears in vortex tangles of random wave fields
that are solutions of the Helmholtz equation [27]. This may be
an indication that for one-time small-scale quantities, quantum
turbulent tangles can be interpreted simply as random vortices.

The large curvature fluctuations and the torsion fluctuation
about a zero mean are evidence of KWs at all scales
propagating on quasiplanar vortex rings. By exploiting the
accuracy of the tracking algorithm we are able to directly
detect KWs on those rings. Competing theories have been put
forward to statistically predict a power-law KW spectrum in
the form of $n_k \sim k^{-\alpha}$ (here k is the Kelvin wave number)
and explain the energy transfer through KW scales. Vinen
et al. considered strong nonlinear interactions and derived by a
scaling argument the exponent $\alpha_V = 3$ [28]. On the other hand,
assuming weak nonlinearity (small amplitude KWs compared
to their respective wavelengths), Kozik and Svistunov [29] and
L'vov and Nazarenko [12] obtained the exponents $\alpha_{KS} = 17/5$
and $\alpha_{LN} = 11/3$, respectively, considering two different orders
of interaction. We can compute the KW spectrum of a ring R by
applying a Gaussian kernel of width $\alpha L$ in order to establish
the configuration of the unperturbed ring $R_{ump}$. This can be
used to define the KWs on it as $R_{KW}(s) = R(s) - R_{ump}(s)$,
where $s \in [0, L]$ is the arc-length parametrization of the ring.
Being $R_{KW}$ a periodic set of three signals (one for each
spatial dimension), the KW spectrum is then defined as
$n_k = |\tilde{R}_{KW}(k)|^2 + |\tilde{R}_{KW}(-k)|^2$, where $\tilde{R}_{KW}(k)$ is the Fourier
transform of $R_{KW}(s)$. In [10] we checked that this procedure
is able to capture well the KWs superimposed on a ring. Here
we compute the KW spectrum averaging over the spectra of
the 50 largest rings such that it has small fluctuations and it
always spans over two Kelvin wavelength decades. For the
Gaussian filter, we use the value $\alpha = 0.1$; varying this fraction
weakly modifies the large-scale values of the spectrum, but
the data in the inertial range remain unchanged. The KW
spectra are shown for different times in Fig. 5(a). It is evident
that all accessible KW modes get populated at early times
due to reconnection events that trigger the cascade [30]. We
observe KW spectra exhibiting power laws with an exponent
independent of time where the best scaling is appreciated at
the time where the rings are the longest ($4 \le t \le 7$). To get
the best estimation of the power-law exponent, we repeated
the Taylor-Green decay in a simulation box twice larger; in
this new configuration the scaling range spans almost two
wave-number decades. In Fig. 5(b) we show the spectrum
at $t \sim 5$: the observed power-law exponent is close to the
weak-wave turbulence predictions and seems to agree with
the L'vov and Nazarenko $\alpha_{LN} = 11/3$ one. This can be better
appreciated by looking at the compensated spectra with respect
to $\alpha_{LN}$ and $\alpha_{KS}$ shown in the inset. This finding supports the
result in favor of L'vov and Nazarenko's prediction previously
obtained while studying the KW oscillations about a perfect
EVOLUTION OF A SUPERFLUID VORTEX FILAMENT ...

PHYSICAL REVIEW E 93, 061103(R) (2016)

<figure><img src="image_6.png" /><figcaption>FIG. 5. (a) Temporal evolution of KW spectra (averaged over the 50 longest rings). Resolution 256³. (b) KW spectrum at $t \sim 5$ (averaged over the 50 longest rings) for run at resolution 512³. The dashed line displays the $k^{-11/3}$ scaling. The inset displays the respective $k^{11/3}$ (solid blue) and $k^{17/5}$ (dashed red) compensated spectra.</figcaption></figure>

straight line in the GP model [31]. We highlight that although the weak-wave turbulence prediction for the KW spectrum is formally derived for KWs on an isolated straight vortex line using the VF model, it remarkably turns out to be valid in a dense turbulent tangle also driven by the GP model. This

is certainly due to the fact that the predicted KW spectrum was found for the longest rings. Small rings quickly lose their energy by phonon radiation and exchange momentum with sound waves. Both contributions are important to understand dissipation of superfluids at very low temperature and further studies are still needed to fully comprehend the relevance of such mechanisms.

Tracking vortices in GP turbulence opens up a new way for studying and understanding the topological configuration and properties of quantum vortex tangles. Although unlikely, we show that rings can link creating a local (in time and space) fluctuation of the center-line helicity. It will be of great interest to repeat a similar analysis setting where the mean helicity of the flow is not zero, like the ABC flow introduced in [22] where linking and self-linking processes could be substantially enhanced. Overall, the results presented in this work confirm that some predictions traditionally associated to superfluid liquid helium become important in weakly interacting BECs at low temperature described by the GP model. Nowadays BEC experimentalists are able to create and track few vortices in harmonic traps [32,33]. A controlled experimental setting with a turbulent BEC, such as the one presented in this work, has yet to be achieved but it should be realizable in the near future.

G.K., D.P., and A.V. were supported by the Royal Society and the Centre National de la Recherche Scientifique (CNRS) through the International Exchanges Cost Share Scheme (Ref. IE150527). Computations were carried out on the Mésocentre SIGAMM hosted at the Observatoire de la Côte d'Azur and on the High Performance Computing Cluster supported by the Research and Specialist Computing Support service at the University of East Anglia. Isosurfaces in Fig. 1 and in the movie provided in the Supplemental Material have been rendered using VAPOR [34].

[1] U. Frisch, *Turbulence: The Legacy of AN Kolmogorov* (Cambridge University Press, Cambridge, UK, 1995).

[2] R. Donnelly, *Quantized Vortices in Helium II* (Cambridge University Press, Cambridge, UK, 1991), Vol. 3.

[3] K. W. Schwarz, *Phys. Rev. B* **38**, 2398 (1988).

[4] A. W. Baggaley and C. F. Barenghi, *J. Low Temp. Phys.* **166**, 3 (2011).

[5] G. P. Bewley, M. S. Paoletti, K. R. Sreenivasan, and D. P. Lathrop, *Proc. Natl. Acad. Sci. USA* **105**, 13707 (2008).

[6] L. Pitaevskii and S. Stringari, *Bose-Einstein Condensation* (Oxford University Press, New York, 2003), Vol. 116.

[7] E. Henn, J. Seman, G. Roati, K. Magalhaes, and V. Bagnato, *Phys. Rev. Lett.* **103**, 045301 (2009).

[8] C. F. Barenghi, R. L. Ricca, and D. C. Samuels, *Phys. D* (Amsterdam, Neth.) **157**, 197 (2001).

[9] W. F. Vinen, *Phys. Rev. B* **64**, 134520 (2001).

[10] A. Villois, D. Proment, H. Salman, and G. Krstulovic, *arXiv:1604.03595*.

[11] W. F. Vinen, *Proc. R. Soc. London, Ser. A* **242**, 493 (1957).

[12] V. L'vov and S. Nazarenko, *JETP Lett.* **91**, 428 (2010).

[13] C. Nore, M. Abid, and M. Brachet, *Phys. Fluids* **9**, 2644 (1997).

[14] C. Rorai, J. Skipper, R. M. Kerr, and K. R. Sreenivasan, *arXiv:1410.1259*.

[15] See Supplemental Material at http://link.aps.org/supplemental/10.1103/PhysRevE.93.061103 for a movie of full time evolution.

[16] G. Krstulovic and M. Brachet, *Phys. Rev. E* **83**, 066311 (2011).

[17] P. M. Walmsley and A. I. Golov, *Phys. Rev. Lett.* **100**, 245301 (2008).

[18] A. W. Baggaley, C. F. Barenghi, and Y. A. Sergeev, *Phys. Rev. B* **85**, 060501 (2012).

[19] A. C. White, C. F. Barenghi, N. P. Proukakis, A. J. Youd, and D. H. Wacks, *Phys. Rev. Lett.* **104**, 075301 (2010).

[20] S. R. Stalp, L. Skrbek, and R. J. Donnelly, *Phys. Rev. Lett.* **82**, 4831 (1999).

[21] M. W. Scheeler, D. Kleckner, D. Proment, G. L. Kindlmann, and W. T. M. Irvine, *Proc. Natl. Acad. Sci. USA* **111**, 15350 (2014).

[22] P. Clark di Leoni, P. D. Mininni, and M. E. Brachet, *arXiv:1602.06880*.
[23] Although linked rings are present in the flow, the probability of finding them is very small. It can be estimated counting all the linked rings and it is at most of the order $10^{-4}$ for the Taylor-Green flow.

[24] D. Kleckner, L. H. Kauffman, and W. T. M. Irvine, Nat. Phys. (2016), doi:10.1038/nphys3679.

[25] P. Clark di Leoni, P. D. Mininni, and M. E. Brachet, Phys. Rev. A 92, 063632 (2015).

[26] L. Kondaurova, V. L'vov, A. Pomyalov, and I. Procaccia, Phys. Rev. B 89, 014502 (2014).

[27] A. J. Taylor and M. R. Dennis, J. Phys. A: Math. Theor. 47, 465101 (2014).

[28] W. F. Vinen, M. Tsubota, and A. Mitani, Phys. Rev. Lett. 91, 135301 (2003).

[29] E. Kozik and B. Svistunov, Phys. Rev. Lett. 92, 035301 (2004).

[30] D. Kivotides, J. C. Vassilicos, D. C. Samuels, and C. F. Barenghi, Phys. Rev. Lett. 86, 3080 (2001).

[31] G. Krstulovic, Phys. Rev. E 86, 055301 (2012).

[32] G. Lamporesi, S. Donadello, S. Serafini, F. Dalfovo, and G. Ferrari, Nat. Phys. 9, 656 (2013).

[33] S. Serafini, M. Barbiero, M. Debortoli, S. Donadello, F. Larcher, F. Dalfovo, G. Lamporesi, and G. Ferrari, Phys. Rev. Lett. 115, 170402 (2015).

[34] See https://www.vapor.ucar.edu.