Dynamic multiscaling in stochastically forced Burgers turbulence

We carry out a detailed study of dynamic multiscaling in the turbulent nonequilibrium, but statistically steady, state of the stochastically forced one-dimensional Burgers equation. We introduce the concept of interval collapse time, which we define as the time taken for a spatial interval, demarcated by a pair of Lagrangian tracers, to collapse at a shock. By calculating the dynamic scaling exponents of the moments of various orders of these interval collapse times, we show that (a) there is not one but an infinity of characteristic time scales and (b) the probability distribution function of the interval collapse times is non-Gaussian and has a power-law tail. Our study is based on (a) a theoretical framework that allows us to obtain dynamic-multiscaling exponents analytically, (b) extensive direct numerical simulations, and (c) a careful comparison of the results of (a) and (b). We discuss possible generalizations of our work to higher dimensions, for the stochastically forced Burgers equation, and to other compressible flows that exhibit turbulence with shocks.


I. INTRODUCTION
Studies of statistically homogeneous and isotropic turbulence, one of the most important examples of a non-equilibrium steady state (NESS), lie at the interfaces between non-equilibrium statistical mechanics, fluid dynamics, and spatiotemporal chaos.In this NESS, the turbulent fluctuations of the flow velocity span wide ranges of spatial and temporal scales; and the correlation or structure functions, which characterize these fluctuations, display power laws in space and time [1][2][3][4][5][6][7][8].It has been suggested several times that these power-law forms are like those seen in correlation functions in equilibrium critical phenomena [9][10][11][12][13].Perhaps the simplest example of a power-law form in such turbulence is the scaling of the energy spectrum, E(k) ∼ k −α , for wave numbers k in the inertial range η −1 k L −1 I , where η and L I are, respectively, the dissipation and energyinjection length scales [1,14].The phenomenological theory, proposed by Kolmogogorov in 1941 (K41), yields the exponent α K41 = 5/3, but the intermittency of turbulence leads to multifractal corrections to this and to other exponents (see below and, e.g., Refs.[1,[14][15][16]).To understand this multifractality we must go well beyond [1,[14][15][16] the theoretical framework that is used to explain power-law correlation functions at equilibrium critical points [9][10][11].
We elucidate dynamic multiscaling in the NESS of the stochastically forced Burgers equation.We show that the theoretical methods that have been used to study such multiscaling of turbulence in incompressible hydrodynamical PDEs cannot be used for turbulence in the stochastically forced Burgers equation (often referred to as Burgers turbulence) [21][22][23][24][25], which is compressible.
Therefore, we introduce an infinity of time scales that can be used to characterize the spatiotemporal evolution of Burgers turbulence in d spatial dimensions.Our study is based on (a) a theoretical framework, for d = 1, that allows us to obtain dynamic multiscaling exponents analytically and (b) extensive direct numerical simulations (DNSs) for d = 1 and d = 2.To the best of our knowledge, this is the first time that such multiscaling exponents have been obtained analytically for turbulence in a nonlinear hydrodynamical PDE.Before we present the details of our study, we give a qualitative overview of our work and its significance in the light of earlier studies of related problems.
We can extract time scales from turbulent flows in a variety of different ways.The K41 theory [26] suggests simple dynamic scaling, i.e., all time scales are characterized by the same dynamic exponent θ K41 = 2/3.However, from DNSs and a heuristic understanding of various models of turbulence, we have learned that incompressible turbulence exhibits dynamic multiscaling, i.e., different time scales are linked to length scales via different dynamic exponents [4-8, 20, 27].
L'vov, Podivilov, and Procaccia [4] had proposed that the characteristic times of eddy-velocity correlation functions of various orders are linked to the eddy sizes through multiple dynamic exponents, which can be related to the equal-time structure function exponents by linear bridge relations.Thereafter, Mitra and Pandit [5] not only confirmed this from their shell-model DNS, but also showed that the dynamic exponents depend specifically on how the time scales are defined.In addition, for the passive-scalar problem, Mitra and Pandit [28] showed that simple dynamic scaling is obtained if the velocity field is of the type in the Kraichnan model [29,30]; but bona fide dynamic multiscaling, with an infinity of nontrivial dynamic exponents, is obtained only if the advecting velocity itself exhibits dynamic multiscaling.We note that the analysis of such passive-scalar turbulence is much simpler than its fluid-turbulence counterpart because the passive-scalar equation is linear in the scalar concentration.For fluid turbulence, we must confront the nonlinearity of the incompressible Navier-Stokes equations.
Typical shell models, which are coupled ordinary differential equations (ODEs) designed to mimic turbulence in hydrodynamical PDEs, are defined in a logarithmically discretised wave-number (Fourier) space, with shells, labelled by the wave number k, and the complex, scalar shell velocity u k .Given their simplicity, such shell models cannot yield reliable flow fields; however, they have been remarkably successful in obtaining turbulent states with spectral and multifractal properties akin to those that are obtained for turbulence in their parent hydrodynamical PDEs (see, e.g., Refs.[14,31]).The interaction between velocities in shell models is in marked contrast to the interaction of Fourier modes of the velocities in hydrodynamical PDEs.
For instance, in the incompressible Navier-Stokes equation, the nonlinear term couples all velocity Fourier modes to each other, whereas shell-model velocities u k interact only with velocities in nearest-or next-neighbor shells, so small-k velocities do not affect large-k velocities directly.By contrast, in hydrodynamical PDEs, small eddies (largek modes) are advected directly by large ones (small-k modes), because of the coupling of all Fourier modes to each other, whence we get a sweeping effect [32] that yields eddy lifetimes that are linearly proportional to the eddy size, i.e., the dynamic exponent is unity.This sweeping effect, which masks the underlying dynamic multiscaling of turbulent flows, is a manifestation of Taylor's hypothesis.
For incompressible turbulence, Belinicher and L'vov [33] had suggested that quasi-Lagrangian (QL) velocities, calculated in the reference frame of a Lagrangian particle or tracer, should be free from sweeping effects.This was shown explicitly by Ray, Mitra, Perlekar, and Pandit [7], who quantified the dynamic multiscaling of forward-cascade turbulence in the incompressible two-dimensional (d = 2) Navier-Stokes (NS) equation.
Biferale, Calzavarini, and Toschi [27] used QL velocities to obtain dynamic multiscaling for turbulence in the incompressible threedimensional (d = 3) Navier-Stokes equation.Ray, Mitra, Perlekar, and Pandit [7] also showed that sweeping effects can be suppressed by friction, which removes energy from a flow at all spatial scales.
The characterization of dynamic multiscaling in incompressible-fluid turbulence via the QL approach cannot be used when we consider turbulence in a compressible fluid.The Lagrangian particles, which we require to define QL fields, get trapped in shocks, which form in compressible turbulent flows.We show this explicitly by DNSs for turbulence in the onedimensional (d = 1) stochastically forced Burgers equation.Before we show how to overcome this difficulty, it is useful to recall some elementary results for Burgers turbulence.
The Burgers equation, the simplest compressible hydrodynamical PDE, which was introduced for the study of fluid equations in the low-viscosity limit and then for pressure-less gas dynamics [34,35], has the same nonlinear term as the NS equation; so it is often used as a testing ground for statistical theories of turbulence [21,22].The d-dimensional Burgers equation can be solved via the Hopf-Cole transformation [35][36][37]; in the inviscid limit, this yields a maximum principle for a velocity potential [21,22,38].In cosmology, the Burgers equation is used to model the formation and distribution of large-scale structures in the universe [39], under the adhesion approximation [40].
Studies of Burgers turbulence (sometimes referred to as Burgulence) comprise investigations of the statistical properties of (a) solutions of the Burgers equation with random initial data or (b) solutions of the Burgers equation with stochastic forcing [25,41,42]; we concentrate on the latter.
The stochastically forced Burgers equation can be derived by taking a spatial derivative of the Kardar-Parisi-Zhang (KPZ) equation [43,44], which models the height profile of a growing interface.For the physically relevant forms of the stochastic forcing used in the KPZ equation, correlation functions display simple scaling of equaltime and time-dependent correlation functions.
We consider turbulence in the d = 1 Burgers equation with a zero-mean, white-in-time, Gaussian random force whose variance, in Fourier space, scales as ∼ k −β .This stochastically forced equation yields a NESS with properties that are akin to those in d = 3 NS turbulence.In particular, one-loop renormalization-group (RG) calculations yield a K41-type energy spectrum E(k) ∼ k −5/3 for β = 1, for k in the inertial range [2, 3, 23-25, 41, 42, 45].DNSs show that velocity structure functions (see below) in this NESS display multiscaling [2, 3, 23-25, 41, 42, 45].It has been conjectured that this multiscaling could be a numerical artifact [25] which could wane with increasing spatial resolution of the DNS, and be replaced by bifractal scaling.
Given the challenges that we have outlined above, the study of dynamic scaling in the stochastically forced Burgers equation is not as well developed as it is in incompressible fluid turbulence [4-8, 20, 27].An earlier DNS-based study of certain time-dependent, Eulerian velocity structure functions for this model yielded a single dynamic exponent of unity [2] that was attributed to the sweeping effect.We go beyond this Eulerian study and the QL investigations of incompressible fluid turbulence [4-8, 20, 27].
As we have noted above, Lagrangian tracers get trapped at shocks in compressible turbulent flows, because of which the QL transformation might not be adequate for the removal of sweeping effects.For turbulence in the d = 1 stochastically forced Burgers equation, we overcome this difficulty by using a pair of tracers separated initially by a Lagrangian interval of length .We then compute the interval-collapse time, τ col , which we define as the time at which this pair collapses to a point at a shock.We find that τ col depends on both and the location of the interval.Hence we compute, for each value of , the probability distribution function (PDF) of τ col and extract a hierarchy of time scales, T p col ( ), from its orderp moments.We make the dynamic-scaling Ansätze, T p col ( ) ∼ z p col , and obtain therefrom the intervalcollapse exponents z p col for different values of p.We find from our high-resolution DNS that z p col is not a linear function of p.This indicates dynamic multiscaling and intermittency.We develop a theory for the pdependence of z p col , which indicates that z p col should be a piecewise linear function of p, i.e., we obtain a particular case of dynamic multiscaling, namely, dynamic biscaling.Our theory also yields an analytical expression for the PDF of τ col .The PDF is non-Gaussian and has a power-law tail; we compare this with the results of our DNS.We also show from our DNS that the QL approach is indeed unable to suppress sweeping effects.Furthermore, we examine the extension of our interval-collapse framework for examining dynamic scaling in the stochastically forced Burgers equation to dimensions d > 1 and the challenges in carrying out such a study.
The remainder of this paper is organized as follows: In Sec.II we define the stochastically forced Burgers equation and outline the numerical methods that we use.In Sec.III, we investigate the statistical properties of the interval-collapse times, τ col , and show how to calculate the interval-collapse exponents z p col .In Sec.IV, we explore dynamic multiscaling via QL velocity structure functions.Finally, in Sec.V, we discuss the significance of our results and propose ways of extending our approach to d > 1 and compressible turbulent flows with shocks.

II. MODEL AND NUMERICAL METHODS
The d = 1 stochastically forced Burgers equation that we consider is where u(x, t) is the fluid velocity at position x and time t, ν is the kinematic viscosity, and f is a zero-mean, Gaussian white-in-time random force whose Fourier components f (k, t) satisfy where k is the the wave number.We choose β = 1, because, at the level of a one-loop RG, this choice of β yields a K41-type energy spectrum , for k in the inertial range, where û(k, t) is the Fourier transform of u(x, t), and the angular brackets denote a time average over the NESS [2,24,25,41].
Here, η = (ν 3 / ) 1/4 is the dissipation scale beyond which viscous losses are significant, is the mean energy dissipation rate, and is the integral length scale.In our DNSs of Eqs. ( 1) and (2), we use periodic boundary conditions in a domain of length L = 2π and N collocation points.To achieve high spatial resolutions we use N = 2 16 (Run R1) and N = 2 20 (Run R2).We employ a standard pseudospectral method with the 2/3-dealiasing rule [46].For time-stepping we use the implicit Euler-Maruyama scheme [47,48].We generate the stochastic force f (k, t) in Fourier space, with a high-wave-number cutoff at k c = N/8.We have confirmed that our results remain unchanged even if we FIG. 1. (Color online) Space-time plots of representative tracer trajectories (run R1) for initially equi-spaced tracers; we choose 1000 tracers for ease of visualization here.The vertical axis gives the nearest-neighbor tracer-pair number (the tracers are numbered in ascending order and the pair number i denotes two tracers that are initially at 2πi/1000 and 2π(i + 1)/1000); the horizontal axis denotes the time (in units of iteration number).The color bar denotes the separation between the tracers within each of the nearest-neighbor tracer pairs.As time progresses, the separation between the tracers in such pairs decreases.Bottom panel: Expanded versions of the tracks of the tracers (distinguished by different colors) enclosed within the redbordered rectangular boxes in the top panel.
evolve the unforced equation by using the second-order exponential time-differencing Runge-Kutta method and then add the forcing term, f (k, t) √ δt, to the velocity field at the end of every time step of step size δt.The parameters for our DNS runs are given in Appendix A. Most of our DNSs and data analysis have been carried out on a GPU cluster with NVIDIA Tesla K20 accelerators.Once the system reaches its NESS, the energy spectrum E(k) shows inertial-range scaling over two decades in k; a typical snapshot of a steadystate velocity profile and a compensated plot of E(k), for run R2, are shown in the Appendix A.
After our system reaches its NESS, we introduce N p equi-spaced Lagrangian particles (or tracers), whose equations of motion are where X i and U i are the instantaneous position and velocity, respectively, of the i−th tracer, and δ(x) is the Dirac delta function.We solve (3) by using the forward-Euler method.The δ(x − X i ) factor is implemented by linear interpolation.As time progresses, the tracers cluster at the shocks, as we show in Fig. 1.

FIG. 2.
Examples of the possible cases of collapse of an interval of initial length r from our DNSs: Panel (A) The collapsing interval contains a shock at t = 0 and it collapses at that shock at t = τ .Panel (B) The interval has no shock at t = 0 and collapses, at t = τ , at a shock which either appears within it (upper row) or merges with one of its ends (lower row) at time t = t * .In all figures, t, t * and τ have been non-dimensionalized with TL.

III. INTERVAL-COLLAPSE TIMES
We define τ col ( ) as the time taken for an interval, with initial length , to collapse to a point at a shock, with t = 0 the time at which we seed the flow with tracers.We consider many such intervals at t = 0 and address the following questions: • What is the PDF of τ col ( )?
• How does this PDF depend on ?
We start by investigating the moments of this PDF.

A. Dynamic Scaling: Theoretical Considerations
We extract a hierarchy of time scales, T p col , by normalizing order-p moments of the PDF of τ col by the the large-eddy turnover time, T L = L I /u rms , with u rms the root-mean-square velocity, as follows: • x denotes an average over different spatial locations x and r ≡ /L is the non-dimensionalized interval length.The velocity difference, δu(r), across inertial range separations η/L r 1 in a turbulent flow scales as δu(r) ∼ v 0 r h , where v 0 ≡ u 0 /u rms , the largescale fluctuating velocity field is u 0 , and h is the scaling exponent of the velocity fluctuations.Theoretical considerations suggest [25] that the NESS of Eqs. ( 1) and (2) exhibits the following bifractal scaling, for equal-time velocity structure functions: This bifractal scaling is obtained as follows: For an interval of length r, (a) δu(r, x) ∼ v 0 r h , if the interval does not contain any shock; (b) δu(r, x) ∼ v 0 , independent of r, if it does contain a shock.The value of the exponent h depends on β; for β = 1, h = 1/3.In the limit r 1, an interval of length r can have at most one shock with a probability p r ∝ r [49].If we substitute these expressions for δu(r, x) in S p (r), we obtain the result for ζ p in Eq. ( 5).The DNS results of Ref. [25] yield multiscaling; but it has been suggested in this study that such multiscaling might be an artifact that should be replaced by the bifractal scaling [Eq.( 5)] in the limit of infinite spatial resolution.The exponent z p col versus p from runs R1 (pink triangles) and R2 (blue triangles).The solid line shows our theoretical prediction, Eq. ( 18).Inset: The local slope m(r) as a function of r, the color and symbol coding being the same as that in (a).The means and standard deviations of m(r), calculated over the shaded region, yield the exponents z p col and their error bars, respectively.
We generalize these theoretical arguments to obtain exact expressions for the dynamic scaling exponents of the NESS of Eqs. ( 1) and (2).At any instant of time t, consider two Lagrangian particles at the two ends, x and x + y, of an interval of length y.Then where ẏ = dy/dt.The time of collapse of an interval τ (r, x) ≡ τ col (r, x)/T L is obtained by integrating Eq. ( 6): In the small-r limit, we have the two cases A and B, which we describe below and illustrate in Fig. 2.
(A) There is a shock within the collapsing interval of size r at t = 0 (see panel (A) in Fig. 2).In this case δu(y, x) ∼ v 0 is independent of y, for all y ∈ [0, r].Hence, from Eq. ( 7) we obtain The probability that an interval of size r contains a shock is proportional to r (this probability is ar, with a a constant that does not affect the power laws that we obtain below), so (B) There is no shock in the interval at t = 0 (see panel (B) in Fig. 2).As we have noted below Eq. ( 5), δu(r, x) ∼ r h for this case.However, either a shock forms within the collapsing interval (first row in panel (B) of Fig. 2) or one of the ends of the interval gets trapped at a shock (second row in panel (B) of Fig. 2), at some time t * .The probability p s1 of finding an interval of initial size r, which collapses as in the first row of panel B in Fig. 2, is p s1 ≡ br, where b is a constant that does not affect our results for exponents.Similarly, the probability of finding an interval, which collapses as in the second row of panel B in Fig. 2, is p s2 ≡ 1 − (a + b)r.We must now distinguish between the following two sub-cases.
(B1) t * is very close to τ (r, x), i.e., t * ∼ τ (r, x).Furthermore, δu(y, x) ∼ y h , for all y ∈ [η, r], where η → 0 in the inviscid limit; hence, by using Eq. ( 7), we obtain Hence, the contribution from this case to T p col (r) is w 1 and w 2 are the weights associated with the respective types of interval collapse (see above).The first term on the right-hand side is the dominant one in the small-r limit.
(B2) t * is significantly smaller than τ , i.e., t * < τ , so case (A) follows after the time t * , i.e., where r 1 is the size of the interval when the shock forms.Then, by integrating Eq. ( 7) from t = 0 to t * , during which time δu(y, x) ∼ y h , we obtain By solving for r 1 from Eq. ( 13), substituting it in Eq. ( 12), expanding in powers of r, and, in the small-r limit, neglecting terms of order r h compared to unity, we obtain Hence, Here, we have incorporated the contributions from the two modes of interval collapse (cf.Eq. ( 11)), expanded the p-th power of τ from Eq. ( 14), and kept only the leading-order term in the small-r limit; w 3 and w 4 are the weights of the respective leading-order contributions.
We combine the contributions from cases (A), (B1), and (B2) [Eqs.(9), (11), and (15)] to obtain where W 1 , W 2 , W 3 and W 4 , the respective weights of these contributions, do not affect the values of the dynamic-scaling exponent, because, for any positive p, the second or last terms are the dominant ones (r 1).Consequently, the dynamic-scaling exponents, defined by are where we have set h = 1/3.Thus, the bifractal scaling of equal-time velocity structure functions [Eq.( 5)], in the d = 1 model of Burgers tuurbulence that we consider, is accompanied by bifractal dynamic scaling that is embodied in Eq. (18).

B. Dynamic Scaling: Direct Numerical Simulations
We now use data from our DNSs to calculate T p col (r); and we present log-log plots of it versus r, for 1 ≤ p ≤ 6, in Fig. 3(a).For all these values of p, T p col (r) shows power-law scaling over almost a decade and a half in r.From these power-law regions [in the blue-shaded rectangle in Fig. 3(a)], we extract the interval-collapse exponents z p col via a local-slope analysis of the curves.In Fig. 3(b), we plot z p col as a function of p, with pink and blue triangles for z p col from runs R1 and R2, respectively.In the inset of Fig. 3(b), we plot the local slopes, m(r), whose means and standard deviations across the blue-shaded portion of the inset yield z p col and its error bars, respectively.The black lines indicate the bifractal predictions of Eq. (18).Note that the values of z p col from runs R1 and R2 lie within error bars of each other; moreover, the exponents from the higher-resolution run R2 lie slightly closer to the bifractal-scaling prediction than those from run R1.The deviation from this prediction is most pronounced near the transition value, p = 3/2.
Our numerical results for z p col [Fig.3(b)] indicate dynamic multiscaling, much as the numerical results in Ref. [25] suggest multiscaling of the equal-time velocity structure function exponents ζ p .However, as suggested in Ref. [25], it behooves us to ask whether this is bona fide multiscaling or a numerical artifact.We remark that the small deviation of our DNS results for z p col from the bifractal-scaling predictions [Eq.(18)] are comparable to what has been observed [25] for the equal-time velocity-structure-function exponents ζ p .It is possible that this small deviation could decrease as we increase the resolution of our DNS.However, given the tiny difference between the results for runs R1 and R2, we might well have to go to DNS resolutions that are computationally prohibitive before our results for z p col approach the bifractal result of Eq. ( 18).Note that our DNS run R2 is, to date, the highest-resolution DNS of Eqs. ( 1) and (2).
We have also checked whether the interval-collapsetime exponents z p col depend on using different representative initial conditions, from the NESS of Eqs.
(1) and ( 2), into which we introduce N p equi-spaced tracers to start the runs with particles.We find that, given a large-enough spatial resolution, this dependence is small, and it lies within the errors bars that we have given for z p col .We show this explicitly in Fig. 8 in Appendix B with four representative DNSs with N = 2 16 collocation points.

C. Probability Distribution Functions (PDFs) of interval-collapse times
We turn now to the PDFs, Φ(τ col ), of τ col for different values of r; for notational convenience we suppress the argument r of τ ≡ τ col /T L in this subsection.

PDFs: Theoretical Considerations
We use the relation where P(v 0 ) is the PDF of v 0 ; for the Burgers equation with an external force that is limited to small k, this PDF is known to be a Gaussian [50].In our DNSs of Eqs. ( 1) and ( 2), we find this PDF to be a Gaussian too (see Appendix C), i.e., We now examine the forms of Φ(τ ) for the three cases (A), (B1), and (B2) that we have considered in Sec.

III A.
(A) There is a shock in the interval of size r at t = 0.The probability of finding such an interval is p r ∝ r.From Eq. ( 8), v 0 ∼ r/τ .Hence, for a given value of r, and in the small-r limit, (B) There are no shocks within the interval at t = 0, as in Sec.III A. The interval can collapse in one of the two ways illustrated in panel (B) of Fig. 2 and can occur with probability p s1 ≡ br or p s2 ≡ 1 − (a + b)r, depending on its mode of collapse (see above).We consider the following sub-cases: (B1) t * ∼ τ and v 0 ∼ r 1−h /τ , so, to leading order in r, in the limit r 1, (B2) Here, t * < τ and v 0 ∼ r/(τ − t * ), so, to leading order in r, for r 1, where the summation is over all possible values of t * , with weights w(t * ), and Θ, the Heaviside step function, ensures that, for a given value of t * , the contributions for all τ ≤ t * is zero.
By combining Eqs. ( 21)-( 23) and substituting h = 1/3, we get where µ 1 , µ 2 , and µ 3 are the weights of three contributions to Φ(τ ).We can recover the values z p col , given in Eq. ( 18), by calculating the moments of Φ(τ ) [Appendix D], and thus check our derivation of Φ(τ ).In summary, the contributions to the orderp moments of Φ(τ ), from the first two terms on the right-hand-side (RHS) of Eq. ( 24), scale as r 2p/3 and r p+1 , respectively.The leading-order contribution from the last term, Φ B2 , scales linearly with r and this comes from the binomial expansion of the term τ p after making the substitution, s = r/(τ − t * ).Note that Eq. ( 24) explicitly demonstrates that there is no unique dynamic scaling exponent; the arguments of the two exponentials suggest dynamic exponents of 2/3 and 1, respectively; whereas, the power-law factors multiplying the exponentials suggest dynamic exponents of 1/3, 1 and 1/2, respectively.Furthermore, this PDF is not a Gaussian.

PDFs: Direct Numerical Simulations (DNSs)
We now compare the analytical form of Φ(τ ) [Eq. (24)] with the results that we obtain from our DNSs for run R2.An accurate determination of the tails of Φ(τ ) is difficult because of binning errors and insufficient sampling of rare events in these tails.Hence, we use the rank-order method [25], which circumvents binning errors, to calculate the small-τ and large-τ limits of the cumulative PDF (CPDF), Q S (τ ), and the complementary CPDF (c-CPDF), Q L (τ ), of τ , respectively, which are defined as follows: In the τ → 0 limit, the contribution from Φ B2 (τ ) is negligible, so, for fixed r 1, Q S (τ ) scales as where erfc(x) is the complementary error function.In the second part of Eq. ( 26), we use the leading-order asymptotic form of erfc(x) in the limit of x → ∞.Furthermore, r 4/3 τ 2 log τ r 2/3 , so, in the limit of τ → 0, we get In Fig. 4(a), we plot |log 10 Q S | −1/2 versus τ /r 2/3 for different values of rN , with N the number of collocation points.In the inset, we plot |log 10 Q S | −1/2 versus τ for different values of rN .The plots in Fig. 4(a) collapse, fairly well, onto the dashed straight line [Eq.( 27)], thereby providing numerical support for our τ → 0 result for Q S (τ ).
In the large-τ limit, we expand the RHS of ( 24) in a Taylor series, retain only the leading-order contribution for r 1, and arrive at the following form for the c-CPDF: In Fig. 4(b), we present log-log plots of Q L (τ ) versus τ /r 2/3 for different values of rN ; the inset shows loglog plots of Q L (τ ) versus τ .We observe a reasonable collapse of the c-CPDFS onto a curve, which is close to, but steeper than, the straight dashed line that represents Eq. ( 28).This discrepancy arises because of the higher-order terms, which we have neglected in the Taylor expansion of the RHS of Eq. ( 24), and which have the same sign as that of the dominant term, so we expect the curves of Q L to be steeper for moderately large values of τ /r 2/3 .The corrections because of these higher-order terms reduce on increasing τ or decreasing r.Consequently, in Fig. 4(a), we observe that, for small values of r, the plots of Q L (τ ) (blue and red curves) tend to align with the dashed black line with slope −1, at large values of τ /r 2/3 ; this alignment is not so good as r increases.
In some cases they employ QL structure functions to uncover dynamic multiscaling, which is masked by the sweeping effect in the Eulerian framework.We show how to apply this QL approach to the study of dynamic scaling in the stochastically forced d = 1 Burgers equation [Eqs.(1) and ( 2)].
The order-p time-dependent structure function is defined as where δu(x, r, t) = |u(x + r, t) − u(x, t)| and S p (r) is the order-p equal-time structure function [Eq.( 5)b].By definition, F p (r, 0) = 1 for all r.From F p (r, t) we extract integral time scales of orders p and degree 1, as follows: where the second part defines χ I p , the order-p integraltime-scale exponent, by using the dynamic scaling Ansatz.The multifractal model of turbulence allows us to derive the following bridge relations [5][6][7] between χ I p and the equal-time exponents ζ p : If we use Eq. ( 5)c, we obtain: for p < 3 ; We calculate S p (r), F p (r, t) and T I p (r) from our DNS by using both the Eulerian velocity u(x, t) and the QL velocity where R(t) is the position of a Lagrangian particle, at time t, which was at position x 0 at t = 0 (see Appendix A for details ).In Fig. 5, we plot F Eu p and F QL p versus t, for different values of r, for the order p = 4.
We extract the Eulerian integral time scales, T I, Eu p (r), from Eq. ( 30) by integrating F Eu p (r, t) via the trapezoidal rule.The larger the time delay t, the more unreliable are the numerical data for these structure functions.Therefore, we use a time τ up as the upper limit of the integral, with τ up the time at which F Eu p decreases to a certain fraction λ.We consider both the short-and late-time estimates for T I, Eu Next we calculate the integral time scales with λ = 0.65.We observe that the corresponding dynamic exponents χ I, Eu p [green diamonds in Fig. 6(b)] lie close to, and within error bars of, unity for all p.These results remain unchanged for 0.6 < λ < 0.7.
Similar Eulerian results were obtained previously by Hayot and Jayaprakash [2].By using a different class of time-dependent structure functions, they defined a dynamic exponent z and observed that z is indeed unity at intermediate time scales, whereas, at short times, it is not unity.They suggested that this was a manifestation of the Taylor hypothesis or the sweeping effect: The advection of the small eddies by large ones yields a linear relation between spatial and temporal scales and thence a dynamic exponent z = 1.
We now calculate the QL structure functions F QL p , extract T I, QL p (r) [see the lower panel of Fig. 6(a)], and thence the exponents χ I, QL p , in the same way as we did above for their Eulerian counterparts.In particular, we use (a) 0.93 < λ < 0.96 and (b) 0.6 < λ < 0.7.In FIG. 5. Plots of (a) F Eu p (r, t) and (b) F QL p (r, t) versus t for p = 4 and rN/L = 100 (blue circles), 200 (red squares), 300 (yellow diamonds) and 400 (violet triangles).In both figures, the red and green shaded areas demarcate the regimes 0.93 < Fp(r, t) < 0.96 and 0.6 < Fp(r, t) < 0.7, respectively.Across each of these regimes, the integral-time-scale exponents remain unchanged (within numerical error bars).Thus, QL structure functions, which eliminate sweeping effects when we study dynamic multiscaling in incompressible turbulence [4-7, 20, 27, 51], do not remove sweeping effects in the d = 1 Burgers turbulence that we study.The qualitative reason for this is that tracers get trapped in shocks and, thereafter, they move with the shock.We conjecture that this inability of QL structure functions to remove sweeping effects holds in all models of compressible turbulence, once shocks are formed.

V. CONCLUSION AND DISCUSSIONS
We have carried out a detailed study of dynamic multiscaling in the NESS of the stochastically forced d = 1 Burgers equation [Eqs.( 1) and ( 2)].In particular, we have shown that there is not one but an infinite number of time-scales, associated with a single length scale, so this NESS displays dynamic multiscaling.We have proposed a theoretical framework to calculate these exponents.As mentioned earlier, to the best of our knowledge, this is the first time that such dynamic multiscaling exponents have been obtained analytically for turbulence in any nonlinear hydrodynamical PDE.We have validated our theoretical findings with the results of our high-resolution direct numerical simulations.We have also quantified the behaviors of the the tails of the CPDFs of τ col for different .Furthermore, we have explicitly shown that the QL methods [4-7, 20, 27], which have been used to study such multiscaling of turbulence in incompressible flows are not suitable for Burgers turbulence because it is compressible.
Our DNS results for z p col deviate slightly from the bifractal prediction [Eq.( 18)], which seems to indicate dynamic multiscaling, much like what is observed in Ref. [25] for the equal-time exponents ζ p .As suggested in Ref. [25], this deviation might possibly be a numerical artifact, rather than bona fide dynamic multiscaling, and could well decrease on increasing the resolution of our DNS by a few orders of magnitude.
It should be possible to use the multifractal model of turbulence to derive a bridge relation that links the collapse-time exponents z p col to the equal-time structure function exponents ζ p .The multifractal model has been used, earlier, to derive the integral-time bridge relation Eq. ( 31), from Eq. ( 30) in incompressible fluid turbulence [4-7, 20, 27, 51].
We outline a similar derivation for z p col in our d = 1 Burgers turbulence model in Appendix F, where we also discuss its limitations.
To generalize our study to the d−dimensional case, we consider the stochastically forced Burgers equation with the zero-mean, white-in-time stochastic force f (x, t).The covariance of its Fourier components satisfies The subscripts i and j denote the Cartesian components of the force.This force must be chosen to be curl free to ensure that the flow is not vortical.The choice β = d yields a K41-type inertial-range energy spectrum at the level of a one-loop RG.In the turbulent NESS for this model [Eqs.(34) and (35)] shocks are extended structures (lines if d = 2 and surfaces if d = 3).
For introductions to such shocks, we refer the reader to Refs.[22,52].The natural generalization of the interval-collapse times, which we have used above for d = 1 Burgers turbulence, are the collapse times of Lagrangian d−simplices, i.e., triangles and tetrahedra in d = 2 and d = 3, respectively, with tracers at their vertices.We have carried out preliminary DNSs to study dynamic multiscaling in the NESS of the d = 2 stochastically forced Burgers equation.We find that we must distinguish between the collapse of a triangle to one with vanishing area, in one of the following ways: (i) all three vertices collapse to a point; (ii ) any two vertices come together; or (iii) all three vertices become colinear, with none of them overlapping.We have verified that the collapse times for cases (i), (ii), and (iii) capture the dynamic multiscaling of the NESS of Eqs. ( 34) and ( 35), but with dynamic exponents that depend on which of the cases [(i)-(iii)] we consider.The collapse of Lagrangian tetrahedra in d = 3 requires even more cases than that for triangles in d = 2.The elucidation of the dependence of these exponents on initial conditions in this NESS is a challenge because it requires high-resolution DNSs (for d = 1 Burgers turbulence see Fig. 8 in Appendix B).The details of our studies of dynamic multiscaling in the NESSs of Eqs. ( 34) and ( 35), for d = 2 and d = 3, will be given elsewhere.
Our study has obvious implications for the elucidations of dynamic scaling in compressible turbulence (CT) with shocks.We will explore this in future work.We note that the statistics of velocity fluctuations in the stochastically forced d = 1 Burgers equation that we consider is similar to that in the fully detailed d = 1 model of CT with unit Mach number, investigated in Ref. [53].Of course, the study of CT with shocks is extremely relevant in the engineering of combustion systems and supersonic jets and in astrophysical studies of, e.g., the solar wind, interstellar gases, and molecular clouds [54][55][56][57][58][59][60][61][62][63][64][65][66].The application of our interval-collapse times, and their generalizations to d > 1 hold promise for studies of the spatiotemporal statistics of the compressible turbulent flows mentioned above and those that are obtained in DNSs of compressible hydrodynamical equations.It would be interesting to explore whether times such as τ col , and their generalization to heavy inertial particles, could be related to characteristic times in processes that require the agglomeration of masses or chemical species advected by highly compressible turbulent flows in astrophysical turbulence (e.g., dust in molecular clouds) and in turbulent combustion [67,68].In fact, in Refs.[67,68], the authors have looked at the possible effects of clustering in weakly compressible flows.Similar studies for clustering because of shocks will be very interesting.
It would also be fascinating to explore, in detail, the relation of our work with statistical-mechanical studies, carried out in the context of the KPZ equation, which investigate the behavior of passive random walkers following the wrinkles of a rough surface whose fluctuations are determined by this equation [69][70][71][72].(B1) The contribution from the second term on the RHS of Eq. ( 24) is By making the substitution y = s 2 /2, the integrals over each term in Eq. (D11) can be be written in terms of the lower incomplete gamma function γ(n, x) = ´x 0 t n−1 e −t dt.We use the asymptotic form of γ(n, x), in the limit of x → 0, γ(n, x) ∼ x n n , because r 1.We find that the leading-order contribution from I 1 scales as r 2 .By making the same substitution in Eq. (D12), we find that each one of the integrals in Eq. (D12) can be written in terms of the upper incomplete gamma function Γ(n, x) = ´∞ x t n−1 e −t dt.By using the asymptotic properties of the resulting terms in the limit of r 1, we find that the r−dependent leading-order contribution scales linearly with r, with logarithmic corrections.As a result,  8), (10), and ( 14), and T L ∼ L/v 0 , τ ≡ τ col ( )/T L , r ≡ /L, τ 1 for all collapsing intervals (recall that t * < τ ), C 1 C 2 for all p > 0, whence we have In this way, we are able to extract from P (τ col ) the same scaling behavior of the moments of τ col as we had obtained from our theoretical arguments given in the main paper.

FIG. 3 .
FIG. 3. (a)Log-log plots of T p col (r) versus r for p = 1 (dark blue circles), p = 2 (red squares), p = 3 (yellow diamonds), p = 4 (violet inverted triangles), p = 5 (green pentagrams), and p = 6 (light blue triangles) from run R2; we extract the exponents z p col by calculating the local slopes of these graphs across the shaded region.(b) The exponent z p col versus p from runs R1 (pink triangles) and R2 (blue triangles).The solid line shows our theoretical prediction, Eq. (18).Inset: The local slope m(r) as a function of r, the color and symbol coding being the same as that in (a).The means and standard deviations of m(r), calculated over the shaded region, yield the exponents z p col and their error bars, respectively.

p
(r) by employing large and small values of λ, respectively.For λ = 0.95 we plot T I, Eu p (r) for different values of p in the upper panel of Fig. 6(a); and from the local slopes of these graphs, we evaluate the dynamic exponents χ I, Eu p and display them via blue triangles in Fig. 6(b).These exponents closely match with the multifractal-bridge-relation prediction, which we derive by substituting the Eulerian equal-time structure function exponents ζ Eu p (see Appendix E) into Eq.(31) and which is indicated by the shaded region in Fig. 6(b).The width of this shaded region indicates the error bars in the equal-time exponents that are used in the bridge relation.The values of χ I, Eu p remain unchanged, within error bars, for 0.93 < λ < 0.96.

FIG. 6 .
FIG. 6. Plots for run R2: (a) Log-log plots of the integral time scales, T I, Eu p (r) and T I, QL p (r), versus r.In both panels, the blue circles and red squares correspond to orders p = 2 and p = 3, respectively, with λ = 0.95 (see text); points obtained with λ = 0.65 are indicted by yellow diamonds and green triangles for orders p = 2 and p = 3, respectively.The pink shaded area indicates the regime in which we carry out local-slope analyses.(b) Blue and violet triangles are the integral-time-scale exponents, χ I, QL p

1 FIG. 9 .
FIG. 9. Plot of the normalized c-CPDF QL(|v0|) of the magnitude of v0 versus |v0|, represented in blue; the black dashed curve is the c-CPDF of a standard normal distribution.The two curves match closely with each other, thereby justifying our assumption in the main text.