Intermittency in the not-so-smooth elastic turbulence

Elastic turbulence is the chaotic fluid motion resulting from elastic instabilities due to the addition of polymers in small concentrations at very small Reynolds (\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${{{{{{{\rm{Re}}}}}}}}$$\end{document}Re) numbers. Our direct numerical simulations show that elastic turbulence, though a low \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${{{{{{{\rm{Re}}}}}}}}$$\end{document}Re phenomenon, has more in common with classical, Newtonian turbulence than previously thought. In particular, we find power-law spectra for kinetic energy E(k) ~ k−4 and polymeric energy Ep(k) ~ k−3/2, independent of the Deborah (De) number. This is further supported by calculation of scale-by-scale energy budget which shows a balance between the viscous term and the polymeric term in the momentum equation. In real space, as expected, the velocity field is smooth, i.e., the velocity difference across a length scale r, δu ~ r but, crucially, with a non-trivial sub-leading contribution r3/2 which we extract by using the second difference of velocity. The structure functions of second difference of velocity up to order 6 show clear evidence of intermittency/multifractality. We provide additional evidence in support of this intermittent nature by calculating moments of rate of dissipation of kinetic energy averaged over a ball of radius r, εr, from which we compute the multifractal spectrum.


I. INTRODUCTION
Turbulence is a state of irregular, chaotic and unpredictable fluid motion at very high Reynolds numbers (Re), which is the ratio of typical inertial forces over typical viscous forces in a fluid.It remains one of the last unsolved problems in classical physics.Conceptually, the fundamental problem of turbulence shows up in the simplest setting of statistically stationary, homogeneous and isotropic turbulent (HIT) flows: What are the statistical properties of velocity fluctuations?More precisely, consider the (longitudinal) structure function of velocity difference across a length-scale r : x O e p M h s 7 P L z K w Q l v y B F w + K e P W P v P k 3 T p I 9 a L S g o a j q p r s r S A T X x n W / n M L K 6 t r 6 R n G z t L W 9 s 7 t X 3 j 9 o 6 T h V D J s s F r H q B F S j 4 B K b h h u B n U Q h j Q K B 7 W B 8 O / P b j 6 g 0 j + W D m S T o R 3 Q o e c g Z N V Z q N L B f r r h V d w 7 y l 3 g 5 q U C O e r / 8 2 R v E L I 1 Q G i a o 1 l 3 P T Y y f U W U 4 E z g t 9 V K N C W V j O s S u p Z J G q P 1 s f u m U n F h l Q M J Y 2 Z K G z N W f E x m N t J 5 E g e 2 M q B n p Z W 8 m / u d 1 U x N e + x m X S W p Q s s W i M B X E x G T 2 N h l w h c y I i S W U K W 5 v J W O N M a U D W k f W 5 Z K G q L 2 0 + m j Y 3 J q l S 7 p R c q O N G S q / r 1 I a a j 1 K A z s Z k j N Q C 9 6 E / F f 7 2 k W M J 9 u e j d + y m W c G J R s F t 5 L B D E R m T R B u l w h M 2 J k C W W K 2 / 8 J G 1 B F m b F 9 5 W 0 x 3 m I N y 6 R + X v K u S p f V i 2 L 5 N q s o B 8 d w A m f g w T W U 4 R 4 q U A M G C C / w C m / O s / P u f D i f s 9 U V J 7 s 5 g j k 4 X 7 8 X N 5 V / < / l a t e x i t > O 7 s k j e S L P z o P z 4 r w 6 b + P R K W e y s 0 F + w H n / A i y n o w I = < / l a t e x i t > ⇠ k 2.3 < l a t e x i t s h a 1 _ b a s e 6 4 = " K e i z g 3 P y o q n h e 5 6 T 8 f 3 s / P u f D i f s 9 U V J 7 s 5 g j k 4 X 7 8 X N 5 V / < / l a t e x i t > k < l a t e x i t s h a 1 _ b a s e 6 4 = " K e i z g 3 P y o q n h e 5 6 T 8 f 3  Here, u(x) is the velocity field as a function of the coordinates x and the symbol ⟨•⟩ denotes averaging over the statistically stationary state of turbulence.Here and henceforth we use the Einstein summation convention, repeated indices are summed.The p-th order structure function S p is the p-th moment of the probability distribution function (PDF) of velocity differences -if we know S p for all p then we know the PDF.Typically, energy is injected into a turbulent flow at a large length scale L, while viscous effects are important at small length scales η, called the Kolmogorov scale, and dissipate away energy from the flow.In the intermediate range of scales S p (r) ∼ r ζp where scaling exponents ζ p are universal, i.e. they do not depend on how turbulence is generated.The dimensional arguments of Kolmogorov give ζ p = p/3, which also implies that the shell-integrated energy spectrum (distribution of kinetic energy across wavenumbers) E(k) ∼ k −5/3 , where k is the wavenumber.Experiments and direct numerical simulations (DNS) over last seventy years have now firmly established that the ζ p (p) is a nonlinear convex function -a phenomenon called multiscaling or intermittency.Even within the Kolmogorov theory, turbulence is non-Gaussian because the odd order structure functions (odd moments of the PDF of velocity differences) are not zero.Intermittency is not merely non-Gaussianity, it implies that not only a few small order moments but moments of all orders are important in determining the nature of the PDF.We often write ζ p = p/3 + δ p , where δ p are corrections due to intermittency.A systematic theory that allows us to calculate ζ p starting from the Navier-Stokes equation is the goal of turbulence research.
Turbulent flows, both in nature and industry, are often multiphase, i.e. they are laden with particles, may comprise of fluid mixtures, or contain additives such as polymers.Of these, polymeric flows are probably the most curious and intriguing: the addition of high molecular weight (about 10 7 ) polymers in 10-100 parts per million (ppm) concentration to a turbulent pipe flow reduces the friction factor (or the drag) up to 5-6 times (depending on concentration) [3][4][5].Evidently, this phenomena, called turbulent drag reduction (TDR), cannot be studied in homogeneous and isotropic turbulent flows; nevertheless, polymer laden homogeneous and isotropic turbulent (PHIT) flows have been extensively studied theoretically [6][7][8], numerically [2,[9][10][11][12][13][14][15][16][17][18][19][20][21][22][23][24], and experimentally [1,[25][26][27][28], to understand how the presence of polymers modifies turbulence, following the pioneering work by Lumley [29] and Tabor and de Gennes [30].The simplest way to capture the dynamics of polymers in flows is to model the polymers as two beads connected by an overdamped spring with a characteristic time scale τ p .A straightforward parameterization of the importance of elastic effects is the Deborah number De ≡ τ p /τ f , where τ f is some typical time scale of the flow.In turbulent flows, such a definition becomes ambiguous because turbulent flows do not have a unique time scale, rather we can associate an infinite number of time scales even with a single length scale [31][32][33].In such cases, a typical timescale used to define De is the large eddy turnover time of the flow, τ L [2].The phenomena of PHIT appear at high Reynolds and high Deborah numbers.
Research in polymeric flows turned into a novel direction when it was realized that even otherwise laminar flows may become unstable due to the instabilities driven by the elasticity of polymers [34,35].Even more dramatic is the phenomena of elastic turbulence (ET) [36], where polymeric flows at low Reynolds but high Deborah numbers are chaotic and mixing, with a shell-integrated kinetic energy spectrum E(k) ∼ k −ξ .It is still unclear whether this exponent is universal or not -experiments and DNS in two dimensions have obtained 3 ≤ ξ ≤ 4, and theory [37] sets a lower bound with ξ > 3. A three-dimensional DNS of decaying homogeneous, isotropic turbulence with polymers additives (modelled as discrete dumbbells) also revealed an exponent ξ ≈ 4 at late times (with a mild De dependence), when turbulence had sufficiently decayed and elastic stresses were dominant, likely marking the onset of ET [20].In summary, as shown in Fig. (1), HIT (in Newtonian turbulence) appears at large Re and zero (and small) De; PHIT appears at large Re and intermediate De number, while ET appears at small Reynolds and large Deborah numbers.
Recently, experiments [1] and DNS [2] revealed an intriguing aspect of PHIT: The energy spectrum showed not one but two scaling ranges, a Kolmogorov-like inertial range at moderate wave numbers and a second scaling range with E(k) ∼ k −2.3 resulting purely due to the elasticity of polymers [2].This is illustrated in the gray shaded region in Fig. (1).Even more surprising is the observation that both of these ranges have intermittency correction δ p which are the same.This hints that even at low Re, where elastic turbulence (ET) appears, intermittent behavior may exist.In this paper, based on large resolution DNS of polymeric flows at low Reynolds number, we show that this is indeed the case.

II. MODEL
We generate a statistically stationary, homogeneous, isotropic flow of a dilute polymer solution by the DNS of the Navier-Stokes equations coupled to the evolution of polymers described by the Oldroyd-B model: Here, u is the incompressible solvent velocity field, i.e. ∂ β u β = 0, p is the pressure, S is the rate-of-strain tensor with components S αβ ≡ (∂ α u β + ∂ β u α )/2, µ f and µ p are the fluid and polymer viscosities, ρ f is the density of the solvent fluid, τ p is the polymer relaxation time, and C is the polymer conformation tensor whose trace C γγ is the total endto-end squared length of the polymer.To maintain a stationary state, we inject energy into the flow using an Arnold-Beltrami-Childress (ABC) forcing, i.e., The injected energy is ultimately dissipated away by both the Newtonian solvent (ε f ) and polymers (ε p ).The total energy dissipation rate, ⟨ε T ⟩, is given by: where We show typical snapshots of the two energy dissipation rates on two dimensional slices of our three dimensional DNS in Fig. (2).Details on numerical schemes and simulations are discussed in the Methods section.The Newtonian ABC flow shows Lagrangian chaos in the sense that the trajectories of tracer particles advected by such a flow have sensitive dependence on initial condition [38].Hence we expect that a polymer advected by the flow will go through a coil-stretch transition for large enough τ p .The back reaction from such polymers may give rise to elastic turbulence.The energy spectra of the Newtonian flows (and those for our non-Newtonian flows with small τ p ) do not show any power-law range, and drops-off rapidly in wavenumber k, see Fig. (S1a) in the Supplementary Material.Beyond a certain value of τ p , the flow becomes chaotic, and the resulting flows with De ≳ 1 are able to sustain elastic turbulence.Henceforth we focus only on the flows that show ET.
Note an important difference between ET and usual Newtonian HIT.In the latter, the Kolmogorov length (or time) scale is defined as the scale where the inertial and viscous effects balance each other.Although we continue to use the same definition -Re λ and τ K are calculated from the Newtonian DNS -these scales lose their usual meaning because ET appears at small Re at scales where the inertial term is negligible.The η we obtain is, as expected, quite close to the scale of energy injection.Therefore, we use the box-size L which is also the scale of energy injection, as our characteristic length scale.

III. RESULTS
We present our results for three different Deborah number flows with De = 1, 3, and 9 and Taylor scale Reynolds number Re λ ≈ 40.Let us begin by looking at the (shell-integrated) fluid energy spectrum where û(m) is the Fourier transform of the velocity field u(x).We show the spectra for the three De numbers in Fig. (3a).The spectra E(k) show power-law scaling over almost two decades when plotted on a log-log scale.
Clearly, E(k) ∼ k −ξ with ξ = 4 independent of the Deborah number.Note that in DNS of decaying PHIT, ξ goes from 2.3 to 4 (and beyond as turbulence decayed) as time progressed [20].While ξ = 2.3 was recently confirmed for PHIT via both DNS [2] and experiments [1], we now show via DNS that ET is, in fact, a stationary state marked by ξ = 4 which is sustained by purely elastic effects, for a large enough polymer elasticity.We have verified, using representative DNSs, that the scaling exponents of ET remains the same if we use the FENE-P model for polymers or a different forcing scheme [39], that is not white-in-time.We have also checked that reducing the resolution to N 3 = 512 3 reproduces the same spectra.Finally, we have turned off the advective nonlinearity in (2a) and also obtained the same spectra, thereby confirming that the turbulence we obtain is purely due to elastic effects.We plot all these spectra in Fig. (S1) of the Supplementary Material.
We also define the energy spectrum associated with polymer degrees of freedom as where the matrix B with components B αγ is the (unique) positive symmetric square root of the matrix C, defined by [40,41].We obtain E p (k) ∼ k −χ with χ = 3/2 as shown in log-log plot of E p (k) in Fig. (3b).Note that the scaling range of E p (k) is somewhat smaller than that of E(k).In the statistically stationary state of ET the effect of the advective nonlinearity must be subdominant.Hence, at scales smaller than the scale of the external force, the viscous term in the momentum equation must balance the elastic contribution [37].Using a straightforward scaling argument, described in detail in the Supplementary Material, section A2 we obtain : which is satisfied by the values of ξ and χ we obtain.For further confirmation we calculate all the contributions to the scale-by-scale kinetic energy budget in Fourier space (see the Supplementary Material, section A1).As expected, the contribution from the advective term in (2a) is negligible.Earlier theoretical arguments [37] have suggested ξ > 3 which has also been observed in experiments [42][43][44] ξ ≈ 3.5 over less than a decade of scaling range.We obtain ξ ≈ 4 which satisfies the inequality and agrees with shell-model simulations [45].Earlier theoretical arguments [37,46] had also assumed the same balance in the momentum equation that we have, but in addition had assumed scale separation and a large scale alignment of polymers in analogy with magnetohydrodynamics, obtaining ξ = χ + 2, which is not satisfied by our DNS.Note further that experiments often obtain power-spectrum as a function of frequency and they can be compared with power-spectrum as a function of wavenumber (typically obtained by DNS) by using the Taylor "frozen-flow" hypothesis [47].In the absence of a mean flow and negligible contribution from the advective term it is not a priori obvious that the Taylor hypothesis should apply to ET.We have confirmed from our DNS that a frequency dependent power-spectrum obtained from time-series of velocity at a single Eulerian point also gives ξ = 4 (see the Supplementary Material, Fig.

(S1b)).
A. Second order structure function Next, we consider the second order structure function, S 2 (r), which is the inverse Fourier transform of E(k).This requires some care.As a background, let us first consider the case of HIT (Newtonian homogeneous and isotropic turbulence).Let us ignore intermittency and concentrate on scaling a-la Kolmogorov.The second order structure function is expected to have the following form Shaded region shows the standard deviation on the exponents computed from 18 snapshots.The two red + symbols mark the exponents ζ2 and ζ4 obtained with an alternate forcing scheme [39].That they lie well within error bars goes on to show that the results are independent of the large scale forcing.(f) Probability distribution function of δ 2 u(r) for four different values of r at De = 9.The distributions are non-Gaussian at small separations, while they become closer to a Gaussian (shown as a black dash-dotted curve) for large r.The corresponding cumulative distribution functions, computed using the rank-order method, are shown in the Supplementary Material, section C.
The range of scales L > r > η is the inertial range.Let us remind the reader that the behaviour S 2 ∼ r 2 for small enough r follows from the assumption that the velocities are analytic functions of coordinates, which must always hold for any finite viscosity, however small.We call S 2 ∼ r 2 the trivial scaling.The strategy to extract the exponent ζ 2 from DNS is to run simulations at higher and higher Reynolds number, which means smaller and smaller η to obtain a significant inertial range from which ζ 2 can be extracted.In Kolmogorov theory ζ 2 = 2/3, we call this the non-trivial scaling.The theory of ET is much less developed than that of HIT.Nevertheless, we may assume that the velocities must still be analytic functions.Hence for ET the following must hold As this scaling follows directly as a consequence of analyticity of velocities, we again call this the trivial scaling.If there is a non-trivial scaling of second order structure function in ET -here we are not talking about intermittency corrections to a non-trivial scaling but just the existence of the non-trivial scaling exponent -it may show up in the following manner )). Usual straightforward power counting implies that S 2 (r) ∼ r ξ−1 ∼ r 3 (see section A2 of Supplementary Material for details), while we obtain S 2 ∼ r 2 .This paradox is resolved by noting that, in the limit r → 0, r 3 is subdominant to r 2 , hence S 2 (r) ∼ r 2 as r → 0 for any velocity field whose spectra E(k) ∼ k −ξ with ξ > 3, see e.g., Ref. [48,Appendix G].This is also known from the direct cascade regime of two-dimensional turbulence (with Ekman friction) where E(k) ∼ k −γ with γ > 3 and S 2 (r) ∼ r 2 , see e.g., Refs.[49, 50, page 432].This suggests that S 2 (r) satisfies ( 9) with ζ 2 = ξ − 1 ≈ 3. To test this we plot the compensated second order structure function S 2 (r)/r 3 as a function of r in the Supplementary Material Fig. (S3c).We detect no range at small or intermediate r where this non-trivial scaling holds.Now we consider the possibility that S 2 (r) ∼ Ar 2 + Br 3 + h.o.t.Here, the symbol h.o.t denotes higher order terms in r.We use a trick [51,52] to extract the subleading term which scales with the non-trivial scaling exponent: the idea is to remove the analytic contribution by considering the second difference of velocities: We plot Σ 2 (r) in Fig. (3c).We find that Σ 2 (r) shows a significant scaling range as r → 0 with the non-trivial scaling exponent ζ 2 ≈ 3.This implies that, in ET the velocity fluctuations across a length scale r can be expanded in an asymptotic series in r as where G αβ and H αβ are (undetermined) expansion coefficients, and h ≈ ζ 2 /2 = 3/2.The use of Σ 2 is necessary to extract the subleading contribution.
To appreciate the importance of this result, let us revisit the Kolmogorov theory of turbulence: in the limit r → 0 at a finite viscosity µ f , ⟨δu α (r)⟩ ∼ r since velocity gradients are finite.But if we first take the limit ν → 0 and then r → 0 (ν ≡ µ f /ρ f is the kinematic viscosity) ⟨δu α (r)⟩ ∼ r h with h ≈ 1/3.The velocity field is rough.In contrast, ET is, by definition, a phenomenon at a finite viscosity (small Reynolds number), thus, the limit ν → 0 does not make sense -the velocity field is always smooth.But the non-trivial nature of ET manifests itself in the first subleading term in the expansion (11), and this is best revealed not by the velocity differences, but by the second difference of velocity.This is the first important result of our work.

B. Intermittency based on velocity differences
The crucial lesson to learn from the previous section is that in ET, to uncover the non-trivial scaling of velocity differences we must use the second differences of velocity rather than the usual first difference.Other than this peculiarity, the rest of this section follows the standard techniques [47] used to study intermittency/multifractality.
We define the p-th order structure function of the second difference of velocity across a length scale r as: We show a representative plot of Σ p 's for all integer p = 1, ... To obtain reasonable error bars on ζ p , we have proceeded in the following manner: first we find a suitable scaling regime for each order by visual inspection; next, in these chosen ranges, we find the local slopes of the log-log plot of Σ p (r) vs r, to obtain ζ p as a function of r: ζ p (r) = (∆ log Σ p (r))/(∆ log r).This process is repeated for multiple time snapshots (two successive snapshots are separated by at least one eddy turnover time) of the velocity field data.The mean value over the set of exponents thus obtained is the exponent ζ p in Fig. (3e) and the standard deviation sets the error bar which are shown as a shaded region.Clearly, ζ p is a non-linear function of p.This unambiguously establishes the existence of intermittency in ET.Furthermore, we have confirmed that the structure functions S p of even order up to 6 grow as r p for small r, see the Supplementary Material section B3 for discussion.Thereby we confirm, following the prescription in Ref. [53], that the structure functions of all order are analytic.The structure functions begin to depart from this analytic scaling at a scale that depends very weakly on De (if at all) but this scale decreases as p increases.We also use another forcing scheme [39] and calculate the exponents ζ p for p = 2 and 4, marked as two + symbols in red colour in Fig. (3e).Within errorbars they agree with the values we have obtained suggesting that the ζ p are universal.
Let us again emphasize that intermittency is a fundamental property of structure functions, both S p and Σ p .The use of Σ p is merely to help us extract the exponents ζ p .

PDF of velocity differences
Another way to demonstrate the effects of intermittency is by looking at the PDF of velocity differences across a length scale.From the structure function we have obtained intermittent behavior for scales r/L < 1.Thus, we expect the PDF of velocity differences to be close to Gaussian for r/L ≈ 1 and to have long tails (decaying slower than Gaussian) for r/L < 1.This indeed is the case, as is shown in Fig. (3f) where we plot the PDFs of δ 2 u(r) for different separations r (for De = 9).The tails of the distribution of δ 2 u decay much slower than Gaussian, thereby clearly demonstrating intermittency.
Note that also the PDF of the usual velocity differences δu(r) is non-Gaussian, see the Supplementary Material Fig. (S6).But the PDF of δ 2 u falls off slower than δu at the same r, in other words larger fluctuations are more likely to appear in the second difference of velocity -it is more intermittent.This non-Gaussianity of probability (c) The multifractal spectra of the fluid dissipation field calculated from the scaling of the energy dissipation rate εr calculated over a cube of side r.The black dash-dotted line shows the spectrum for Newtonian HIT [56].
distributions can be quantified by the Kurtosis (also called Flatness) defined by for the first and second difference of velocity, respectively.For Gaussian distributions the Kurtosis is 3.We find K 2 ≈ 3 as r → L, i.e., the PDFs (of δ 2 u) are close to Gaussian for large separations.From the scaling behaviour of structure function we obtain K 2 (r) ∼ r ζ4−2ζ2 as r → 0, which is consistent with K 2 (r) ∼ r −1.63 obtained from the distribution of second differences shown in Fig. (4a).Furthermore, we find that the Kurtosis is independent of De.The gray shaded region marks the range used to compute the scaling exponents for K 2 (r).We obtain: −1.6 ± 0.3, −1.6 ± 0.1, and −1.6 ± 0.1 for De = 1, 3 and 9 respectively.This is further evidence in support of universality of intermittency in ET.The Kurtosis of first difference K(r) is also close to a Gaussian as r → L, but grows much slower than K 2 (r) as r → 0 as shown in Fig. (4b).In Newtonian HIT at high Re, the most recent DNS [54,55] shows ζ 4 − 2ζ 2 ≈ −0.12 so that K ∼ r −0.12 .Hence, the intermittency we obtain in ET is more intense than what is observed in HIT.We also calculate the cumulative PDF (CDF) of δ 2 u by rank-order method, thereby avoiding the usual binning errors that appear while calculating PDFs via histograms.In section C1 of the Supplementary Material we show that rescaling the abscissa of the CDFs by the root-mean-square value of δ 2 u does not collapse the CDFs for different r, i.e., the PDFs are not Gaussian.
Altogether the PDFs of δ 2 u provide us three additional evidences in support of intermittency in ET.

C. Intermittency based on dissipation
In HIT (high Re, Newtonian turbulence) there are two routes to study intermittency: one is through structure functions and another is through the fluctuations of the energy dissipation rate [47] -the PDF of ε f deviates strongly from a log-normal behaviour [57].We now take the second route for ET, in which case there are two contributions to the total energy dissipation -ε f and ε p .In Fig. (5a) and Fig. (5b) we plot the PDFs of the logarithm of ε f and ε p , respectively.We find that the former decays slower than a Gaussian, i.e., the PDF itself falls off slower than a log-normal, whereas the latter decays faster than a Gaussian.The fact that the PDF of ε p falls off much faster than that of ε f can even been seen by comparing Henceforth, following the standard analysis pioneered by Meneveau and Sreenivasan [56] for HIT, we study the scaling of the q-th moment of the viscous dissipation averaged over a cube of side r, ⟨ε q r ⟩ ∼ r λq , where (14a) Here the symbol ⟨•⟩ r denotes averaging over a cube of side r.The Legendre transform of the function λ q gives the multifractal spectrum (also called the Cramer's function) f (α): where singularities in the dissipation field with exponent α − 1 lie on sets of dimension α.We plot the f (α) spectrum for ET in Fig. (5c).There are minor differences between the multifractal spectrum for De = 3 and 9 on one hand and De = 1 on the other hand.The clear collapse of the multifractal spectra at large De hints towards a universal multifractality in ET in the limit of large De.For comparison, we also plot, in Fig. (5c) the multifractal spectrum for HIT as a black dash-dotted curve [56].In HIT the intermittency model based on velocity are closely connected to the intermittency models based on dissipation [47].The development of such a formalism for ET, although important, is not considered in this work.

IV. DISCUSSION AND SUMMARY
We note that the phenomenon of elastic turbulence has no Newtonian counterpart -in the absence of the polymers this phenomenon disappears.Nevertheless, as HIT is the model of turbulence that has been studied in great detail we have used it as an illustrative example to compare with ET.Such comparison must be done with care.In HIT, the theory of Kolmogorov helps us understand the simple scaling of the energy spectrum, although a systematic derivation starting from the Navier-Stokes equation is still lacking.The key insight of Kolmogorov theory is that the energy flux across scales, due to the nonlinear advective term, is a constant.In practice, the flux is a fluctuating quantity, where its mean value determines the simple scaling prediction ζ p = p/3, while the fluctuations of the flux is the reason behind intermittency.The fluctuations of the flux shows up as fluctuations of the energy dissipation rate (because the advective term conserves energy) which is multifractal.
Elastic turbulence was first discovered at the start of this century.Almost all studies of ET, so far, have concentrated on understanding the scaling of the energy spectrum.A theory at the level of Kolmogorov theory for HIT is still lacking.Nevertheless, it is clear that the mechanism of ET is very different from HIT.In the latter, it is the nonlinear advective term that is responsible for turbulence, while in the former the advective term is expected to be subdominant, it is the stress from the polymers that must balance the viscous dissipation in the range of scales where ET is found.We show that this is indeed the case.A consequence of this balance is that the scaling exponents of E(k) and E p (k) are related to each other by (6).In ET there is not one but two possible mechanisms of energy dissipation.A particularly intriguing result we obtain is that only one of them, ε f , shows intermittent behaviour, since the energy dissipation rate due to the polymers is not intermittent, and its logarithm remains sub-Gaussian.
In summary, we have shown that both the velocity field and the energy dissipation field in ET are intermittent/multifractal.But this multifractality is very different from the multifractality seen in HIT.In HIT, in the limit of viscosity going to zero, the velocity field is rough.In contrast, the velocity field in ET is smooth at leading order, and roughness and the multifractal behavior appears due to the sub-leading term.Consequently, although the velocity difference across a length scale is intermittent, it is necessary to use the second difference of velocity, to properly reveal the intermittency.Finally, note that in HIT, the multifractal exponents are expected to be universal, i.e., they are independent of the method of stirring and the Reynolds number (in the limit of large Reynolds number).In ET the multifractality appears at small Re and large De > 1.All the evidences from our DNSs suggest that intermittency in ET is also universal with respect to Deborah number, method of stirring and choice of model of polymers, although significant future work with high resolution DNSs are necessary to provide conclusive evidence.

Methods
We solve eqns.2a, 2b using a second order central-difference scheme on a L = 2π tri-periodic box discretized by N 3 = 1024 3 collocation points, such that L/N = ∆ ≈ 0.05η, where η ≡ (ν 3 / ⟨ε f ⟩) 1/4 is the Kolmogorov dissipation length scale and ν ≡ µ f /ρ is the kinematic viscosity.Integration in time is performed using the second order Adams-Bashforth scheme with a time step ∆t ≈ 10 −5 τ K , with τ K ≡ ν/ ⟨ε f ⟩.We use 18 snapshots in our analysis, with successive snapshots separated by ≈ 1.6 × 10 4 τ K .We choose a µ f so as to obtain a laminar flow in the Newtonian case with A = B = C = 1 -this corresponds to the Taylor scale Reynolds number Re λ ≈ 40.Next, we choose µ p such that the viscosity ratio µ f /(µ f + µ p ) = 0.9 -this corresponds to dilute polymer solutions [16], and we vary τ p over two order of magnitudes.The numerical solver is implemented on the in-house code Fujin; see https: //groups.oist.jp/cffu/codefor additional details and validation tests.The very same code has been successfully used on various problems involving Newtonian and non-Newtonian fluids [58][59][60][61].

S
p (r) ≡ ⟨[δu(r)] p ⟩ , (1a) where δu(r) ≡ [u α (x + r) − u α (x)] r α |r| .(1b) < l a t e x i t s h a 1 _ b a s e 6 4 = " Q k 5 F W m 9 2 5 n 9 Z b O I g p a H X H X p 2 l H c = " > A A A B 6 X i c b V D L S g N B E O y N r x h f U Y 9 e B o P g K e y K r 2 P Q i 8 c Y z A O S J c x E F W X G h l O y I X j L L / 8 l r b O q d 1 m 9 u D + v 1 G 7 y O I p w B M d w C h 5 c Q Q 3 u o A 5 N Y B D C E 7 z A q z N 2 n p 0 3 5 3 3 R W n D y m U P 4 B e f j G 2 4 6 j U 8 = < / l a t e x i t > Re < l a t e x i t s h a 1 _ b a s e 6 4 = " u V R L o d K H C v e x Q d N n I / c d z x z t Z a 0 = " > A A A B 6 X i c b V D L S g N B E O y N r x h f U Y 9 e B o P g K e y K r 2 N Q D x 6 j m A c k S 5 i d 9 C Z D Z m e X m V k h L P k D L x 4 U 8 e o f e fN v n C R 7 0 G h B Q 1 H V T X d X k A i u j e t + O Y W l 5 Z X V t e J 6 a W N z a 3 u n v L v X 1 H G q G D Z Y L G L V D q h G w S U 2 D D c C 2 4 l C G g U C W 8 H o e u q 3 H l F p H s s H M 0 7 Q j + h A 8 p A z a q x 0 f 4 O 9 c s W t u j O Q v 8 T L S Q V y 1 H v l z 2 4 / Z m m E 0 j B B t e 5 4 b m L 8 j C r D m c B J q Z t q T C g b 0 Q F 2 L J U 0 Q u 1 n s 0 s n 5 M g q f R L G y p Y 0 Z K b + n M h o p P U 4 C m x n R M 1 Q L 3 p T 8 T + v k 5 r w 0 s + 4 T F K D k s 0 X h a k g J i b T t 0 m f K 2 R G j C 2 h T H F 7 K 2 F D q i g z N p y S D c F b f P k v a Z 5 U v f P q 2 d 1 p p X a V x 1 G E A z i E Y / D g A m p w C 3 V o A I M Q n u A F X p 2 R8 + y 8 O e / z 1 o K T z + z D L z g f 3 1 j 0 j U E = < / l a t e x i t > De Laminar flow Elastic turbulence Polymeric turbulence < l a t e x i t s h a 1 _ b a s e 6 4 = " 2 7 Z B A h y y K d c 4 m R z T d W Q m Y M a M v m 4 = " > A A A B / H i c b V D L S g N B E O z 1 G e M r 6 t H L Y B A 8 h V 3 x d Q x 6 8 Z i A e U C y h N l J J x k y O 7 v M z I p h i T / g V f / A m 3 j 1 X / w B v 8 N J s g e T W N B Q V H V T T Q W x 4 N q 4 7 r e z s r q 2 v r G Z 2 8 p v 7 + z u 7 R c O D u s 6 S h T D G o t E p J o B 1 S i 4 x J r h R m A z V k j D Q G A j G N 5 N / M Y j K s 0 j + W B G M f o h 7 U v e 4 4 w a K 1 W H n U L R L b l T k G X i Z a Q I G S q d w k + 7 G 7 E k R G m Y o F q 3 P D c 2 f k q V 4 U z g O N 9

k
< l a t e x i t s h a 1 _ b a s e 6 4 = " m R l Q Q Z g J t t 8 v C d j B + z Z t b F D T / o 8 = " > A A A C H H i c b V D L S g M x F M 3 U d 3 1 V X b o w W A Q 3 l h m p j 6 X o x q W C V c E Z S y a 9 b U M z m S G 5 I y 3 D L P 0 N f 8 C t / o E 7 c S v 4 A 3 6 H 6 W N h q w c C h 3 P u z b m c M J H C o O t + O Y W p 6 Z n Z u f m F 4 u L S 8 s p q a W 3 9 2 s S p 5 l D j s Y z 1 b c g M S K G g h g I l 3 C Y a W B R K u A k 7 Z 3 3 / 5 g G 0 E b G 6 w l 4 C Q c R a S j Q F Z 2 i l e m n L R + j i 4 J 8 s 1 k y 1 I M 9 8 I y L a u c / 2 q n l e L 5 X d i j s A / U u 8 E S m T E S 7 q p W + / E f M 0 A o V c M m P u P D f B I G M a B Z e Q F / 3 U Q M J 4 h 7 X g z l L F I j B B N j g g p z t W a d B m r O 1 T S A f q 7 4 2 M R c b 0 o t B O R g z b Z t L r i / 9 6 3 W H A e D o 2 j 4 N M q C R F U H w Y 3 k w l x Z j 2 m 6 I N o Y G j 7 F n C u B b 2 f s r b T D O O t s + i L c a b r O E v u d 6 v e I e V g 8 t q + e R 0 V N E 8 2 S T b Z J d 4 5 I i c k H N y Q W q E k 0 f y T F 7 I q / P k v D n v z s d w t O C M d j b I G J z P H w / 9 o w Q = < / l a t e x i t > ⇠ k 4 < l a t e x i t s h a 1 _ b a s e 6 4 = " / h a 4 2 J y 4 P f w t o 3 y y L k G N n P B z X T o = " > A A A B / H i c b V D L S g N B E O y N r x h f U Y 9 e B o P g K e y K r 2 N Q B I 8 J m A c k S 5 i d 9 C Z D Z h / M z I p h i T / g V f / A m 3 j 1 X / w B v 8 N J s g e T W N B Q V H V T T X m x 4 E r b 9 r e V W 1 l d W 9 / I b x a 2 t n d 2 9 4 r 7 B w 0 7 J e d o 3 L l 6 r B 0 e p Z V l C c b Z I v s E I c c k 1 N y Q S 5 J l X D y S J 7 J C 3 m 1 n q w 3 6 9 3 6 G I 3 m r G x n n Y z B + v w B W W O i k w = = < / l a t e x i t > ⇠ k 5/3 < l a t e x i t s h a 1 _ b a s e 6 4 = " / h a 4 2 J y 4 P f w t o 3 y y L k G N n P B z X T o = " > A A A B / H i c b V D L S g N B E O y N r x h f U Y 9 e B o P g K e y K r 2 N Q B I 8 J m A c k S 5 i d 9 C Z D Z h / M z I p h i T / g V f / A m 3 j 1 X / w B v 8 N J s g e T W N B Q V H V T T X m x 4 E r b 9 r e V W 1 l d W 9 / I b x a 2 t n d 2 9 4 r 7 B w 0 r e J 3 o x 3 w y A M f u W R a 1 x 0 7 x G b M F A o u I S k 0 I g 0 h 4 z 3 W h b q h P v N A N + P 0 g I R u G 6 V N O 4 E y z 0 e a q r 8 3 Y u Z p P f B c M + k x v N e T 3 l D 8 1 + u P A s b T s X P S j I U f R g g + H 4 V 3 I k k x o M O m a F s o 4 C g H h j C u h L m f 8 n u m G E f T Z 8 E U 4 0 z W 8 J f c 7 J e d o 3 L l 6 r B 0 e p Z V l C c b Z I v s E I c c k 1 N y Q S 5 J l X D y S J 7 J C 3 m 1 n q w 3 6 9 3 6 G I 3 m r G x n n Y z B + v w B W W O i k w = = < / l a t e x i t > ⇠ k 5/3 < l a t e x i t s h a 1 _ b a s e 6 4 = " / h a 4 2 J y 4 P f w t o 3 y y L k G N n P B z X T o = " > A A A B / H i c b V D L S g N B E O y N r x h f U Y 9 e B o P g K e y K r 2 N Q B I 8 J m A c k S 5 i d 9 C Z D Z h / M z I p h i T / g V f / A m 3 j 1 X / w B v 8 N J s g e T W N B Q V H V T T X m x 4 E r b 9 r e V W 1 l d W 9 / I b x a 2 t n d 2 9 4 r 7 B w 0 V J Z J h n U U i k i 2 P K h Q 8 x L r m W m A r l k g D T 2 D T G 9 5 O / O Y j S s W j 8 E G P Y n Q D 2 g + 5 z x n V R q r d d Y s l u 2 x P Q Z a J k 5 E S Z K h 2 i z + d X s S S A E P N B F W q 7 d i x d l M q N W c C x 4 V O o j C m b E j 7 2 D Y 0 p A E q N 5 0 + O i Y n R u k R P 5 J m Q k 2 m 6 t + L l A Z K j Q L P b A Z U D 9 S i N x H / 9 Z 5 m A f P p 2 r 9 2 U x 7 G i c a Q z c L 9 R B A d k U k T p M c l M i 1 G h l A m u f m f s A G V l G n T V 8 E U 4 y z W s E w a Z 2 X n s n x R O y 9 V b r K K 8 n A E x 3 A K D l x B B e 6 h C n V g g P A C r / B m P V v v 1 o f 1 O V v N W d n N I c z B + v o F 2 p i V W Q = = < / l a t e x i t > E < l a t e x i t s h a 1 _ b a s e 6 4 = " 2 7 Z B A h y y K d c 4

FIG. 1 .
FIG. 1. Polymeric flows.An illustrative sketch of the different regimes of polymer-laden flows: classical Newtonian turbulence (HIT) at large Re and zero (or small) De, polymeric turbulence (PHIT) at large Re and intermediate De, and elastic turbulence (ET) at small Re and large De.The shaded region shows the recently observed elastic scaling regime in addition to the classical Kolmogorov scaling [1, 2].

FIG. 2 .
FIG. 2. Flow visualizations.Two-dimensional slices of the three-dimensional domain showing snapshots of the normalized (a) fluid dissipation field ε f / ⟨ε f ⟩ and of (b) the polymer dissipation field εp/ ⟨εp⟩ in ET for De = 9.

3 FIG. 3 .
FIG. 3. Spectra and structure functions.(a) The fluid energy spectra show a universal scaling E(k) ∼ k −4 independent of De.A steeper than k −3 fall-off of spectrum means the velocity fields are smooth; S2(r) ∼ r 2 for small r, shown in the inset of panel (c).(b) The polymer spectra Ep(k) ∼ k −3/2 follows from the scaling of E(k).(c) Plot of the second order structure function of second differences which scale as Σ2(r) ∼ r 3 .The exponent is same for the different De, although the range of scaling depends weakly on De.The inset shows the analytic scaling of S2(r).(d) Structure function of second differences, Σp, for various orders p for De = 9.(e) The exponents ζp versus p, calculated from the scaling behaviour of Σp.Departure from the straight line ζp = 3p/2 shows intermittency.Shaded region shows the standard deviation on the exponents computed from 18 snapshots.The two red + symbols mark the exponents ζ2 and ζ4 obtained with an alternate forcing scheme[39].That they lie well within error bars goes on to show that the results are independent of the large scale forcing.(f) Probability distribution function of δ 2 u(r) for four different values of r at De = 9.The distributions are non-Gaussian at small separations, while they become closer to a Gaussian (shown as a black dash-dotted curve) for large r.The corresponding cumulative distribution functions, computed using the rank-order method, are shown in the Supplementary Material, section C.

6 FIG. 4 .
FIG. 4. KurtosisThe kurtoses, (a) K2 and (b) K as a function of the scale r for De = 1, 3 and 9.The red dashed line is at ordinate equal to 3. We also show in (a) a line of slope −1.6.The scaling exponent of kurtosis, obtained from fitting the data in the gray shaded region, are: −1.6 ± 0.3, −1.6 ± 0.1, and −1.6 ± 0.1, for De = 1, 3, and 9, respectively.This demonstrates both the non-Gaussian nature of the PDFs and the universality of the exponents with respect to De.The kurtosis of δu, K, grows slower as r → 0 and may not be universal.To compare, we also plot, in (b), corresponding result for Newtonian HIT.Both the kurtoses K, K2 → 3 (shown in dotted-red line) as r → L. This indicates that at large separations the statistics of velocity difference are close to a Gaussian.

((FIG. 5 .
FIG. 5. Dissipation rates.(a) PDFs of the logarithm of the fluid energy dissipation rate ε f for all three De numbers.We denote by µ log ε f/p and σ log ε f/p the mean and variance of the logarithm of energy dissipation rates.The distributions deviate significantly from a log-normal behaviour in both the left and right tails.The right tails coincide for large De, similar to the coincident right tails at large Re in Newtonian HIT.(b) PDFs of the logarithm of the polymer energy dissipation rate εp for all three De numbers.The PDFs of log εp are sub-Gaussian, i.e. decay faster than a Gaussian, indicating εp is not intermittent.(c)The multifractal spectra of the fluid dissipation field calculated from the scaling of the energy dissipation rate εr calculated over a cube of side r.The black dash-dotted line shows the spectrum for Newtonian HIT[56].
To check if it does, we plot S 2 (r)/r 2 for three different De ranging from 1 to 9 in the Supplementary Material Fig.(S3b).At small enough r they all show S 2 ∼ r 2 .As r increases they all depart from this trivial scaling at a length scale ℓ which depends very weakly on De, if at all.This implies that even if a non-trivial scaling for S 2 exists in ET, it may require DNS at impossibly high De to be able to extract ζ 2 .Nevertheless, we have now demonstrated that in ET: (a) S 2 shows trivial scaling at small r, i.e., the velocity field is analytic, and (b) there is departure from the trivial scaling.Does the departure from the trivial scaling show a new scaling range?To explore this possibility we plot on a log-log scale S 2 (r) as a function of r in the Fig. (S3a) of the Supplementary Material, but it is unclear if there is a clear scaling range at intermediate r.Even if there is a non-trivial scaling exponent, it cannot be detected from the data, which is the highest resolution DNS of ET done so far.It often helps to detect a scaling range if we know beforehand what the scaling exponent is.In ET, unlike HIT, there is no theory that tells us what ζ 2 should be, but we do know that the Fourier spectrum of energy behaves like )This requires introduction of a new length scale ℓ which cannot depend on Reynolds number because in ET we are already in the range of small and fixed Reynolds number.The scale ℓ may depend on the Deborah number.