Vortex clustering, polarisation and circulation intermittency in classical and quantum turbulence

The understanding of turbulent flows is one of the biggest current challenges in physics, as no first-principles theory exists to explain their observed spatio-temporal intermittency. Turbulent flows may be regarded as an intricate collection of mutually-interacting vortices. This picture becomes accurate in quantum turbulence, which is built on tangles of discrete vortex filaments. Here, we study the statistics of velocity circulation in quantum and classical turbulence. We show that, in quantum flows, Kolmogorov turbulence emerges from the correlation of vortex orientations, while deviations—associated with intermittency—originate from their non-trivial spatial arrangement. We then link the spatial distribution of vortices in quantum turbulence to the coarse-grained energy dissipation in classical turbulence, enabling the application of existent models of classical turbulence intermittency to the quantum case. Our results provide a connection between the intermittency of quantum and classical turbulence and initiate a promising path to a better understanding of the latter.

Vortices are manifestly the most attractive feature of fluid flows occurring in Nature.They are highly rotating zones of the fluid that often take the form of elongated filaments, of which tornadoes are one prominent example in atmospheric flows.Such structures can travel and interact with other vortex filaments, as well as with the surrounding fluid.In fact, the dynamics of vortex filaments in fluid flows is highly non-trivial, as they can reconnect changing the topology of the flow 1 .Their non-trivial arrangements may lead to very complex configurations and in particular to turbulence, an out-of-equilibrium state characterised by a large scale separation between the scale at which energy is injected and the one at which it is dissipated.In three-dimensional flows, because of the inherently non-linear character of turbulence, energy initially injected at large scales is transferred towards the small scales through a cascade-like process.
In turbulent flows, the typical thickness of a vortex filament is comparable to the smallest active scale of turbulence 2 , itself usually much smaller than the eddies carrying most of the energy content of the flow.Vortex filaments may thus be seen as the fundamental structure of turbulence, whose collective dynamics leads to the multi-scale complexity of such flows.Indeed, depending on their individual intensities and orientations, a set of vortex filaments located within a given spatial region may contribute constructively or destructively to the fluid rotation rate.In fluid dynamics, the rotation rate of a two-dimensional fluid patch is commonly quantified by the velocity circulation around the closed loop C surrounding the patch, where v is the fluid velocity field.Note that, by virtue of Stokes' theorem, the circulation is equal to the flux of vorticity,  = ∇ × v, through the fluid patch.
The above view of vortex filaments as the fundamental unit of fluid flows is particularly appropriate in superfluids, such as low-temperature liquid helium and Bose-Einstein condensates (BECs).Indeed, in such fluids, vortices are well-defined discrete objects about which the circulation is quantised, taking values multiple of κ = ℎ/.Here ℎ is Planck's constant and  is the mass of the bosons constituting the superfluid 3 .Such property arises from their quantum nature, as vortices are topological defects of the macroscopic wave function describing the system.For this reason, vortex filaments in superfluids are called quantum vortices.
One of the most striking properties of low-temperature superfluids is their total absence of viscosity.Despite this fact, quantum vortex reconnections are possible, since Helmholtz' theorem that forbids reconnections in classical inviscid fluids 1 breaks down due to the vanishing fluid density at the vortex core.This picture was first suggested by Feynman in the 50's 4 and later confirmed numerically in the framework of the Gross-Pitaevskii (GP) equation 5 .Since then, quantum vortex reconnections have been observed experimentally in superfluid helium 6 and in BECs 7 .They are characterised by universal scalings laws 8,9 and have been linked to irreversibility, both in experiments 10 and in numerical simulations 11 .In the early vortex filament simulations by Schwarz 12 , it was noticed that quantum vortex reconnections are a key physical process for the development of quantum turbulence, a state described by the complex interaction of a tangle of quantum vortices.Such a state is illustrated by the vortex filaments (in green and yellow) visualised in Fig. 1, obtained from the GP simulations performed in Ref. 13.
Quantum turbulence is characterised by a rich multi-scale physics.At small scales, between the vortex core size (about 1 Å in superfluid 4 He) and the mean inter-vortex distance ℓ (∼ 1 µm), the physics is governed by the dynamics of individual quantised vortices 14 .At such scales, Kelvin waves (waves propagating along vortices) and vortex reconnections are the main physical processes carrying energy along scales 15,16 .In contrast, at scales larger than ℓ, the quantum nature of the superfluid becomes less important and a regime comparable to classical turbulence emerges.Indeed, at such scales, a Kolmogorov turbulent cascade is observed, provided that a large scale separation exists between ℓ and the largest scale of the system.In particular, the scaling law predicted by Kolmogorov's celebrated K41 theory 17 for the kinetic energy spectrum has been observed in superfluid helium experiments 18,19 and in numerical simulations of quantum turbulence [20][21][22] .
Previous studies have suggested that, in quantum turbulence, the emergence of K41 scaling laws is associated to a local polarisation of the vortex tangle 14,[23][24][25][26][27] .In other words, within a given spatial region, the orientations of nearby vortices are not independent, but instead have some degree of correlation.This phenomenon is visible in Fig. 1, where vortex bundles -regions of same-coloured vortex filaments -can be clearly identified.This local polarisation is present even in ideally isotropic flows, and should not be confused with the preferential large-scale orientation of vortices, which typically occurs in anisotropic flows.A classical example of the latter is a rotating cylindrical vessel filled with superfluid helium 4 .
In a recent work 13 , we have shown that the quantitative similarities between classical and quantum turbulence go far beyond the Kolmogorov energy spectrum.Indeed, both systems display the emergence of extreme events that result in the break-down of Kolmogorov's K41 theory -a phenomenon known as intermittency.Our work was motivated by the recent study of Iyer et al. [28] , which suggested that intermittency has a relatively simple signature on the statistics of circulation in classical turbulence.In particular, the moments of the circulation measured over fluid patches of area A ∼  2 follow a power law of the form with scaling exponents λ  that increasingly deviate from the K41 prediction λ K41  = 4/3 as the moment order  increases.By performing simulations of a generalised GP equation, we have shown that the anomalous scaling exponents λ  in the inertial scales of quantum turbulence closely match those observed in classical turbulence 13 .Note that, up until now, most of the advances in the understanding of intermittency have been made in terms of velocity increments.However, despite many theoretical efforts 17,[29][30][31] , there is still no first-principles theory able to explain this phenomenon.The above-cited findings suggest that circulation may provide an alternative path towards a better understanding of turbulence (as first hinted the by pioneering theoretical work of Migdal [32] ), and eventually, to novel circulation-based theories of intermittency 33,34 .
The strong similarity between the statistics of circulation in classical and quantum turbulence is particularly striking given the very different nature of vortices in both types of fluids.This statistical equivalence opens the way for an interpretation of the intermittency of classical turbulent flows in terms of the collective dynamics of discrete vortex filaments carrying a fixed circulation.With this idea in mind, we relate in this work the intermittent statistics of velocity circulation in classical and quantum turbulence.We start by investigating in quantum turbulence how local vortex polarisation, as well as the non-trivial spatial distribution of vortex filaments, affect circulation statistics.We address the following questions: Is it possible to study both effects separately?Do they contribute in the same way to the flow intermittency?We then provide a relation between the spatial distribution of discrete vortices, and the coarse-grained energy dissipation rate in classical turbulence, at the core of existent intermittency models.
In this work, quantum and classical turbulent systems are respectively studied using high-resolution direct numerical simulations of a generalised GP and the incompressible Navier-Stokes (NS) equations.Discrete vortices and their signs are extracted from the GP fields and then analysed.To disentangle the effects of polarisation and spatial vortex distribution, we additionally study a disordered turbulence state.Such state is generated from the discrete vortex data by randomly resetting the sign of each individual vortex while keeping its position fixed.To illustrate the differences between the turbulent (non-disordered) and the disordered turbulence states, we plot in Fig. 2 the kinetic energy spectrum associated to each vortex configuration (see Methods for details on the computation of the spectra from discrete vortices).First, we see that the turbulent case displays a clear  −5/3 range, in agreement with the energy spectra obtained from the full GP and NS fields.Note that, in the case of GP fields, we show the incompressible kinetic energy spectrum, which contains 86% of the total energy of the system -the other components being the compressible, internal and quantum energy 20,22 .Secondly, the K41 scaling disappears once polarisation is artificially suppressed from the tangle, leading to a trivial  −1 scaling range for the disordered state (see Methods for a brief derivation).Note that this same scaling has already been observed in vortex filament simulations, once the vortex tangle has been decomposed into polarised and random components 26 .

Results
A simple discrete model of circulation.Let us first consider a set of  discrete vortices, each of them carrying a circulation κ  , where   = ±1 is the sign of each vortex.From now on we set the quantum of circulation to κ = 1 for simplicity.We propose to model the total circulation of the -vortex collection, Γ  =  =1   , as a biased one-dimensional random walk.Polarisation is naturally introduced by letting each random step   be positively correlated with the instantaneous position Γ −1 , i.e. the total circulation of all previous vortices.
Concretely, we construct inductively the following toy model for the circulation.The sign of the first vortex,  1 , has equal probability of being positive or negative.Then, the sign of vortex  + 1 is positive with a probability  +1 , which we set to depend on the total circulation at step  as  the vertical axis is in arbitrary units.In the GP case, the incompressible part of the kinetic energy is plotted 35 .Also shown are the energy spectra obtained after applying the vortex detection procedure to the GP fields (see Methods for details), both before and after the randomisation of the vortex orientations (turbulent and disordered cases, respectively).Dashed and dash-dotted lines respectively represent the Kolmogorov scaling  −5/3 and the disordered scaling ∈ [0, 1] at each step .For the sake of simplicity, we choose here  () = β (see the Supplementary Information for the general case), where β ∈ [0, 1] is an adjustable parameter that sets the polarisation of the system.When β = 0, one retrieves a standard random walk with scaling |Γ| 2  ∼ .Conversely, for β = 1, one recovers a fully polarised set of vortices behaving as |Γ| 2  ∼  2 .The resulting model is a discrete Markov process, since the probability distribution of Γ +1 only depends on the state {, Γ  } via the probability  +1 .Concretely, the probability P  (Γ) of having Γ  = Γ obeys the master equation Multiplying this equation by Γ 2 , summing over all Γ and, for the sake of simplicity, taking the limit of continuous , one gets a closed equation for the circulation variance, where averages are over all realisations after  steps.For large , this equation predicts the scaling Γ 2  ∼  for β < 1/2 (corresponding to a set of vortices with negligible polarisation), and Γ 2  ∼  2β otherwise.In particular, choosing β = 2/3, one recovers the Kolmogorov scaling by replacing  ∝ A. This relation between the number of vortices  and the loop area A containing them is expected to hold on average under spatial homogeneity conditions, but neglects potentially important inhomogeneities in the spatial vortex distribution that may affect high-order moments of  (see also Fig. 5).Besides, for the -th order moment, the model predicts the self-similar scaling |Γ|   ∼  γ  with γ  = β.More generally, for any suitable function  () defining the probability   of the model, one obtains the linear scaling γ  =  min {max {1/2,  (0)}, 1} (see the Supplementary Information for more details on the calculations).
The toy model introduced above shows in a very simple manner how a specific correlation (or polarisation) is responsible for the emergence of non-trivial scaling laws, as already suggested by previous works on quantum turbulence [24][25][26] .In addition, the model yields self-similar statistics, suggesting that polarisation is not sufficient to reproduce the observed intermittency of circulation in classical 28 and quantum 13 turbulent flows.At this point, we may speculate that the lack of intermittency in the model is likely associated with the missing notion of space.Indeed, on average one expects to have a number of vortices  ∼ A/ℓ 2 crossing a loop of area A, where ℓ is the mean inter-vortex distance.Yet, fluctuations in their spatial distribution -associated to the appearance of vortex clusters and voids -may strongly influence high-order moments.As seen in Fig. 1, such effects clearly take place in turbulent flows, where they are linked to the formation of coherent structures.
Comparison with quantum turbulence data.The ideas hinted at by our toy model can be verified using actual quantum turbulence data.With this aim, we identify all vortex filaments present in our GP simulations (see Methods for details), and compute circulation statistics as a function of the number  of considered vortices.Crucially, groups of vortices are chosen based on their spatial proximity, which is required to preserve the correlation between vortices.On the other hand, with such a conditioning, one may expect the effect of strong spatial fluctuations of the vortex distribution to be somewhat relaxed.In practice, for each two-dimensional cut of the simulation, we consider sets of  neighbouring vortices in order to compute the circulation moments |Γ|   .Then, to improve the statistics, we repeat such measurement for each cut and along the three Cartesian directions.
The resulting second-order moment Γ 2  is shown in Fig. 3(a), along with the moment Γ 2 A measured for different loop areas A (data from Müller et al. [13] ).At small scales, Γ 2 A ∼ A 1 due to the discrete nature of vortices 13 .In contrast, within the inertial range, both moments clearly exhibit the expected Kolmogorov scaling.In particular, Γ 2  ∼  γ 2 with γ 2 = 4/3.This result allows us to use the extended self-similarity (ESS) framework 36 to determine the scaling properties of higher-order moments via the relation

𝑛
. Remarkably, as shown in Figs.3(b-c), the moments display a clear self-similar behaviour with γ  = 2/3, thus obeying Kolmogorov scaling for all orders.The self-similarity is also observed in the normalised probability density functions (PDFs) of Γ for different values of  [Fig.3(d)], which nearly collapse and are close to Gaussian.This behaviour should be contrasted with the non-collapsing PDFs of Γ for different loop areas A [Fig. 3(e)].Note that, in both cases, the chosen values of A and  lay within the inertial range, represented by a grey background in Figs.3(a-c).

Disentangling polarisation and spatial vortex distribution.
The fitted scaling exponents for the turbulent case 2γ turb  , discussed above, are plotted in Fig. 4 (blue filled stars) as a A/ 2 ; n 10 -3 10 -2 10 -1 10 0 10 1 10 2 10 3 function of the moment order .These exponents are compared to the measured values of λ turb  (blue filled circles) obtained according to Eq. ( 2), where averages are performed for different loop areas A. The latter are the same as in Ref. 13.The factor 2 in front of γ turb  comes from considering the relation  ∼ A ∼  2 .As discussed earlier, the moments averaged for different  closely follow the self-similar K41 scaling 2γ turb  ≈ 4/3 (blue solid line), while the λ turb  exponents -affected by the spatial vortex distribution -show signs of intermittency 13 .
To further distinguish the effects of polarisation and spatial vortex distribution on circulation statistics, we perform the following numerical experiment.We recompute the circulation in the quantum turbulent flow, but before doing this, we randomise the sign of each vortex on each analysed twodimensional cut while keeping its position fixed.By doing this, we get rid of the system polarisation, while maintaining the non-trivial spatial distribution of vortices.We refer to this system as disordered turbulence.In our non-intermittent toy model, this setting would correspond to the unpolarised value β = 0, yielding the self-similar circulation scaling |Γ|   ∼  /2 .In Fig. 4, we display the corresponding measured exponents of the disordered state λ diso  and 2γ diso  (red unfilled markers).Remarkably, even after suppressing vor- tex polarisation, λ diso  also presents intermittency deviations.In contrast, the scaling exponents γ diso  satisfy the expected self-similar behaviour 2γ diso  ≈  (red solid line).The previous results suggest that the non-trivial polarisation of vortices, while being responsible for Kolmogorov scalings, has no major influence on the intermittency of the system.Furthermore, they indicate that the latter originates from fluctuations of the spatial distributions of vortices.From our above observations, one may therefore expect the scaling exponents of the circulation to be given by a composition of the polarisation and spatial distribution effects.That is, we may conjecture that the scaling exponents λ  and γ  are related by where  is some yet unknown function.
In order to check this idea, we can try to relate the scaling exponents of the turbulent and disordered turbulent systems.If relationship (5) were to hold true, one should have that λ diso  = λ turb 3 /4 .Using this relation with the measured exponents of the turbulent case, one indeed recovers the intermittency exponents λ diso  of the disordered case, as shown by the green squared markers in Fig. 4.This result strongly highlights the importance of the fluctuations of vortex concentration on the intermittency of circulation.
Spatial vortex distribution and OK62 theory.As a first step towards relating the intermittency of classical and quantum turbulence, we now quantify the spatial distribution of vortices in the latter system.If vortices were homogeneously distributed in space, then the number  of vortices within loops of area A would be expected to follow a Poisson distribution with mean value  A ∝ A. In that case, the moments of  would scale as   A ∼ A  for sufficiently large A. Equivalently, the number of vortices per unit area Z A = /A would follow the trivial scalings Z  A A ∼ 1 for all  > 0.
As shown in Fig. 5(a), this is clearly not the case, indicating that the spatial distribution of vortices is non-trivial in quantum turbulence (as may be inferred from the visualisation of Fig. 1).Indeed, while the first-order moment recovers a constant (consistently with the relation  A ∼ A), higher-order moments of Z A follow a different scaling with a negative exponent -a sign of anomalous behaviour.This is confirmed by the PDFs of  displayed in Fig. 5(b), which are long-tailed and strongly differ from a Poisson distribution (dashed line).
In classical turbulence, it is today well accepted that the intermittency of velocity fluctuations is linked to the emergence of violent events, characterised by strong spatial fluctuations of the kinetic energy dissipation rate ε(x).Such idea led Obukhov and Kolmogorov in 1962 to develop a refined similarity theory of turbulence, commonly referred to as OK62 theory, where such fluctuations are taken into account 17,37,38 , unlike K41 theory which only deals with the mean value of ε(x).This refined theory considers the scale-averaged (or coarse-grained) energy dissipation rate , where B(x, ) is a ball of radius  centred at x.When applied to the spatial velocity increments δ  over a distance , OK62 theory states that the statistics of δ  /(ε  ) 1/3 is self-similar and universal.Most intermittency models use ε  to predict the anomalous scaling of velocity increment statistics 17 .Some early experiments in classical turbulence showed that, when velocity increments are conditioned on the coarse-grained dissipation, their statistics becomes Gaussian 39,40 , proving that the intermittency of velocity fluctuations is hidden behind the distribution of energy dissipation.This observation was later confirmed by numerical simulations 41 .
In the case of low-temperature quantum turbulence, such as the one studied here, energy is taken away from the inertial range and transferred towards small scales by the Kelvin wave cascade and vortex reconnections 9,42 .Furthermore, the velocity field diverges at the vortex core, and thus the definition of the dissipation field is delicate.Nevertheless, we can give a phenomenological interpretation of the dissipation by assuming that the system is well represented as a dilute point-vortex gas.Such a picture was recently used by Apolinário et al. [33] to model the velocity circulation in classical turbulence, and becomes particularly pertinent in quantum fluids.Although the superfluid is inviscid, one can model small scale physics by some effective viscosity ν eff 23,43 , whose value is not important here.This approach allows us to directly estimate the coarse-grained dissipation field by using its classical definition in terms of velocity gradients and a Dirac-like supported vorticity field (see Supplementary Information).Given a disk of radius  crossed by  vortices, a straightforward calculation gives the estimate where ε  is the average of the local dissipation rate ε(x) over the disk, A = π 2 is the disk area, ξ the typical vortex thickness, and κ the quantum of circulation.The number of vortices per unit area Z A = /A would then be the quantum analogous of the coarse-grained dissipation ε  .Remarkably, and similarly to ε  -which is known to exhibit log-normal statistics in classical turbulence 44   To make a stronger connection between classical and quantum turbulence, we recall that the classical coarse-grained energy dissipation rate is a highly fluctuating quantity that presents anomalous scaling laws traditionally denoted by ε   ∼  τ( ) .It follows from Eq. ( 6) that the number of vortices should satisfy

𝑛 𝑝
A ∼  α( ) ∼ A α( )/2 , with α( ) = 2 + τ( ).(7)   Note that, because of homogeneity, τ(1) = 0, which translates as  A ∼ A for the mean number of vortices.In the classical turbulence literature, there are several multi-fractal models for the anomalous exponents τ( ) that are able to reproduce experimental and numerical measurements 17 .Among those, the She-Lévêque model 45 has one adjustable parameter D ∞ corresponding to the fractal dimension of the most singular structures of the system.In the original model, which closely matches existent turbulence measurements 36,[45][46][47] , these structures are assumed to be vortex filaments, hence D ∞ = 1.The combination of prediction (7) with the original She-Lévêque model, represented by the green dashed lines in Fig. 5(a), is in good agreement with our quantum turbulence data for sufficiently large A, although some deviations due to the limited scaling range may be present.
Classical turbulence and conditioned circulation.We now apply some of the previous ideas to classical turbulence.We perform a direct numerical simulation of the Navier-Stokes equations in a statistically steady state at a Taylor-based Reynolds number of Re λ = 510.The simulation is performed using 2048 3 collocation points.We then compute the velocity circulation over planar square loops of area A =  2 , and, following the framework of the OK62 refined similarity hypothesis, we condition its statistics on the coarse-grained dissipation field ε  .The latter is obtained by averaging the local dissipation ε over the interior of each loop.See Methods for details on the numerical simulations and the data analysis.We first consider the unconditioned velocity circulation PDFs, shown in Fig. 6(a).The PDFs display heavy tails (associated with intermittency) which depend on the considered scale /λ, with λ the Taylor micro-scale.This is consistent with the classical turbulence simulations of Iyer et al. 28,48 .The PDF tails are strongly suppressed when the statistics is conditioned on low values of the local coarse-grained dissipation, ε  / ε ∈ [0.5, 1], as seen in Fig. 6(b).The suppression of intermittency is also manifest in Fig. 6(c), where the scaling exponents of circulation are displayed after conditioning on different intervals of ε  .With no conditioning (black crosses), the scaling exponents match those of Iyer et al. [28] , whereas when conditioning on low values of ε  the K41 self-similar scaling is recovered.
Note that the above conditioning is slightly different from the one presented in Fig. 3, as here we are conditioning both on the loop area A and on the value of ε  within such loops.In the case of quantum turbulence, the equivalent would be to study Γ  | A , i.e. to consider only loops of area A having  vortices.Such a double conditioning is very restrictive, as it requires a very large amount of statistics.Nevertheless, we perform a similar analysis, considering loops having a low, average and high number of vortices relative to the mean.The respective scaling exponents are displayed in Fig. 6(d).
We find that, for loops with low and average number of vortices, the self-similar K41 scaling is recovered, whereas for loops having large vortex concentrations the statistics is still intermittent.The lack of self-similarity in regions of high dissipation (in classical flows) or high vortex concentration (in quantum flows) hints at the idea that not all such events contribute equally to circulation statistics.
Can OK62 theory describe circulation intermittency?Considering the relation introduced in Eq. ( 6) and the fact that the number of vortices per unit area follows the same intermittent behaviour as ε  , one could try to apply OK62 theory to relate scaling exponents of circulation λ  with those of dissipation τ( ), as traditionally done for velocity increments.Within this reasoning, Γ ∼ ε 1/3   4/3 , yielding a OK62-based relation λ  = 4/3 + τ( /3).However, such a relation is in strong disagreement with our data (classical and quantum turbulence, see Supplementary Information) and with early Navier-Stokes studies 49 .Nevertheless, this disagreement is not in contradiction with the fact that the anomalous scaling of the number of vortices is well described by standard multifractal dissipation models (see Fig. 5).Indeed, if one considers a vortex dipole (two vortices of same magnitude and opposite sign), their contribution to large fluctuations of the local dissipation field and to velocity increments may be very important.On the other hand, for the circulation, the dipole contribution is exactly zero due to vortex cancellation.This fact suggests that not all extreme dissipation events result in extreme circulation values.In particular, intense circulation events would be correlated to those highly dissipative structures in turbulence which carry a strong vortex polarisation, such as vortex sheets or bundles (at scales  ℓ).Note that, in classical fluids, the idea of vortex filaments organising into groups forming vortex sheets is consistent with the recently-proposed sublayers vortex picture of dissipation 50 .
The previous observations motivate us to introduce a modified OK62 theory for the circulation ("mOK62" in the following), where the most relevant singular structures are not vortex filaments but structures of higher fractal dimension.
To check this idea, we adapt the She-Lévêque model τ SL ( ) [Eq.(8)] by setting D ∞ = 2.2 instead of 1.The chosen dimensionality exactly corresponds to the monofractal fit obtained by Iyer et al. [28] and Müller et al. [13] for the highorder circulation moments ( > 3) in classical and quantum turbulence, and, as suggested in the former work, it may be linked to the effect of wrinkled vortex sheets.Note that, for large , our mOK62 model simplifies to λ  ≈ 10 9  + (3−D ∞ ), which is equivalent to the monofractal fit by Iyer et al. [28] .In Fig. 4, it is shown that the adapted model matches strikingly well the anomalous exponents of circulation both in the turbulent and in the disordered cases for  > 3 (dashed lines), while for  < 3 there are some deviations.
The previous results provide a possible interpretation for the difference between the intermittency of velocity fluctuations and of circulation, based on the different topologies of the dissipative structures contributing to extreme events.We shall notice that an alternative interpretation is also possible, based on the recent works by Apolinário et al. [33] and Moriconi [34] .In this framework, the circulation should scale as , where ν  is Kraichnan's eddy viscosity 51 .The latter is found by assuming that the energy spectrum takes the form E() ∼ ε 2/3  −5/3+α (where α is an intermittency correction), yielding ν  ∼  4/3+α .Note that this phenomenological approach mixes a mean-field approximation for determining ν  with the fluctuations arising from ε 1/2  .Moreover, in its present form, it does not directly account for vortex cancellations.Nevertheless, when combined with the standard She-Lévêque model (with D ∞ = 1), this model provides an expression for the exponents λ  as accurate as our mOK62 model in the turbulent case.There is certainly a need to pursuit further investigations to understand how both models differ and complement each other.

Discussion
In this work, we have attempted at providing an interpretation for the intermittent statistics of velocity circulation in turbulent flows.We have done so by viewing turbulent flows as a polarised tangle of discrete and thin vortex filaments, each carrying a constant circulation.While this view is a priori only appropriate in low-temperature quantum fluids, we expect it to be a very pertinent model of classical turbulence, considering the strong similarities recently unveiled between both systems 13 .
By introducing and solving a simple toy model and by analysing data of GP quantum turbulence simulations, we have shown that, in discrete-vortex systems, the Kolmogorov self-similar scalings result from a partial polarisation of the vortices (in agreement with previous quantum turbulence studies), while the intermittency of circulation statistics is linked to the non-trivial (non-Poissonian) spatial distribution of vortices.In fact, within fluid patches of varying area A in the inertial range of scales, the number of vortices  is found to be the quantum equivalent of the coarse-grained dissipation ε  in classical turbulence, as they both follow the approximately log-normal distribution first hypothesised by the celebrated Obukhov-Kolmogorov OK62 theory for ε  44 .Quantitatively, we show that the intermittency of  is well described by the She-Lévêque model for ε  , confirming the strong equivalence between both observables.
It is important to remark that the quantum turbulence simulations presented in this work have been performed on periodic domains, and are based on the GP equation describing an ideal superfluid at very low temperature.In contrast, most superfluid turbulence experiments using liquid helium are performed in confined systems and at finite temperatures 18,19 , in a regime that may be described by a two-fluid model 52 .Early experimental studies showed that the signature of intermittency on velocity increments is nearly independent of the temperature, matching observations in classical fluids [53][54][55] .These observations were later contradicted by a recent experimental investigation, which showed an enhancement of velocity intermittency in the two-fluid regime compared to classical turbulence 56 , in agreement with previous numerical simulations of related models 57,58 .Compared to velocity increments, we expect the circulation to be a much more robust observable in quantum fluids, as it does not display singular behaviour in the vicinity of vortices 13 .For this reason, measuring the scaling properties of circulation in future experiments may help disambiguate existent contradictions, and provide a clearer answer on the intermittency of finite-temperature quantum turbulence.Recent experiments have made initial attempts at reconstructing Eulerian velocity fields from Lagrangian particle tracking measurements in turbulent superfluid helium 56 .Such a technique could be used in principle to measure the velocity circulation in superfluid helium, although addressing high-order statistics might still be challenging.However, note that such an approach is delicate because, due to the two-fluid nature of finite-temperature superfluid helium, particles may fail to capture important Eulerian flow features 59,60 , and further work is needed to determine its suitability.
Finally, using data from NS and GP simulations, we have confirmed that the classical OK62 theory does not fully account for the intermittency of the circulation in classical and quantum turbulence.We have provided an explanation based on the presumed topology of the turbulent structures that most contribute to extreme circulation events.We have then proposed a modified OK62 description of circulation, where relevant singular structures have a fractal dimension D ∞ ≈ 2.2 associated to vortex sheets 28 .This value differs from the dimensionality D ∞ = 1 of isolated vortex filaments, used in the modelling of velocity increment statistics 45 .Using this idea, we have shown that the intermittency of circulation is well reproduced by a modified version of the She-Lévêque model, bringing support to the vortex sheet interpretation first proposed by Iyer et al. [28] .All the previous ideas were additionally tested by introducing a disordered turbulence state, obtained by artificially suppressing vortex polarisation from a GP numerical simulation.
There are still some questions that remain open for future works.In particular, further investigation on the topology of relevant structures for the intermittency of circulation is required.We have argued that the smallest structures significant for circulation are vortex sheets, as simpler structures are irrelevant due to vortex cancellation.One way of approaching this topic is by use of cancellation exponents [61][62][63] , method that exploits the fact that circulation can take either negative or positive values.An alternative approach, proposed by an anonymous Referee, suggests that the most relevant singular structures for velocity circulation should still be vortex filaments.Further investigations on the fractal dimension of circulation would help develop more accurate models of intermittency.
Our findings hint at the existence of a coarse-grained quantity different from ε  , which may better encapsulate the intermittency of circulation in classical turbulence in the spirit of a OK62-like theory.Furthermore, it may be appropriate to investigate the relevance of quantities such as the local vorticity magnitude (or enstrophy) or the local strain.Such coarse-grained quantity would be expected to display intermittent statistics with extreme values associated to the presence of quasi-two-dimensional structures such as vortex sheets.
More generally, our present results reinforce the strong equivalence between classical and quantum turbulence, and constitute an attempt at providing an explicit connection between the intermittency of both systems.We expect such a connection to provide a possible path to a simplified description of the intermittency of classical turbulence, a highly challenging topic from a modelling standpoint, yet extremely relevant to the understanding of fluid flows occurring in Nature.

Methods
Numerical simulations.We study the dynamics of quantum turbulence in the framework of a generalised GP model where ψ is the condensate wave function describing the dynamics of a compressible superfluid at zero temperature.Here,  is the mass of the bosons, μ is the chemical potential,  0 the particle density and  = 4πℏ 2   / is the coupling constant proportional to the -wave scattering length.The dimensionless parameters χ and γ correspond to the amplitude and order of beyond mean field corrections.The non-local interaction between bosons is given by the potential V I (x − y) which is chosen, together with χ and γ, to reproduce the roton minimum in the excitation spectrum and the equation of state of superfluid helium.Details on the chosen parameters can be found in Ref. 22.The use of a standard or a generalised GP model does not affect the statistics of velocity circulation 13 .The hydrodynamic interpretation of Eq. ( 9) stems from the Madelung transformation ψ = √︁ ρ/ ϕ/ℏ , where ρ is the local density and ϕ the phase of the complex wave function.The velocity field is then given by v = ∇ϕ.Note that ϕ is not defined at the locations where ψ vanishes, which implies that the velocity field is singular along quantum vortices 35 .
The generalised GP equation ( 9) is solved in a three-dimensional periodic cube by direct numerical simulations using the Fourier pseudospectral code FROST, with an explicit fourth-order Runge-Kutta method for the time integration 22 .The quantum turbulent regime is studied in a freely decaying Arnold-Beltrami-Childress (ABC) flow 13,21 with 2048 3 collocation points.To reduce acoustic emissions, the initial condition is prepared using a minimisation process 20 .The box has a size L = 1365ξ and the inter-vortex distance is ℓ ≈ 28ξ, with ξ the healing length.
We also perform direct numerical simulations of the incompressible Navier-Stokes equations v using the Fourier pseudospectral code LaTu 64 in a periodic cubic domain.
The temporal advancement is performed with a third-order Runge-Kutta scheme.Above,  is the pressure field, ν the fluid kinematic viscosity, and f an external forcing stirring the fluid.The latter acts at large scales within a spherical shell of radius |k | ≤ 2 in Fourier space.The turbulent regime is studied once the simulation reaches a statistically steady state.The simulation is performed using 2048 3 collocation points at a Taylor-based Reynolds number of R λ = 510.

Evaluation of circulation and coarse-grained dissipation.
To obtain the circulation from GP and NS simulation data, we take advantage of the spectral nature of both solvers, and compute the circulation from the Fourier coefficients of the velocity fields.Namely, over a given L-periodic 2D cut of the physical domain, we write the circulation over a square loop of side  , centred at a point x = ( , ), as the convolution where ω = (∇ 2D × v) • ẑ is the out-of-plane vorticity field and B  (x) is a square of side  centred at x.The convolution kernel can be written as the product of two rectangular functions, H  (x) = Π( / ) Π( / ), where Π( ) = 1 for | | < 1 2 and 0 otherwise.Note that we have used Stokes' theorem to recast the contour integral (1) as a surface integral of vorticity.The convolution in Eq. ( 12) can be efficiently computed in Fourier space using the Fourier transform of the rectangular kernel, which may be written in terms of the normalised sinc function as H  (  ,   ) = ( /L) 2 sinc(   /2π) sinc(   /2π).
As mentioned earlier, the GP velocity field diverges at vortex locations.To minimise the numerical errors resulting from such singularities, we first resample each two-dimensional cut of the GP wave function field ψ(x) into a very fine grid of resolution 32768 2 , using Fourier interpolation.The velocity field is then evaluated in physical space using the Madelung transformation.This resampling procedure is described in more detail in Ref. 13.
In NS simulations, the above algorithm is also applied to compute the coarse-grained dissipation ε  (x) over squares of side  .Instead of the vorticity, the convoluted quantity is in this case the dissipation field ε(x) = 2ν      , where s(x) = [∇v + (∇v) T ]/2 is the three-dimensional strain-rate tensor.
Vortex detection from GP simulations.For a given two-dimensional cut of a GP velocity field, we identify the signs and locations of the quantum vortices crossing the cut as follows.First, the circulation is computed on a discrete grid following the procedure described above, taking small square loops of side  ∼ ξ ℓ.The result is a discrete circulation field, where each circulation value is either zero if no vortex crosses the small loop centred at that position, or ±κ if a single vortex crosses it.For very small loop sizes, the former case is much more likely than the latter.As a result, the vortex distribution can be sparsely described by storing the locations and signs of the non-zero circulation values.By repeating this procedure over different cuts of the simulation, one can reconstruct the three-dimensional vortex structure, as visualised in Fig. 1.
Energy spectrum computation from discrete vortices.For each twodimensional cut, once the positions r  and the signs   of each vortex crossing the plane are determined, we first compute a regularised two-dimensional vorticity field ω(r) = κ N =1   δ η (r − r  ), where N is the number of vortices on the 2D cut.Here, δ η (r) = exp (− |r | 2 /2η 2 )/2πη 2 , and η is the scale of the regularisation (we have used η = ξ in Fig. 2).Then, the energy spectra are computed by noting that where v and ω are the Fourier transforms at the wavevector k of the velocity field and of ω, respectively.Finally, by averaging over all 2D cuts and integrating over a shell |k| = , the energy spectrum reads where the integral is performed over all angles Ω.Note that the largewavenumber range in Fig. 2 is determined by the regularised Dirac function δ η and has no physical meaning.For disordered turbulence, as there is no correlation between the signs and the vortex positions, it is easy to show that ,       k•(r  −r  ) = N , from where it follows E() ∼  −1 .
with n  = (r − r  )/Δ  .By replacing these expressions in Eq. (B5), we obtain the local energy dissipation of a point-vortex system where B(r  , ξ) is a small ball of radius ξ around vortex r  and we omitted primes on the integration variables to simplify notation.Those balls are excluded to avoid the divergences of point vortices and it is justified by the regularisation of the vorticity.As we will see, the integrals are dominated by such divergences.
For  = , the integral is simpler and becomes where we have used that  2  +  2  = 1, assumed that  ξ and kept only the dominant contribution.For  ≠ , the divergence is milder as in the denominator it intervenes the distance between two vortices    , which is assumed to be much larger than ξ for a diluted system.Such terms contribute with a divergence of order log (ξ)/(   ) 2 .
Finally, using Eq.(B9) and keeping only dominant terms, we obtain at distances much larger than the inter-vortex distance Note that in terms of scaling laws, we could have used Eq.(B3) and ε  ∼ ν eff Ω  to obtain the same result.

C. Suitability of OK62 theory for circulation scaling exponents
The celebrated K41 theory for turbulence describes a self-similar behavior for the scaling exponents of the high-order moments of the velocity increments 17 .This theory can also be adapted to describe the scaling exponents of velocity circulation |Γ|  ∼  λ  , resulting in the scaling λ K41  = 4/3.However, it was observed that the self-similar hypothesis breaks down generating deviations in the scaling exponents, in particular for high-order moments.It was later proposed by Obukhov and Kolmogorov in 1962 a refined similarity hypothesis 37,38 , that can also be applied to the circulation exponents as where τ( ) corresponds to the scaling exponents of the energy dissipation ϵ   ∼  τ( ) .There are different models that describe the behavior of τ( ).In particular, in this work we use the She-Lévêque model 45 that has no adjustable parameters.Figure S1 shows the scaling exponents λ  of the velocity circulation for moments up to order 8 obtained from numerical simulations of classical and quantum turbulence.Both numerical simulations were performed using 2048 3 collocation points, with a Re λ = 510 in classical turbulence and a scale separation of L/ξ = 1365 in quantum turbulence.For high order moments, the scaling exponents deviate from both K41 and the refined OK62 models.These deviation were also observed in low-Reynolds numbers numerical simulations of the Navier-Stokes equation using the log-normal model for the scaling of dissipation 49 .

Fig. 1
Fig. 1 Visualisation of a quantum turbulent vortex tangle.Instantaneous state obtained from GP simulations using 2048 3 collocation points.Green and yellow colours correspond to opposite orientations of the vortex lines with respect to the vertical direction.The inset shows a horizontal two-dimensional cut of the system.See Methods for the vortex identification algorithm.

Fig. 2
Fig. 2 Kinetic energy spectrum in quantum and classical turbulence.Spectra are obtained from simulations of the generalised Gross-Pitaevskii (GP) model and the incompressible Navier-Stokes (NS) equations.Wave numbers  are respectively normalised by the mean inter-vortex distance ℓ and by the Taylor micro-scale λ, while

Fig. 3
Fig. 3 Circulation statistics in quantum turbulence.(a) Circulation variance as a function of the area A of each loop and of the number  of neighbouring vortices.The dotted line represents the scaling A 1 observed at small scales of quantum turbulence 13 .(b-c) Moments and local slopes of |Γ|   as a function of Γ 2  according to the ESS approach, for  ∈ [2, 8].Dashed green lines represent the K41 scaling γ  = 2/3.(d-e) PDFs of the circulation for different (d) numbers of vortices and (e) loop areas.All PDFs are normalised by the respective standard deviations.Dashed red lines represent a unit Gaussian distribution.

Fig. 4
Fig. 4 Velocity circulation scaling exponents.Exponents of the turbulent (in blue) and the disordered turbulence (in red) cases.For each case, the scaling exponents are defined as |Γ|  A ∼  λ  (circles) and |Γ|   ∼  γ  (stars).Error bars indicate 95% confidence intervals.Self-similar predictions for each case are shown as solid lines.Green squares show the relation between the turbulent and disordered systems given by Eq. (5).The blue and red dashed lines show the OK62 prediction combined with the She-Lévêque model (8) with D ∞ = 2.2 (termed "mOK62") for turbulence and disordered turbulence, respectively.
-the normalised PDFs of log(Z A ) almost collapse and are close to Gaussian in the bulk [Fig.5(c)],reinforcing the pertinence of relation(6).

Fig. 5
Fig. 5 Statistics of spatial vortex distribution in quantum turbulence.(a) Moments of the number of vortices normalised by the area that contains them, Z A = /A.Dashed lines correspond to the scaling of Eq. (7), using the She-Lévêque prediction 45 for the anomalous exponents τ( ) with D ∞ = 1.(b) PDFs of the number of vortices contained in loops of varying area A/ℓ 2 .The dashed line corresponds to a Poisson distribution of mean equal to  A at A = 140ℓ 2 .(c) Centred reduced PDFs of log(Z A ) for different values of A/ℓ 2 .A log-normal distribution is shown as reference (red dashed line).

Fig. 6
Fig. 6 Circulation intermittency and OK62 theory in classical and quantum turbulence.Top panels: PDFs of the circulation in classical turbulence as a function of the loop size .(a) Unconditioned PDFs.(b) PDFs conditioned on low values of the local coarse-grained dissipation, ε  / ε ∈ [0.5, 1].The different colours correspond to different loop sizes within the inertial range.(c-d) Scaling exponents of the circulation moments in (c) classical and (d) quantum turbulence.Different colours indicate a conditioning (c) on the local dissipation and (d) on the number of vortices within each loop.The unconditioned exponents are shown with black crosses.The Kolmogorov self-similar scaling is shown as reference (red dashed line).Error bars indicate 95% confidence intervals.
grained energy dissipation ε  is defined by averaging ε(x) on a disk B of radius  containing all the vortices  | 2 |r − r  | 2 d 2 r, (B9) Fig. S1 Scaling exponents of the circulation moments defined as |Γ|  ∼  λ  in numerical simulations of classical and quantum turbulence.As reference, the solid red line displays the Kolmogorov 1941 scaling λ K41  = 4/3 and the dashed black line shows the OK62 scaling defined in (C1) using the She-Lévêque model for the scaling of dissipation (C2).Error bars indicate 95% confidence intervals.