Anomalous scaling in the random-force-driven Burgers equation: A Monte Carlo study

We present a new approach to determine numerically the statistical behavior of small-scale structures in hydrodynamic turbulence. Starting from the functional integral representation of the random-force-driven Burgers equation we show that Monte Carlo simulations allow us to determine the anomalous scaling of high-order moments of velocity differences. Given the general applicability of Monte Carlo methods, this opens up the possibility to address also other systems relevant to turbulence within this framework.


1
The small-scale statistical properties of hydrodynamic turbulence is an old and tantalizing problem [1]. For turbulent flow stirred at large scales and far from the boundaries one expects a universal scaling for the small-scale fluctuations. Indeed, experiment gives strong indications for such universal behavior in Navier-Stokes turbulence [2][3][4][5][6]. The exact values of the scaling exponents however are still under debate. In such a situation it is useful to have a model system at hand that shares some essential properties with the original problem and allows for a clear physical understanding.
The random-force-driven Burgers equation is one such example. It was originally conceived as a one-dimensional model for compressible hydrodynamic turbulence [7] and provides a useful benchmark setting to test new analytical and numerical methods for real-world turbulence [8,9]. Here, u is the velocity, and f a centered random field displaying Gaussian statistics. We will consider the special case where the system is driven by a self-similar forcing that is white in time. The two-point correlation function of the stochastic forcing in Fourier space is given by where the parameter β determines the relative importance of the stirring mechanism at different scales, and the dimensionful constant D 0 measures its strength. For β large and negative the forcing effectively acts at large scales. On the other hand, the kinematic viscosity ν in (1) provides a dissipation scale η and for ν → 0 + the two characteristic scales η, and the system size L separate. The stochastic forcing drives the system into a non-equilibrium steady state, where in the range η ≪ |k| −1 ≪ L the energy flux through wavenumber k behaves as Π(k) ∝ |k| 1+β [10,11].
The case β = −1 corresponds to the physically interesting situation where the flux Π(k) is constant (up to logarithmic corrections) and the interplay of the stochastic forcing and advective term leads to a Kolmogorov energy spectrum E(k) ∝ |k| −5/3 [12,13]. The physical picture behind this scenario is the appearance of shocks with a finite dissipative width (see e.g. Fig. 1). The large fluctuations associated with the negative gradient of the front give the dominant contribution to the anomalous scaling of velocity differences ∆u = u(x+ r) −u(x).
In particular, we have the structure functions |∆u| n ∝ r ζn , and the scaling exponents ζ n = 1 for n ≥ 3 strongly deviate from the Kolmogorov scaling prediction ζ n = n/3 that follows from a naive dimensional analysis [8,9]. These rare fluctuations are strongly non-Gaussian and lead to the known asymptotic left tail of the probability distribution function (PDF) for velocity differences P(∆u, r) [14].
Here, we approach the problem from the functional integral point of view [15][16][17]. The functional integral gives a non-perturbative definition of the field theory and thus, it is ideally suited to study the strong and rare fluctuations present in fully developed turbulence that give the main contribution to the high-order moments of velocity differences. By sampling the associated probability distribution functional via Monte Carlo methods we show that it is possible to determine the scaling behavior of structure functions from first principles. Monte Carlo simulations are directly transferable to other systems of interest and are free of any modeling assumptions. Though not directly competitive with conventional time-advancing methods as, e.g. pseudo-spectral or finite-difference methods, Monte Carlo simulations may provide a unique perspective on such important problems as, e.g.
intermittency in fully developed turbulence [14]. In view of the well-established anomalous scaling behavior of Burgers turbulence [9,18] and the physical picture of the underlying mechanisms for intermittency [14,19,20], this provides an ideal setting to test our method and understand possible systematic effects at finite Reynolds number and system size. We emphasize again that in this paper we are not aiming to complete, or even improve on the accuracy obtained with other methods for Burgers turbulence. We rather want to provide a test of the generally applicable functional integral method and to demonstrate that a very reasonable accuracy can be obtained from this approach, a fact that was highly unclear at the beginning of this project.

The functional integral for the random-force-driven Burgers equation is obtained via the
Martin-Siggia-Rose formalism [15-17, 21, 22] by introducing an auxiliary response field µ.
We have the field theory with the action where D(x − y) is the spatial part of the two-point correlation function (2). In this form the action does not satisfy positivity. To obtain a Gibbs measure that can be sampled by a Markov chain Monte Carlo (MCMC) algorithm we integrate out the auxiliary field. This leaves us with the probability density functional The theory is then defined by placing the field u(x, t) on the sites of a regular space-time lattice Λ, i.e. (x, t) ∈ Λ. This way, we impose a UV cutoff that eliminates the details of those processes occurring deep in the dissipative regime. Then, the measure is given by [du] → (x,t) ∈ Λ du(x, t) and the action in (5) needs to be discretized appropriately. We replace the dynamics (1) with a finite-difference equation with backward-time discretization where ǫ is the lattice spacing in time direction. This ensures the correct dynamics in the continuum limit [23]. For the advective term we take the anti-symmetric spatial derivative where a is the lattice spacing in the spatial direction. With this choice of discretization the problem is amenable to a local over-relaxation algorithm [24]. Starting from an initial configuration {u(x, t), (x, t) ∈ Λ} the set of single-site variables is updated iteratively by 4 the successive application of a transition probability P (u(x, t) → u ′ (x, t)). We use the highquality ranlux (pseudo) random number generator [25] which is essential for large-scale lattice simulations. Specific improvements, e.g. Chebyshev acceleration [26] significantly reduce thermalization and autocorrelation times for the relevant observables in the inertial range.
It is necessary to map the discretized theory to its continuum counterpart and one has to ensure that the parameters are well-defined in the continuum limit. For that purpose the kinematic viscosity is identified with ν =ν a 2 /ǫ whereν is the viscosity in lattice units, and the Reynolds number scales as Re ∝ ν −1 . Furthermore, we have to ensure that the relevant scales of the system are resolved. In particular, we have to ensure that the dissipation scale fits on the lattice, i.e. η = Re −3/4 L a where L is the IR scale present in our system as a consequence of the finite lattice size. One may immediately recognize that this imposes a hard constraint on the realization of lattice simulations -fully developed turbulence requires a large computational effort where the number of lattice sites in the spatial direction scales as ∝ Re 3/4 , for given L. In practice, we are therefore bound to work at non-zero viscosity ν.
Simulations at moderate to high Reynolds numbers require massively parallel architec- The black curve indicates a bifractal scaling behavior.
appearance of multiscaling [18]. While in principle these contributions should be taken into account for the accurate determination of the scaling behavior, in practice it is difficult to distinguish different types of scaling contributions without any further assumptions. We obtain the scaling spectrum (Fig. 2b) where the error bars given are those of the LLS fit in the scaling range. Clearly, the n = 5 data point in Fig. 2b has minimal error which follows simply from our definition of the scaling range. We see that the scaling exponents are close to the bifractal scaling prediction [8,9], and within error bars agrees with the results of [18], obtained at high spectral resolution. As a last remark we want to add that we have not used extended self-similarity (ESS) [29] at any point in our analysis. Though ESS may enlarge the effective scaling range we found that it can suggest a clean scaling behavior even if subleading terms are present.
Since we are dealing with  In the continuum both the action in (5) and the measure are invariant under the set of Galilean transformations x → x + r , u(x, t) → u(x + r, t) + v , r = vt .
To avoid an over counting of physically equivalent field configurations one should eliminate these modes by the Faddeev-Popov procedure [23]. While gauge fixing is unavoidable for generic correlators [30,31] this is not so for velocity differences, as solely considered in this work which are manifestly invariant under transformations (8).
One may also check the statistics for velocity differences directly on the level of the probability distribution functions P(∆u, r). This gives valuable qualitative information on the physical behavior in our simulations of Burgers turbulence. In Fig. 3a we show the PDF of velocity differences for a set of values of the separation r, where we use the dimensionless variable φ = ∆u/[ ∆u 2 ] 1/2 to quantify the fluctuations. At large scales, far from the inertial range we clearly recognize the effects of the random Gaussian forcing (red). In the dissipative region the left tail of the PDF is especially pronounced and captures the strong fluctuations described by the shocks (orange). For separations η ≪ r ≪ L in the inertial range we see that the PDF P(∆u, r), plotted for three different values of r, nicely collapse onto each other (blue). In particular, in the regime where the fluctuations are much smaller than the root-mean-square velocity |∆u| ≪ u rms , the PDF of velocity differences has a universal scaling form where z is the dynamic exponent [32]. In the asymptotic region −∆u/r z ≫ 1 where ∆u < 0 we expect the algebraic scaling P(∆u, r) ∝ (∆u) γ with exponent γ = −4. The relevant region is shown in Fig. 3a (indicated by the arrow) and Fig. 3b. The corresponding scaling prediction with exponent γ = −4 is plotted for comparison as the black line in Fig. 3b.
Though our statistics are not sufficient to give a tight prediction on the scaling exponent, indications for the conjectured scaling behavior can be inferred from Fig. 3b.
At this point we want to give a short remark on some issues that arise when turning to incompressible three-dimensional Navier-Stokes turbulence. It is well-known, that the inclusion of the pressure term is one of the main obstacles in simulations of turbulence, as the requirement of incompressibility introduces non-local interactions. In the functional integral formulation this leads to a non-vanishing Fadeev-Popov determinant that can be treated by standard procedures (see e.g. [33]). In our lattice approach, we have chosen to use a local update over-relaxation algorithm that proved to be quite efficient for our purposes.
The long-range correlations in our system imposed by the forcing however, prohibit any attempt to parallelize in the spatial direction. This poses a severe problem when turning to higher dimensions. For that purpose it is absolutely mandatory to switch to a global update algorithm as, e.g. a Hybrid Monte Carlo algorithm [34] that is usually used in