Circulation in high Reynolds number isotropic turbulence is a bifractal

The turbulence problem at the level of scaling exponents is hard in part because of the multifractal scaling of small scales, which demands that each moment order be treated and understood independently. This conclusion derives from studies of velocity structure functions, energy dissipation, enstrophy density (that is, square of vorticity), etc. However, it is likely that there exist other physically pertinent quantities with uncomplicated structure in the inertial range, potentially resulting in huge simplifications in the turbulence theory. We show that velocity circulation around closed loops is such a quantity. By using a large databases of isotropic turbulence, generated from numerical simulations of the Navier-Stokes equations over a wide range of Reynolds numbers, we show that circulation exhibits a bifractal behavior at the highest Reynolds number considered: space filling for moments up to order $3$ and a mono-fractal with an unchanging dimension of about $2.5$ for higher orders; this change in character roughly at the third-order moment is reminiscent of a"phase transition". We explore the possibility that circulation becomes effectively space filling at much higher Reynolds numbers even though it may technically be regarded as a bifractal. We confirm that the circulation properties depend on only the area of the loop, not its shape; and, for a figure-$8$ loop, the relevant area is the scalar sum of the two segments of the loop.

From Leonardo da Vinci's half-a-millennium old drawings of turbulent motions in the river Arno [1] to their visualizations on modern day computers [2], evidence abounds that turbulent motion comprises organized structures, often evocatively described as "eddies" and "vortices". On the other hand, both phenomenological and analytical theories of turbulence [3,4] have largely focused on multi-point correlators of velocity, whose connection to the physical structures is not always clear. Furthermore, the most obvious multi-point correlators are best described as multifractals [5][6][7][8], which makes the problem very difficult to explore analytically. In this paper, we show that the velocity circulation around closed loops, besides providing a plausible link between vortical structures and statistical objects, has a very simple bifractal structure: space filling for moments roughly up to order 3 and a mono-fractal with an unchanging dimension of about 2.5 for higher orders. We comment on possible Reynolds number effects.
For reference, the circulation around a loop of linear dimension r is defined as where ∂D r denotes the boundary of the loop, u is the velocity, ω ≡ ∇ × u is the vorticity, dl is an elemental length along ∂D r andndA is an elemental area of D r . Migdal [9], who initiated the statistical theory of circulation, made the case that, in the inertial range that is far * krs3@nyu.edu from both the forcing and dissipation scales, the probability density function (PDF) of Γ r depends uniquely on the scaling variable G = Γ 2k r /r 4k−2 , if the loop size r lies within the inertial range, where k is some parameter. Kolmogorov's arguments (henceforth K41) [3] lead to k = 3/2, with the circulation moments scaling as Γ p r ∼ r 4p/3 . Migdal also argued that the PDF depends only on the area circumscribed by a simply-connected loop, and not on its actual shape, as long as it is entirely contained in the inertial range. This is the area rule. Another result (call it the figure-8 area rule) is that, when the loop is a figure-8 with a common vertex, the circulation statistics depend only on the scalar sum of the areas of the component segments of figure-8, rather than their vectorial sum as one traverses the loop.
This theory was followed up soon after by a small number of experimental and numerical papers [10][11][12][13][14]. Because they were all limited to low Reynolds numbers, the verification of the area rules was stymied by the modest extent of the inertial range; and, obviously, they could not evaluate the high-Reynolds-number properties. Further, the experimental flows were not homogeneous, which thus complicated the inferences drawn therefrom. It is thus not a great surprise that these earlier studies did not agree quantitatively among themselves. One common inference of these early studies is that circulation is highly intermittent, just as the velocity increments are, and display multifractal scaling with no unique value of k, making no new insights possible. We assess these properties persuasively here by taking recourse to direct numerical simulation (DNS) data of statistically stationary, homogeneous and isotropic turbulence in a periodic box over a wide range of Reynolds numbers [15], including the highest of any so far [16]. We also show that the scaling variable G = Γ 3 r /r 4 , corresponding to k = 3/2 in Migdal's expression [9], stemming from K41, shows a very good collapse in the inertial range.
Migdal did not have conclusive thoughts on the scaling of circulation moments. We show that circulation at the highest Reynolds numbers of this study is a bifractal, to which we have already made an allusion; this simple behavior contrasts a highly intermittent, multifractal structure typical of low Reynolds numbers. This result goes a considerable distance in addressing the question: What statistical variable in turbulence would be most apt to study, among the infinite number of possibilities?
For later purposes, we provide the following definitions. The longitudinal velocity increment is defined as ∆ r u ≡ u(x + r) − u(x), where the velocity component u and the separation distance r are both taken in the same direction. The inertial range is defined as that range where the normalized third-order velocity structure function (∆ r u) 3 /r ǫ , ǫ being the global mean value of the energy dissipation rate, is equal to the exact theoretical value of −4/5 [17]. This result has been evaluated in Ref. [18], for the same data as the present.

I. DATA
The DNS data used in this work have been acquired by solving the incompressible Navier-Stokes equations, where u is the solenoidal velocity field (∇ · u = 0), p is pressure, ρ is fluid density, ν is the kinematic viscosity and f is the forcing term that maintains a stationary state [19,20]. We use Fourier pseudo-spectral calculations [21] on a periodic domain of size (2π) 3 with an explicit second order Runge-Kutta integration in time. A combination of phase-shifting and truncation is used to reduce aliasing errors, where the highest resolved wavenumber k max = √ 2N/3 and N is the number of grid points in one direction. Typical spatial resolution, expressed by k max η, was fixed around 1.5 in earlier simulations [22]. Recently [16], it has been pointed out that the spatial and temporal resolution are more stringent at higher Reynolds numbers. For some of the data analyzed here, data were obtained with improved resolution. Table I lists some flow parameters of interest.
One other comment is useful. Circulation Γ r around a loop of side r is calculated using the second equality in Eq. 1, which follows from the Stokes theorem, as the two-dimensional local average of vorticity using the algorithm given in Ref. [23]. Statistical averages were taken over the whole N 3 simulation along the three Cartesian directions. Statistics of Γ r were also calculated using the loop integration of Eq. 1, by means of cubic splines for improved accuracy, and excellent confirmation of the area integral results were obtained for r/η > 5. I. Isotropic DNS database: N 3 is the number of points on a L 3 0 grid with L0 = 2π units, R λ ≡ u ′ λ/ν is the Taylor scale Reynolds number where u ′ is the root-meansquare velocity fluctuation, λ ≡ u ′ / (∂u/∂x) 2 is the Taylor microscale, ν is the kinematic viscosity, L ≈ L0/5 is the integral scale, η ≡ (ν 3 / ǫ ) 1/4 is the Kolmogorov scale, ∆x = L0/N is the grid spacing. Results have been averaged over a time span of at least 10 large-eddy time scales (L/u ′ ), except for the 16, 384 3 data for which the averaging time is much shorter.

II. RESULTS
A. Area rules Figure 1 shows circulation traces for inertial range separations, normalized by its standard deviation, along the length of the simulation domain of size L 0 , at the Taylor scale Reynolds number, R λ = 1300. Two different inertial separations r corresponding to the lower end and the middle of the inertial range plateau in − (∆ r u) 3 /r ǫ , are shown here. These particular signals do not show frequent excursions to very high values (unlike highly intermittent phenomena). Figure 2 shows the PDFs of Γ r for a number of rectangular loops of the same area but differing aspect ratios. The data collapse to a very good accuracy. Even though the result will have to be confirmed for loops of different shapes, this figure supports the expectation that only the area of the loop, not its shape, decides the PDF of circu-  The data follow the scalar area law, as evidenced by the A 2/3 scaling in the inertial range. The linear area dependence in the small-A regime is shown, and the tendency to saturation for large areas is also apparent.
lation. It is thus sometimes convenient to use the symbol Γ A to describe the circulation around a loop of area A; henceforth, when we speak of Γ r , it means circulation around a square of side r.
The inset to 3 shows a figure-8 loop with two different squares touching at a common vertex. If K41 is valid, it readily follows that the standard deviation of circulation On the other hand, if one traverses along the loop in the direction of the arrows marked on the loop, the areas circumscribed by the two squares will have different signs. If the resulting vector area is the right quantity to use, the standard deviation Γ 2 A 1/2 for the figure-8 loop should scale as   3 shows convincingly that the mean-square circulation varies as A 2/3 , where A is the scalar sum of the areas of the two loops, for most of the range. For small A, the variance is clearly linear in A as expected from Taylor expansion, whereas it saturates for large A because of many cancellations.

B. Application of Kolmogorov's similarity argument
As already stated, straightforward application of K41 shows that Γ r ∼ r 4/3 . Migdal [9] argued that the PDF of Γ r decays as some power of Γ r 3 /r 4 . We show the PDF of Γ 3 r /r 4 ǫ in Fig. 4 for various inertial range separations. The PDFs show good collapse, with some deviations in the negative PDF tails; the two insets, highlighting the positive and negative tails of the PDF, are not power laws. They can be fitted nominally by unequal stretched exponentials, as described in the caption to Fig. 4. The mean and the mean square of this distribution yield the third-order and sixth-order moments of circulation. These and other moments are computed separately in Sec. III A.
The results so far are already interesting and tantalizing because they contain the suggestion that K41 may have currency for circulation. This would be quite unlike, say, for velocity increments or enstrophy density. We now explore this feature further by computing the scaling properties.

A. Exponents
We first evaluate various even order moments of Γ r and show moments Γ p r for p = 2, 4, 6 and 8 for R λ = 1300 as functions of r in Fig. 5. They all display proper power laws, Γ p r ∼ r λp , in the inertial range, as reinforced by the constancy of the local slopes, shown in the inset of Fig. 5.
The odd moments of Γ r do not display equally clean power laws because they have negligible intensity via cancellation, which leads to poor convergence. Their scaling improves if one considers absolute values of circulation, |Γ r | p ∼ r λ |p| (while, obviously, those of even orders remain unchanged). Absolute values enable us to define scaling exponents for fractional p as well, up to (but not including) p = −1, for which the moments diverge [24]. The practice of using absolute moments for odd orders is justified, at least a posteriori, as long as the exponents so obtained are monotonic when plotted together with evenorder data. This monotonic behavior is not quite true for velocity structure functions even for very high Reynolds numbers typical of atmospheric flows (see Ref. [25] at R λ = 10, 340) but seems to hold for circulation for which the absolute moment data do not zigzag with respect to even-order data.
In Fig. 6 we plot the scaling exponents λ |p| for circulation as a function of the power index p for all orders at R λ = 1300. The exponents seem to organize themselves into two straight lines, one below about 3 and the other above it. The low-order data can be fitted asymptotically by the K41 line, λ |p| = 4p/3, but one needs a different line with a smaller slope for higher orders. The latter line (with the slope of 7/6 and an intercept of 1/2) can  6. Scaling exponents as a function of moment order p; R λ = 1300. The data can be fitted by two separate lines, one below order 3 and one above it. The low order data can be fitted by the expression 4p/3 (dashed line), which results from the application of K41. The high-order data can be fitted (solid line) by a fractal model (Eq. 4) with D ≈ 2.5. Least square fits are shown. Typical error bar shown for the tenth moment is subsumed by the symbol thickness. The difference between the K41 line and the measured exponents is indicated by δλ |p| , which will be discussed in Sec. IV.
be expressed in terms of the β-model [26], with D ≈ 2.5 as the self-similarity dimension [26,27]. We have so far considered results for the two highest Reynolds numbers of the dataset computed by us. It is instructive to examine how the results depend on the Reynolds number, and assess the asymptotic state.

B. Probability density functions
We have already seen the PDFs in two forms: The PDF for circulation around various shapes of loops within the inertial range in Fig. 2 and those for Γ 3 r /r 4 ǫ in Fig. 4. We now plot the PDFs of Γ r for several values of r in the inertial range, normalized by their own standard deviations. The PDFs collapse on each other for magnitudes below about 5 standard deviations (see inset (a)) and differ mostly in the tails. This shows that we do not have perfect self-similarity. The PDFs also depart strongly from the Gaussian distribution indicated by the dashed line, both at the tails and the core region (see inset (a)). The peaks of the PDFs gradually drop towards the Gaussian value of about 0.4 through the inertial range, as can be seen more explicitly in inset (b), but the tails are of the stretched exponential type ∼ exp(−βΓ α(r) r ). The stretching exponent α(r) is plotted in inset (c) of Fig. 7 as a function of the inertial range separation r; we have α(r) varying weakly with r, supporting the notion that circulation has weak intermittency [28]. We note that the variation of α(r) in the inertial range is gentler than   that of the corresponding stretching exponents of velocity increments, which are known to rapidly change with scale [29]. FIG. 9. The peak magnitude of logarithmic local slope of circulation flatness in the inertial range, plotted against R λ . The solid line is a least-square fit. If extrapolated, the minimum of the logarithmic derivative will reach zero, corresponding to constant flatness in Fig. 8, at R λ ≈ 1900. Data of Ref. [12] at R λ = 216 ( ), is shown for comparison. Inset shows the logarithmic local slopes of circulation flatness as a function of r/η, up to r ∼ O(L) for R λ = 1300 (▽) and R λ = 240 (•).

C. Flatness
The circulation flatness is defined as Figure 8 shows that F (r) varies through the inertial range (η ≪ r ≪ L) but its evolution with Reynolds number is the more interesting point; the figure also compares F (r) at three different Reynolds numbers. At R λ = 240, F (r) increases smoothly from the Gaussian value of 3, for r ∼ O(L), to the dissipation range limit of ω 4 / ω 2 2 , where ω is the vorticity component normal to the plane of circulation. With increasing R λ , however, F (r) appears to develop a plateau in the inertial range, indicating that it may be approaching a constant, independent of scale. The inset of Fig. 8 compares the circulation flatness to those of the longitudinal and transverse velocity increments at R λ = 1300. Here, the transverse increment where the separation distance r is transverse to the velocity component v. Even at R λ = 1300, the flatness factors of both ∆ r u and ∆ r v smoothly grow with decreasing scale, showing that the velocity increments are highly intermittent, whereas F (r) displays the tendency towards constancy, suggesting that Γ r at high Reynolds numbers is only weakly intermittent. It should be stressed that, at lower R λ , all flatness factors of Γ r , ∆ r u and ∆ r v increase rapidly with decreasing scale in the inertial range, as already shown in Refs. [12,14]. Our point is that the flatness F (r) has the potential to become a constant in the inertial range of r as the Reynolds number increases further. If so, this will be the intermittency-free limit in which we expect the logarithmic local slope of flatness d[log F (r)]/d[log r] → 0 in the inertial range. To quantify this approach to the intermittency free limit, we plot the logarithmic local slope of F (r) as a function of spatial separation in the inset of Fig. 9 for two different Reynolds numbers. For R λ = 1300, the local slopes are closer to zero over a wider range of inertial scales than for R λ = 240. The peaks of the flatness local slopes decrease linearly with R λ as shown in Fig. 9, suggesting that it will be zero at R λ ≈ 1900. This result should be regarded merely as an indication and will be discussed further in Sec. IV A below.

IV. DISCUSSION
Small-scale turbulence is known to be characterized by extreme events in time and space. Such intense events can be seen in local velocity increments, which exhibit intense spatial fluctuations, resulting in scaling exponents that vary nonlinearly with respect to the moment order. This is the phenomenon of intermittency that necessitates the superposition of infinitely many scale-invariant configurations to describe the exponents. Much effort has been expended in quantifying the nonlinear trend of the intermittency exponents of velocity moments, with varying degrees of success [30][31][32][33]. As noted in Ref. [34], intermittency corrections from the velocity moments serve as an upper estimate, since the velocity field is infrared divergent. Here, using the largest simulations of isotropic turbulence to-date, we have shown that the structure of velocity circulation is much simpler at higher Reynolds numbers. In contrast to those of velocity increments, the circulation statistics become less intermittent with increasing Reynolds numbers, with the exponents having an approximately bifractal structure: at the level of low order moments, circulation is essentially a space-filling quantity whereas, for moments of order 3 and higher, it appears to have a self-similar dimension D ≈ 2.5.
These findings brighten the prospect of a simplified turbulence theory that can be used to describe smallscale turbulence without having to resort to multifractal models. We have inferred that circulation, though arising from highly convoluted vortex structures in space and time, is essentially space filling for low-order moments and can give rise to a bifractal scaling. It will be interesting to obtain circulation statistics in more realistic cases of anisotropic turbulence, to see if a statistical theory based on vortex filaments can shed new light on topics such as anisotropy effects on small-scale universality.
Two questions appear worth discussing in some detail. (a) Since the data reveal that there is a Reynolds number dependence of the circulation properties, it is worth asking whether they have reached its asymptotic state even at the highest Reynolds number considered here, and, if not, make an educated guess on that state. This point is essentially an expansion of the tentative result deduced from Fig. 9. (b) Since circulation is very closely related to vorticity, it is natural to ask why the simplicity that is apparent in circulation does not translate to enstrophy density, which is known to be a strongly multifractal quantity [35]. We will discuss these questions in that same order.
A. The asymptotic state Figure 10 shows the behavior of δ λp , the difference between the measured scaling exponents and the corresponding Kolmogorov values, for moment orders 4, 6, 8 and 10, as functions of the Reynolds number. For each order, this difference seems to approach K41 roughly as a power law. The rates of approach vary inversely with the order of the moment. We cannot speculate about the behaviors at infinitely large Reynolds numbers, but may expect that this behavior will persist up to some higher Reynolds number. The fact that the approach to the K41 values is slower for the high-order moments suggests that, in principle, the scaling will remain a bifractal for all finite Reynolds numbers, with the "phase transition" point moving to a higher p with increasing Reynolds numbers. In practice, however, moments may fall on the K41 line for sufficiently high values of p that the bifractal behavior may essentially yield place to a space-filling monofractal.

B. Circulation, enstrophy and velocity increments
Since circulation is very closely related to vorticity, it is natural to ask why the simplicity that is apparent in circulation does not translate to locally averaged enstrophy density, which is known to be a strongly multifractal quantity [16,35]. The corresponding question is relevant also with respect to velocity increments. We address this issue here briefly.
FIG. 11. The function ψ(r) (left hand side of inequality 17) as a function of normalized spatial separation r/η for different Reynolds numbers. The curves are seen to collapse, for r ≪ L (L is the integral scale) at the different Reynolds numbers shown, reasonably well. The limit ψ(r → 0) = 1/3 is marked by the horizontal dashed line. Solid black line indicates the corresponding K41 slope.
Under homogeneity, the first moments of dissipation, enstrophy and their local averages are related as ν Ω r = ǫ r = ǫ , using which, inequality (14) (corresponding to p = 1) can be written as Using p = 1 and the inertial range result (∆ r u) 3 = (−4/5)r ǫ , we end up with the same inequality as Eq. 17. Inequality 17 is more general, since it shows that ψ(r) is bounded for all r. As a check, we note that in the dissipative limit r → 0, we have lim r→0 ψ(r) = lim r→0 Γ 2 r ν r 4 ǫ = ω 2 ν ǫ = 1 3 where the relation Ω = 3 ω 2 uses isotropy. It is interesting to consider the double limit of ψ(r), namely Since the mean dissipation is known not to depend on viscosity [2,[36][37][38], the term within square brackets that is r-independent, tends to 0 in the asymptotic limit. This means that in the small-r limit we have, Γ 2 r /r 4 = ω 2 → ∞ for high Reynolds numbers. It follows that vorticity can attain very high moment values in turbulent flows even when those of circulation are bounded. This observation does not necessarily imply that circulation (around finite contours) is always bounded, but it does suggest that it is a more tractable quantity to work with, as compared to vorticity. Figure 11 shows that the DNS data are quite consistent with the dissipative limit. Furthermore, ψ(r) at different Reynolds numbers shows good collapse for r ≪ L, where L is the integral scale.
These arguments show that the simplicity in the structure of circulation does not contradict the known multifractal properties of velocity increments or enstrophy density. We thus believe that we have covered new ground here.
V. ACKNOWLEDGMENTS