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

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 quantised 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 studied showing that linked rings may appear during the decay. The tracking also allows for determining the statistics of small-scales 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.

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 quantised 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 studied showing that linked rings may appear during the decay.The tracking also allows for determining the statistics of small-scales 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.
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 characterised essentially by two main ingredients: the lack of dissipation and the evidence that vortex circulation takes only discrete values multiple of the quantum of circulation [2].Superfluid examples which are routinely created in laboratories are superfluid liquid Helium (He II) and Bose-Einstein condensates (BECs) made of dilute Alkali gases.Here the superfluid phase is usually modelled via a complex field describing the order parameter of the system and quantised vortices appear as topological defects where the superfluid density vanishes.
In three spatial dimensions those defects organise 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 in the superfluid 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 field is constant everywhere and the vortex structure is precisely 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 order of 1 Å = 10 −10 m.Moreover, since Helium 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 quantised vortices and inherently possesses vortex reconnections is the Gross-Pitaevskii (GP) equation that describes the evolution of the superfluid order parameter ψ.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 in the GP equation is order of the healing length ξ, 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 re-arrange themselves following reconnections into complex tangles with non-trivial 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 modelling is needed (unlike the VF model).However GP does not provide direct information on vortices.
In this Letter we apply a novel numerical algorithm [10] to track accurately the configuration of a turbulent vortex tangle evolving accordingly to the GP model.We focus on the evolution and decay of the tangle.Firstly, 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 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 the L'vov-Nazarenko theoretical prediction [12].
The GP model for the condensate wave-function ψ is given by where m is the mass of the bosons and g = 4πa 2 /m, with a the s-wave scattering length.Madelung's transformation ψ(x, t) = ρ(x, t)/m exp [i m φ(x, t)] relates the wave-function ψ to a superfluid of density ρ(x, t) and velocity v = ∇φ.The quantum of circulation about the ψ = 0 vortex lines is Γ = h/m.When Eq.( 1) is linearised about a constant value ψ = ψ0 , the sound velocity is given by c = (g| ψ0 | 2 /m) 1/2 with dispersive effects taking place at length scales smaller than the healing length ξ = ( 2 /2m| ψ0 | 2 g) 1/2 .In the simulations presented in this Letter, the mean density is fixed to the unity and the physical constants in Eq.( 1) are determined by the values of ξ and c = 1.The quantum of circulation results in Γ = 4πc ξ/ √ 2. Numerical integration of Eq.( 1) is performed using a standard pseudo-spectral code.We integrate the so-called Taylor-Green flow [13] with no enforced symmetries at resolutions 256 3 and 512 3 with ξ = 2π/256 and ξ = 2π/512 respectively (see Appendix B).With these units the largest eddy turnover time is 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 easily spotted by plotting the low-value iso-surfaces of the density field as displayed in Fig. 1a.Low-density regions corresponding to vortex lines are plotted in red, while density fluctuations (sound) are rendered in light blue.We track the vortex lines forming the tangle with a recently-developed algorithm [10].Vortex lines are followed using the pseudovorticity field [14] and the exact vortex position is obtained by a Newton-Raphson method.The algorithm is robust and accurate as it takes full advantage of the spectral resolution of the field.It allows for identifying 2 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 modelling is needed (unlike the VF model).However GP does not provide direct information on vortices.
In this Letter we apply a novel numerical algorithm [10] to track accurately the configuration of a turbulent vortex tangle evolving accordingly to the GP model.We focus on the evolution and decay of the tangle.Firstly, we show that after the onset of turbulence, the vortex line density satisfies the Vinen's decay law [11] with a coe cient that is in agreement with He II.Di↵erent 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 the L'vov-Nazarenko theoretical prediction [12].
The GP model for the condensate wave-function is given by where m is the mass of the bosons and g = 4⇡a~2/m, with a the s-wave scattering length.Madelung's transformation (x, t) = p ⇢(x, t)/m exp [i m ~ (x, t)] relates the wave-function to a superfluid of density ⇢(x, t) and velocity v = r .The quantum of circulation about the = 0 vortex lines is = h/m.When Eq.( 1) is linearised about a constant value = ˆ 0 , the sound velocity is given by c = (g| ˆ 0 | 2 /m) 1/2 with dispersive e↵ects taking place at length scales smaller than the healing length In the simulations presented in this Letter, the mean density is fixed to the unity and the physical constants in Eq.( 1) are determined by the values of ⇠ and c = 1.The quantum of circulation results in = 4⇡c ⇠/ p 2. Numerical integration of Eq.( 1) is performed using a standard pseudo-spectral code.We integrate the so-called Taylor-Green flow [13] with no enforced symmetries at resolutions 256 3 and 512 3 with ⇠ = 2⇡/256 and ⇠ = 2⇡/512 respectively (see Supplemental Material for details).With these units the largest eddy turnover time is 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 easily spotted by plotting the low-value iso-surfaces of the density field as displayed in Fig. 1a.Low-density regions corresponding to vortex lines are plotted in red, while density fluctuations (sound) are rendered in light blue.We track the vortex lines forming the tangle with a recently-developed algorithm [10].Vortex lines are followed using the pseudo- vorticity field [14] and the exact vortex position is obtained by a Newton-Raphson method.The algorithm is robust and accurate as it takes full advantage of the spectral resolution of the field.It allows for identifying separately each single line forming the tangle.Figure 1b shows the corresponding tracked tangle displaying in di↵erent colors all the 496 vortex rings (see Supplemental Material for a movie showing the full evolution and details of the algorithm).
During the decay, vortices radiate phonons that populate the 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 [15].In superfluids such a decay is modelled by the Vinen equation [11] for vortex line density L: where 2 is a constant of order of unity.Its solutions manifest a L ⇠ t 1 behaviour at long times: this powerlaw decay has been named quantum turbulent decay and measured in He II experiments [16] and BS numerical simulations [17].In Fig. 2a we show the temporal evolution of 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 L obtained by computing the ratio between the volume of points such that ⇢(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 [18].Even if this technique is able to capture the qualitative behaviour of L, it fails to grasp at long times the power-law predicted by Vinen's equation.This is show in the inset of Fig. 2a where the measured L(t) = (L(t) 1 L(t 0 ) 1 ) 1 , setting t 0 = 17, is compared to Vinen's prediction.We can intuitively explain this discrepancy by reasoning that the vortex core size separately each single line forming the tangle.Figure 1b shows the corresponding tracked tangle displaying in different colors all the 496 vortex rings.
During the decay, vortices radiate phonons that populate the 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 [15].In superfluids such a decay is modelled by the Vinen equation [11] for vortex line density L: where χ 2 is a constant of order of unity.Its solutions manifest a L ∼ t −1 behaviour at long times: this powerlaw decay has been named quantum turbulent decay and measured in He II experiments [16] and BS numerical simulations [17].In Fig. 2a we show the temporal evolution of 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 L obtained by computing the ratio between the volume of points such that ρ(x) < 0.2 and the corresponding surface of a perfect two-dimensional vortex profile.This latter method is a standard technique within GP to estimate L numerically [18].Even if this technique is able to capture the qualitative behaviour of L, it fails to grasp at long times the power-law predicted by Vinen's equation.This is show in the inset of Fig. 2a where the measured ∆L(t) = (L(t) −1 − L(t 0 ) −1 ) −1 , setting t 0 = 17, is compared to Vinen's prediction.This discrepancy can be explained by the presence of sound waves and small scale KWs not detected by the estimation.The tracking data also allow for determining the numerical constant χ 2 = 0.65.This value is in remarkable agreement with experimental values measured in He II in the low temperature limit [11].(proportional to the uniform condensate state) varies in time due to the fact that more and more sound excitations are created by the superfluid decay, altering the estimation of L by fixing the (non time-dependant) density threshold.The tracking data also allow for determining the numerical constant 2 = 0.65.This value is in remarkable agreement with experimental values measured in He II in the low temperature limit [11].
From the movie and snapshots provided in the Supplemental Material, it is clear that the complexity of tangle first increases and then decreases.This observation is intuitively explained by the idea that vortex lines simultaneously stretch, bend and coil during reconnection events.The complexity of tangle can be measured by computing the changes in some of its algebraic and topological quantities, such as the total number of rings, the average crossing number, the total linking Lk and writhe W r [8].Explicit definitions are given in the Supplemental Material.The total number of rings and the average crossing number vary, as expected, in a similar manner to the total vortex length in Fig. 2a (See Fig. 3 in the Supplemental Material).Here we focus on the central line helicity H c / 2 = Lk + W r [19].This quantity is the analogous to helicity in classical fluid dynamics, an important inviscid invariant.The linking number Lk 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 knots and KWs [20].Figure 2b shows the temporal evolution of these three quantities.Initially, Lk = W r = 0, as expected for the Taylor-Green flow.Surprisingly, during the evolution Lk becomes non-zero, indicating the presence of linked rings, such as the ones displayed in Fig.2c [21].This is remarkable as in the GP model su ciently simple vortex configurations usually decay by reducing their linking number [22].Once the decay is established, no linked rings are present and only writhe contributes to H c .The center-line helicity fluctuates about a zero mean, an indication of the presence of KWs.KWs have already been indirectly observed in the Taylor-Green flow during the turbulent stage [23], 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.We explore the time behaviour of the probability density functions (PDFs) of the curvature  and torsion ⌧ of the entire set of vortices in the system.In Fig. 3a we present the PDF of curvature, normalised by its mean value, at di↵erent stages.The temporal evolution of the mean curvature hi and its rms value  rms are also displayed in the inset.We can observe that hi 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 excep- From Fig. 6 in Appendix A, it is clear that the complexity of tangle first increases and then decreases.This observation is intuitively explained by the idea that vortex lines simultaneously stretch, bend and coil during reconnection events.The complexity of tangle can be measured by computing the changes in some of its algebraic and topological quantities, such as the total number of rings, the average crossing number, the total linking Lk and writhe W r [8].Explicit definitions are given in Appendix C. The total number of rings and the average crossing number vary, as expected, in a similar manner to the total vortex length in Fig. 2a (See Fig. 3 in the  Appendix C).Here we focus on the central line helicity . This quantity is the analogous to helicity in classical fluid dynamics, an important inviscid invariant.The linking number Lk 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 knots and KWs [20].Figure 2b shows the temporal evolution of these three quantities.Initially, Lk = W r = 0, as expected for the Taylor-Green flow.Surprisingly, during the evolution Lk becomes non-zero, indicating the presence of linked rings, such as the ones displayed in Fig. 2c [21].This is remarkable as in the GP model sufficiently simple vortex configurations usually decay by reducing their link- (proportional to the uniform condensate state) varies in time due to the fact that more and more sound excitations are created by the superfluid decay, altering the estimation of L by fixing the (non time-dependant) density threshold.The tracking data also allow for determining the numerical constant 2 = 0.65.This value is in remarkable agreement with experimental values measured in He II in the low temperature limit [11].
From the movie and snapshots provided in the Supplemental Material, it is clear that the complexity of tangle first increases and then decreases.This observation is intuitively explained by the idea that vortex lines simultaneously stretch, bend and coil during reconnection events.The complexity of tangle can be measured by computing the changes in some of its algebraic and topological quantities, such as the total number of rings, the average crossing number, the total linking Lk and writhe W r [8].Explicit definitions are given in the Supplemental Material.The total number of rings and the average crossing number vary, as expected, in a similar manner to the total vortex length in Fig. 2a (See Fig. 3 in the Supplemental Material).Here we focus on the central line helicity H c / 2 = Lk + W r [19].This quantity is the analogous to helicity in classical fluid dynamics, an important inviscid invariant.The linking number Lk 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 knots and KWs [20].Figure 2b shows the temporal evolution of these three quantities.Initially, Lk = W r = 0, as expected for the Taylor-Green flow.Surprisingly, during the evolution Lk becomes non-zero, indicating the presence of linked rings, such as the ones displayed in Fig.2c [21].This is remarkable as in the GP model su ciently simple vortex configurations usually decay by reducing their linking number [22].Once the decay is established, no linked rings are present and only writhe contributes to H c .The center-line helicity fluctuates about a zero mean, an indication of the presence of KWs.KWs have already been indirectly observed in the Taylor-Green flow during the turbulent stage [23], 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.We explore the time behaviour of the probability density functions (PDFs) of the curvature  and torsion ⌧ of the entire set of vortices in the system.In Fig. 3a we present the PDF of curvature, normalised by its mean value, at di↵erent stages.The temporal evolution of the mean curvature hi and its rms value  rms are also displayed in the inset.We can observe that hi 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 excep- ing number [22].Once the decay is established, no linked rings are present and only writhe contributes to H c .The center-line helicity fluctuates about a zero mean, an indication of the presence of KWs in this turbulent tangle, as indirectly observed in [23].
We now study statistical properties of some geometrical quantities of the vortex filaments.We explore the time behaviour of the probability density functions (PDFs) of the curvature κ and torsion τ .In Fig. 3a we present the PDF of curvature, normalised by its mean value, at different stages.The temporal evolution of the mean curvature κ and its rms value κ rms are also displayed in the inset.We can observe that κ increases rapidly at early stages and then almost saturates.The rms value of the curvature presents the same tendency with the exception of peaks.These are evidences of reconnection events where high values of curvature are found in localised regions.It is worth noticing that the PDFs, rescaled by their mean curvature, exhibit a relatively good collapse to a self-similar form.This latter observation indicates a power-law behaviour ∼ κ 1 at small curvature values, while an exponentially-decaying tail is present at large curvature values.A similar behaviour has also been observed within the VF model [24].In Fig. 3c 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 an universal power-law behaviour of τ −3 at all times, meaning that the second and higher moments of the torsion diverge during the decay.fect straight line in the GP model [29].We highlight that although the weak-wave turbulence prediction for the KW spectrum was 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 loose their energy by phonon radiation and also 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 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 helicity.It will be of great interest to repeat a similar analysis in a GP setting where the mean helicity of the flow is not zero, like the ABC flow introduced in [20] where linking and selflinking processes could be substantially enhanced.Overall, the results presented in this Letter 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 [30,31].A controlled exper- ing appears in vortex tangles of random wave fields that are solutions of the Helmholtz equation [25].This might be an indication that for one-time small-scale quantities, quantum turbulent tangles could be interpreted as random vortices.
The large curvature fluctuations and the torsion fluctuation about a zero mean suggest also the presence of KWs at all scales propagating on quasi-planar 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 ∼ k −α (here k is the Kelvin wavenumber) and to explain the energy transfer through KW scales.Assuming small amplitude KWs (weak non-linearity), Kozik&Svistunov [26] and L'vov&Nazarenko [12] obtained the exponents α KS = 17/5 and α LN = 11/3 respectively considering two different orders of interaction.We compute the KW spectra of the 50 largest rings during the evolution of the tangle applying a Gaussian kernel in order to establish the unperturbed ring (see Appendix D).The spectra, averaged over the rings, are shown for different times in Fig. 4a.It is evident that all accessible KW modes get populated at early times due to reconnection events that trigger the cascade [27].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 ≤ t ≤ 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 wavenumber decades.In Fig. 4b we show the spectrum at t ∼ 5: the observed power-law exponent agrees with the L'vov&Nazarenko α LN = 11/3 prediction.This can be better appreciated by looking at the compensated spectra with respect to α LN and α KS showed in the inset.This finding corroborates the result in favour of L'vov&Nazarenko's prediction previously obtained while studying the KW oscillations about a perfect straight line in the GP model [28].We highlight that although the weak-wave turbulence prediction for the KW spectrum was 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 loose their energy by phonon radiation and also 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 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 helicity.It will be of great interest to repeat a similar analysis in a GP setting where the mean helicity of the flow is not zero, like the ABC flow introduced in [20] where linking and selflinking processes could be substantially enhanced.Overall, the results presented in this Letter 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 [29,30].A controlled experimental setting with a turbulent BEC, such as the one presented in this Letter, has yet to be achieved but it should be realisable in the near future.
GK, DP and AV were supported by the cost-share Royal Society International Exchanges Scheme (ref. IE150527) in conjunction with CNRS.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.

Appendix A. Vortex Tracking
We have recently developed a robust and accurate algorithm to track vortex lines in the Gross-Pitaevskii equa-tion (GP) with arbitrary geometries in a periodic domain.The full details of the algorithm and the case studies to check its validity can be found in [10].We recall here briefly the basic ideas.A quantised vortex line is defined by the nodal lines of the wavefunction.In three dimensions this corresponds to a line defined by Re[ψ(x, y, z)] = Im[ψ(x, y, z)] = 0 (3) The algorithm is based on a Newton-Raphson method to find zeros of ψ and on the knowledge of the pseudovorticity field W = ∇Re[ψ] × ∇Im[ψ], always tangent to the line, to follow vortex lines (in the spirit of Rorai et al. [14]).
Starting from a point x 0 where the density |ψ| 2 is below a given small threshold (therefore very close to a vortex), we define the orthogonal plane to the vortex line using W(x 0 ).The plane is then spanned by the two directors û1 and û2 as illustrated in Fig. 5 ter approximation for the true point x v where the vortex lies on the plane is then given by x 1 = x 0 + δx.Here the increment δx is obtained using the Newton-Raphson formula: 0 = ψ(x 0 + δx) ≈ ψ(x 0 ) + J(x 0 )δx, (4) where J(x 0 ) is the Jacobian matrix expressed as The increment can be therefore found using δx = −J −1 (x 0 ) • (Re[ψ(x 0 )], Im[ψ(x 0 )]) T .The Jacobian matrix is a non-singular 2 × 2 matrix so its inverse is computed with its explicit formula.We underline that the method will require in general the evaluation of the Jacobian (5) at intermesh points.However, making use of the spectral representation of ψ, we can compute spatial derivatives of the field ψ efficiently using fast Fourier transforms.This process can be iterated until the true vortex location x v is determined upon a selected convergence precision.
To track the next vortex point we use as next initial guess x 0 = x v + ζW, which is obtained evolving along ω ps by a small step ζ.The process is repeated until the entire line is tracked and closed.Once the line is completely tracked the process is repeated with another line until the whole domain has been fully explored.The algorithm has been extensively described and validated in [10] using different case tests.As a final remark, we underline that the algorithm takes full advantage of the spectral resolution when a field is generated integrating the Gross-Pitaevskii with a pseudo-spectral method.A comparison between the low-value isosurfaces of the density field and tracked vortices is displayed in Fig. 6 for different times of the Taylor-Green decay.

FIG. 1 .
FIG. 1. (Color online) a) Isosurface of density field at t = 12.Low-density regions corresponding to vortex lines are plotted in red, while density fluctuations (sound) are rendered in light blue.b) Corresponding vortex tracking.Di↵erent colors correspond to di↵erent vortices.Resolution 256 3 .

FIG. 1 .
FIG. 1. (Color online) a) Isosurface of density field at t = 12 rendered with VAPOR software.Low-density regions corresponding to vortex lines are plotted in red, while density fluctuations (sound) are rendered in light blue.b) Corresponding vortex tracking.Different colors correspond to different vortices.Resolution 256 3 .
FIG. 2. (Color online) a) Temporal evolution of the vortex line density: tracked data are plotted using blue circles, volume estimation in solid red line.The inset shows the long time decay of L (see text), together with Vinen's prediction L Vinen = [ 2 2⇡ (t t 0 )] 1 , with = 0.65 (solid black line).b) Temporal evolution the total linking Lk, writhe W r and center-line helicity H c .c) Visualisation of two linked rings at t = 21.Resolution 256 3 .

3 FIG. 3 .
FIG. 3. (Color online) a) PDFs of curvature  normalised by their respective mean values hi at di↵erent times (same legend as (b)).The inset displays the temporal evolution of the mean and rms values of . b) PDFs of torsion ⌧ at di↵erent times.The inset emphasises their ⌧ 3 power-law tail.Resolution 256 3 .

FIG. 2 .
FIG. 2. (Color online) a) Temporal evolution of the vortex line density: tracked data are plotted using blue circles, volume estimation in solid red line.The inset shows the long time decay of ∆L (see text), together with Vinen's prediction ∆LVinen = [χ2 Γ 2π (t − t0)] −1 , with χ = 0.65 (solid black line).b) Temporal evolution the total linking Lk, writhe W r and center-line helicity Hc. c) Visualisation of two linked rings at t = 21.Resolution 256 3 .
FIG. 2. (Color online) a) Temporal evolution of the vortex line density: tracked data are plotted using blue circles, volume estimation in solid red line.The inset shows the long time decay of L (see text), together with Vinen's prediction LVinen = [ 2 2⇡ (t t0)] 1 , with = 0.65 (solid black line).b) Temporal evolution the total linking Lk, writhe W r and center-line helicity Hc. c) Visualisation of two linked rings at t = 21.Resolution 256 3 .

3 FIG. 3 .
FIG. 3. (Color online) a) PDFs of curvature  normalised by their respective mean values hi at di↵erent times (same legend as (b)).The inset displays the temporal evolution of the mean and rms values of . b) PDFs of torsion ⌧ at di↵erent times.The inset emphasises their ⌧ 3 power-law tail.Resolution 256 3 .

FIG. 3 .
FIG. 3. (Color online) a) PDFs of curvature κ normalised 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 emphasises their τ −3 power-law tail.Resolution 256 3 .

FIG. 5 .
FIG. 5. Sketch of the plane on which the Newton-Raphson method is implemented.

FIG. 6 .
FIG. 6. (Color online) (left) Isosurfaces of density field at different times rendered by VAPOR (https://www.vapor.ucar.edu).Low-density regions corresponding to vortex lines are plotted in red, while density fluctuations (sound) are rendered in light blue.(right) Corresponding vortex tracking.Different colours correspond to different vortices.Resolution 256 3 .
The same scal-