The universality of dynamic multiscaling in homogeneous, isotropic Navier–Stokes and passive-scalar turbulence

We systematize the study of dynamic multiscaling of time-dependent structure functions in different models of passive-scalar and fluid turbulence. We show that, by suitably normalizing these structure functions, we can eliminate their dependence on the origin of time at which we start our measurements and that these normalized structure functions yield the same linear bridge relations that relate the dynamic-multiscaling and equal-time exponents for statistically steady turbulence. We show analytically, for both the Kraichnan model of passive-scalar turbulence and its shell model analogue, and numerically, for the Gledzer–Ohkitani–Yamada (GOY) shell model of fluid turbulence and a shell model for passive-scalar turbulence, that these exponents and bridge relations are the same for statistically steady and decaying turbulence. Thus, we provide strong evidence for dynamic universality, i.e. dynamic-multiscaling exponents do not depend on whether the turbulence decays or is statistically steady.


Introduction
The elucidation of the universal scaling properties of equal-time and time-dependent correlation functions in the vicinity of a critical point was one of the most important achievements of statistical mechanics over the past forty years. The analogous systematization of the power laws and associated exponents that govern the behaviours of structure functions in a turbulent fluid, or in a passive-scalar advected by such a fluid, is a major challenge in the areas of nonequilibrium statistical mechanics, fluid mechanics and nonlinear dynamics. The power-law behaviour of equal-time structure functions has been studied in detail over the past few decades [1]; and, especially in the case of passive-scalar turbulence [2], significant progress has been made in understanding the multiscaling of equal-time structure functions. The nature of multiscaling of time-dependent structure functions has been examined recently [3]- [7] but only for the case of statistically steady turbulence. We develop here the systematics of the multiscaling of timedependent structure functions for the case of decaying fluid and passive-scalar turbulence 6 .
To set the stage for our discussion of time-dependent structure functions in turbulence, it is useful to begin by recalling some well-known results from critical phenomena [9,10]. At a 3 critical point for a spin system in d-dimensions, the equal-time, two-spin correlation function g and its spatial Fourier transformg assume the following power-law scaling forms: Heret ≡ (|T − T c |)/T c , T and T c are the temperature and the critical temperature, respectively, h ≡ H/k B T c , H is the external field, k B is the Boltzmann constant, the spins are separated by the vector r[r = |r|], k is the wavevector, k = |k|, ν, , η are critical exponents, and G andG are scaling functions. Away from the critical point such correlation functions decay exponentially; the associated correlation length ξ c diverges in the vicinity of the critical point; e.g. as ξ c ∼t −ν , if h = 0. Time-dependent correlation functions also assume scaling forms in the vicinity of the critical point and the characteristic relaxation time τ diverges as suggested by the dynamicscaling Ansatz [9] τ ∼ ξ z c , which introduces the dynamic-scaling exponent z.
The generalization of such a dynamic-scaling Ansatz to the case of homogeneous, isotropic turbulence is our prime concern here. The power-law behaviours of equal-time structure functions, in the inertial range (to be defined later), in turbulence are reminiscent of the algebraic dependence of critical-point correlation functions on r . However, there are important differences between the two that must be appreciated before we embark on a systematization of time-dependent structure functions in turbulence. We begin with the increments of the longitudinal component of the velocity δu (x, r, t) ≡ [u(x + r, t) − u(x, t)] · (r/r ) and passive-scalar δθ (x, r, t) ≡ [θ(x + r, t) − θ(x, t)], respectively; here u(x, t) and θ(x, t) denote, respectively, the velocity of the fluid and the passive-scalar density at the point x and time t, and the subscript the longitudinal component. The order-p, equal-time structure functions, for the fluid (superscript u) and passive-scalar (superscript θ ) fields, are defined as follows: the angular brackets indicate averages over the steady state for statistically steady turbulence or over statistically independent initial configurations for decaying turbulence; for stochastic differential equations, like the Kraichnan model (section 2.1), the angular brackets denote an average over the statistics of the noise; and the power laws, characterized by the equal-time exponents ζ u p and ζ θ p , hold for separations r in the inertial range η d r L, where η d is the Kolmogorov dissipation scale and L the large length scale at which energy is injected into the system.
Kolmogorov's phenomenological theory [1,11,12] of 1941 (K41) suggests simple scaling, with ζ u,K41 p = p/3, but experimental and numerical evidence favours equal-time multiscaling with ζ u p and ζ θ p nonlinear, convex, monotone-increasing functions of p. For the simplified Kraichnan model [2], [13]- [15] of passive-scalar turbulence (section 2.1) multiscaling of equaltime structure functions can be demonstrated analytically in certain limits. The analogue of the 4 K41 theory for passive-scalar turbulence is due to Obukhov and Corrsin [16,17]; if the Schmidt number Sc ≡ ν/κ 1, then their theory yields K41 exponents for the passive-scalar case; here ν is the kinematic viscosity of the fluid and κ is the diffusivity of the passive scalar.
We now turn our attention to the order-p, time-dependent structure functions, F u p (r, {t 1 , . . . , t p }) ≡ [δu (x, r, t 1 ) · · · δu (x, r, t p )] and F θ p (r, {t 1 , . . . , t p }) ≡ [δθ (x, r, t 1 ) . . . δθ (x, r, t p )] , and the dynamic scaling exponents associated with the time scales which can be extracted from such structure functions. We will discuss these time scales and their scaling exponents in detail in sections 2 and 4. A straight-forward extension of simple, K41 scaling to time-dependent structure functions implies that the dynamic exponents z K41 p = 2/3 for all p. This naïve extension fails for two reasons: (a) it does not distinguish between the temporal behaviours of structure functions of Eulerian, Lagrangian, and quasi-Lagrangian (section 2) velocities or passive-scalar densities; and (b) it does not account for the multiscaling of structure functions. These difficulties have been overcome to a large extent for statistically steady turbulence [3]- [7], [18] as we summarize below. There is consensus now that Eulerian structure functions display simple scaling with only one dynamicscaling exponent z E = 1 because of the sweeping effect: the mean flow, or the flow caused by the largest eddy, advects small eddies, so spatial separations r in (3) are related linearly to temporal separations τ via the mean-flow velocity [18]. By contrast, it is expected that Lagrangian [18] or quasi-Lagrangian [3]- [7], [19] time-dependent structure functions should show nontrivial dynamic multiscaling. The task of extracting well-averaged time-dependent Lagrangian or quasi-Lagrangian structure functions from a direct numerical simulation (DNS) of the Navier-Stokes (NS) equation is a daunting one [20]: a dynamic exponent has been extracted from a full Lagrangian study [18] only for order p = 2. Thus the elucidation of dynamic multiscaling has relied on predictions based on generalizations of the multifractal formalism [3]- [7] and on numerical studies of shell models [4]- [7]. These studies show that, if dynamic multiscaling exists, time-dependent structure functions must be characterized by an infinity of time scales and associated dynamic multiscaling exponents [6]. Furthermore, the dynamic exponents depend on how we extract time scales from time-dependent structure functions; e.g. for fluid turbulence, time scales obtained from integrals (superscript I and subscript 1) and second derivatives (superscript D and subscript 2) of order-p time-dependent structure functions yield the different dynamic exponents z I,u p,1 and z D,u p,2 . Finally, the different dynamic multiscaling exponents are related by different classes of linear bridge relations to the equal-time multiscaling exponents. For a careful discussion of these issues we must of course define time-dependent structure functions. The details necessary for this paper are given in sections 2 and 4.
The dynamic multiscaling of time-dependent structure functions described briefly above applies to statistically steady turbulence. Does it have an analogue in the case of decaying turbulence, since time-dependent structure functions must, in this case, depend on the origin of time t 0 at which we start our measurements? This question has not been addressed hitherto. Weshowherehowtoansweritindecayingfluidandpassive-scalarturbulence 6 .Inparticular, we propose suitable normalizations of time-dependent structure functions that eliminate their dependence on t 0 ; we demonstrate this analytically for the Kraichnan version of the passivescalar problem and its shell-model analogue and numerically for the GOY shell model [1,21,22] for fluids and a shell-model version of the advection-diffusion (AD) equation. In these models we then analyse the normalized time-dependent structure functions for the case of decaying turbulence like their statistically steady counterparts [6,7]. This 5 requires a generalization of the multifractal formalism [1] that finally yields the same bridge relations between dynamic and equal-time multiscaling exponents as for statistically steady turbulence [6,7]. For the Kraichnan version of the passive-scalar problem we show analytically that simple dynamic scaling is obtained. This is because (section 3) the advecting velocity is random and white in time. In addition, we find numerically for shell models of fluid and passive-scalar turbulence that dynamic-multiscaling exponents have the same values for both statistically steady and decaying turbulence; so, in this sense, we have universality of the multiscaling of time-dependent structure functions in turbulence. The equal-time analogue of this universality has been discussed in [23].
The remaining part of this paper is organized as follows. In section 2, we introduce the models we use and give the details of our numerical simulations. Section 3 presents our analytical studies of decaying turbulence in the Kraichnan model and its shell-model analogue. Section 4 shows how to generalize the multifractal formalism to allow for time-dependent structure functions in decaying turbulence and how to obtain bridge relations between dynamic and equal-time multiscaling exponents in this case. In section 5 we present the results of our numerical studies of dynamic multiscaling in the GOY shell model for fluid turbulence and for shell models of a passive-scalar field advected by a turbulent velocity field. Section 6 is devoted to a discussion of our results in the context of earlier studies; we also suggest possible experimental tests of our predictions.

Models and numerical simulations
We have used several models to study time-dependent structure functions in fluid and passivescalar turbulence. These range from the NS and AD equations to simple shell models; the latter are well suited for our extensive numerical studies. It is useful to begin with a systematic description of these models.
Fluid flows are governed by the NS equation (4) for the velocity field u(x, t) at point x and time t, augmented by the incompressibility constraint (5), since we restrict ourselves to low Mach numbers: ∇.u = 0.
Here ν 0 is the kinematic viscosity, P the pressure, the density ρ is taken to be 1, and f the external force, which is absent when we consider decaying turbulence. If and v are, respectively, characteristic length and velocity scales of the flow, the Reynolds number Re ≡ ( v/ν 0 ) provides a dimensionless measure of the strength of the nonlinear term in (4) relative to the viscous term; for the case of decaying turbulence it is convenient to use the Reynolds number for the initial state, i.e. Re with v the root-mean-square (rms) velocity of the initial condition and the system size (the linear size of the simulation box in a DNS). Given the incompressibility condition (5), the pressure can be eliminated from (4) and related to the velocity by a Poisson equation. The equation for the velocity alone is most easily written in terms of the spatial Fourier transformũ(k, t) of u(x, t) and it can be shown easily thatũ(k, t) is affected directly by all other Fourier modes. This is the mathematical representation of the sweeping effect in which the largest eddies (i.e. modes with small k ≡ |k|) directly advect the smallest eddies (i.e. large-k modes); such direct sweeping lies at the heart of the Taylor hypothesis [24] and leads 6 eventually to trivial dynamic scaling for time-dependent structure functions of Eulerian fields with a dynamic exponent z E = 1 for fluid turbulence. As we have mentioned above, nontrivial dynamic multiscaling is expected if we use Lagrangian or quasi-Lagrangian velocities. The Lagrangian formulation is well known [25]; the quasi-Lagrangian [3,19] one uses the following transformation for any Eulerian field ψ(x, t): whereψ is the quasi-Lagrangian field and R(t; r 0 , 0) is the position at time t of a Lagrangian particle that was at point r 0 at time t = 0. The AD equation for the Eulerian passive-scalar field θ(x, t) is ∂θ ∂t where κ is the passive-scalar diffusivity and, if we consider decaying passive-scalar turbulence, the external force f θ is set to zero. The advecting velocity field u should be obtained, in principle, by solving equations (4) and (5). By using equation (6) we get the quasi-Lagrangian version of the AD equation (7): Direct numerical simulations of equations (4) and (5) or equation (7), though feasible, have not yet provided data that are averaged well enough to yield reliable time-dependent structure functions of quasi-Lagrangian velocity [20] or passive-scalar fields. Time-dependent Lagrangian structure functions have been obtained [18] only for order p = 2. Thus a firstprinciples DNS study of dynamic multiscaling in fluid or passive-scalar turbulence is not possible at the moment. However, significant progress has been made in statistically steady turbulence by studying dynamic multiscaling in simplified models like the Kraichnan model for passive-scalar turbulence and shell models for fluid and passive-scalar turbulence. We discuss these models below since our studies of decaying turbulence will be based on them.

The Kraichnan model (model A)
The Kraichnan model for passive-scalar turbulence [2], [13]- [15] begins with the AD equation (7) but replaces the NS velocity field by one in which each component u i (x, t) of the velocity is a zero-mean, delta-correlated, Gaussian random variable with the covariance The Fourier transform of D i j (r) has the form where q is the wavevector, L the characteristic large length scale, η d the dissipation scale, d the dimension of the system, and ξ a tunable parameter. In the limits L → ∞ and η d → 0, of relevance to turbulence, we have, in real space, with where D 0 ∼ C 1 L ξ ; C 1 and D 1 are dimensional constants. We refer to (7)-(12) as model A to distinguish it from other models that we use. For 0 < ξ < 2, this model shows multiscaling of order-p, equal-time passive-scalar structure functions, as can be shown analytically, in certain limits [2]. However, for the case of statistically steady turbulence, this model exhibits simple dynamic scaling [6].

Shell models
We will also use some shell models for fluid and passive-scalar turbulence. These models are highly simplified representations of the NS or the AD equations (4)- (7) and are, therefore, far more tractable numerically than (4)- (7). Nevertheless, shell models retain enough properties of their parent equations to make them useful testing grounds for the multiscaling of structure functions in turbulence. Shell models are defined on a logarithmically discretized Fourier space in which complex scalar variables (e.g. the velocity u n or passive scalar θ n ) are associated with the shells n and scalar wavevectors k n = k 0 λ n ; typically λ = 2, k 0 = 1/16; and the boundary conditions are that the shell variables vanish if n < 1 or n > N (we use N = 22). These models consist of coupled, nonlinear, ordinary differential equations (ODEs) that specify the temporal evolution of the shell variables u n and θ n . Shell-model ODEs are similar to the Fourier-space versions of their parent partial differential equations: (a) their dissipative terms are linear in one of the shell variables and quadratic in k n ; (b) their analogues of advection terms are linear in k n and bilinear in the shell variables; e.g. for fluid turbulence a representative term is of the form ik n u n u n , with n = n ; and (c) they conserve the shell-model analogues of the energy, helicity, etc, in the absence of dissipation and forcing. However, variables in a given shell are influenced directly only by their nearest-and next-nearest-neighbour shell variables; by contrast, Fourier transformations of the NS and the AD equations couple every Fourier mode to every other Fourier mode, leading to the sweeping effect mentioned above. Thus direct sweeping is absent in shell models, so they are often thought of as an approximate, quasi-Lagrangian representation of their parent equations.
For studies of decaying turbulence one can envisage several initial conditions. We have used initial conditions of two types: (a) in the first (type-I), we drive the system to a statistically steady turbulent state by forcing the first shell (n = 1); we then turn off the force and allow the turbulent state and the associated energy spectrum to decay freely; our measurements are made in this decaying state; and (b) for the case of fluid turbulence we use a second initial condition (type-II) in which all the energy is concentrated in the first few shells with small k n , i.e. large length scales; we then allow the system to evolve without any force; the energy cascades to large values of k n till the energy spectrum becomes similar to that in forced turbulence; this spectrum then decays slowly with time and the measurements we report are made during this stage of evolution.

Model B
We use the shell-model analogue of the Kraichnan model introduced in [26] in which the equation for the passive-scalar variable θ n is where the asterisks denote complex conjugation, a n = k n /2, b n = −k n−1 /2 and c n = k n+1 /2; f n is an additive force that is used to drive the system to a steady state; the boundary conditions are The advecting velocity variables are taken to be zero-mean, white-in-time, Gaussian random complex variables with covariance where C 2 is a dimensional constant. We refer to equations (13) and (14) as model B.
In our numerical simulations of this model we first obtain a statistically steady turbulent state by forcing the first shell with a random, Gaussian, white-in-time force and by using initial conditions (superscript 0): θ 0 n = k 1/2 n e iϑ n , for n = 1, 2 and u 0 n = 0 for 3 n N ; ϑ n is a random phase angle distributed uniformly between 0 and 2π. The force is then switched off and measurements are made as the turbulence decays. We use a weak, order-one, Euler scheme to integrate the resulting Ito form [26] of (13) with an integration time step δt = 2 −24 , number of shells N = 22, diffusivity κ = 2 −14 and ξ = 0.6.
For such a passive-scalar shell model the order-p, equal-time, structure function and its exponent are defined via it is natural, therefore, to define the time-dependent version of S θ p (k n ) as follows: The power-law dependence on the right-hand side of (15) is obtained for k n in the inertial range.
In our numerical calculations we use extended self-similarity (ESS) to extract the exponent ratios ζ θ p /ζ θ 2 .

Model C
The most commonly used shell-model analogue of the NS equation is the GOY model [1,21,22]: d dt + νk 2 n u n = i a n u n+1 u n+2 + b n u n−1 u n+1 + c n u n−1 u n−2 * + f n .
The coefficients a n = k n , b n = −δk n−1 and c n = −(1 − δ)k n−2 are chosen in a manner that conserves the shell-model analogues of energy and helicity in the inviscid, unforced limit; an external force f n drives the system to a steady state. We use the standard value δ = 1/2; the boundary conditions are u −1 = u 0 = 0 and u N +1 = u N +2 = 0. We will refer to this as model C.
We use two different kinds of initial conditions in our study of decaying fluid turbulence in this model. For type-I initial conditions we first drive the system to a statistically steady turbulent state with an external force f n = (1 + i) × 5 × 10 −3 δ n,1 . The force is then switched off and the shell velocities at this instant are taken as the initial condition. The turbulence then decays. Our structure-function measurements are made during this period of decay. To obtain the second type of initial condition, the energy is initially concentrated in the first few shells by choosing the following initial (superscript 0) velocities: u 0 n = k 1/2 n e iϑ n , for n = 1, 2 and u 0 n = k 1/2 n e −k n 2 e iϑ n , for 3 n N , with ϑ n a random phase angle distributed uniformly between 0 and 2π . This energy then cascades down the inertial-range scales without significant dissipation until it reaches dissipation-range scales at cascade completion. The energy dissipation-rate per unit mass shows a peak, as a function of time, roughly at cascade completion [9,27], and the energy 9 spectrum E(k) and the structure function (18) show well-developed inertial ranges. These decay very slowly in time, so at each instant exponents can be determined from plots of E(k) and the structure functions. Therefore, for initial conditions of this type, we wait for cascade completion before making measurements of structure functions.
We employ the slaved Adams-Bashforth scheme [28,29] to integrate the GOY-model equations with a time step δt = 10 −4 . In our numerical simulations the viscosity ν = 10 −7 and the total number of shells N = 22; this provides us with a large inertial range from which exponents can be obtained reliably.
For this model, the order-p, equal-time structure function and its exponent are defined as follows: the associated time-dependent structure function is where the power-law dependence on the right-hand side of (18) holds for k n in the inertial range.
(For statistically steady turbulence, the time-dependent structure function has no dependence on t 0 , so, without loss of generality, t 0 can be taken to be 0.) A direct determination of ζ u p from (18) is not very accurate because of an underlying 3-cycle in the static version of the GOY shell model [30]. The effects of this 3-cycle can be filtered out to a large extent by using the modified structure function u we use u p (n) in our numerical calculation of ζ u p . We recall that (20) is obtained from the more general form of the structure function (18) by studying the scaling of the triple moments of the velocity increments which are free from the period-3 oscillations. We refer the reader to [30] for a complete derivation of (20). We measure time in terms of the initial large eddy-turnover time t L ≡ 1/(u rms k 1 ) and the initial rms velocity u rms ≡ n |u 0 n | 2 1/2 .

Model D
A turbulent velocity field does not have the simple statistical properties assumed in models A and B. To overcome this we study the shell model of [31], hereafter referred to as model D, in which the advecting velocity u n is a solution of the GOY shell model (17). The passive-scalar shell variables θ n obey where a n = k n , b n = −k n−1 /2 and c n = −k n+1 /2; f n = (1 + i) × 5 × 10 −3 δ n,1 is an additive force that drives the system to a steady state; the boundary conditions are For this model, we start with type-I initial conditions, i.e. we force both the coupled equations (17) and (21) till a statistically steady turbulent state is obtained and then switch off the force. The shell variables at this instant of time are taken as the initial condition; and then the turbulence is allowed to decay.
We employ a second-order Adams-Bashforth scheme to integrate model D with a time step δt = 10 −4 and set the diffusivity κ = 5 × 10 −7 so the Schmidt number ν/κ = 1/5; and N = 22 as in model C. The definitions of structure functions for model D are the same as those for model B, i.e. equations (15) and (16). We limit the effects of the 3-cycle (mentioned above for model C) in our numerical evaluations of ζ θ p by using the modified structure function θ p (n), which is defined as

Analytical results for models A and B
In this section, we present our analytical results for time-dependent, passive-scalar structure functions for the Kraichnan model (model A) and its shell-model analogue model B. For model A, we obtain results for both Eulerian and quasi-Lagrangian structure functions. We find, in particular, that time-dependent structure functions for these models can be factorized into a part that depends on the time origin t 0 and a part that depends on t but is independent of t 0 . This important result motivates a similar factorization hypothesis that we propose, and verify numerically, for models C and D in subsequent sections.

Model A
Consider first the Eulerian version of the AD equation (7). We assume that a turbulent statistical steady state has been established because of an external force f θ . We turn off this force at time 0. A spatial Fourier transform of (7) now yields where the tildes denote spatial Fourier transforms, we sum over repeated indices, and the statistics ofũ j (q, t) are specified by equations (9) and (10). The Fourier-space, second-order correlation functionF θ 2 (k, t 0 , t) ≡ θ (−k, t 0 )θ(k, t 0 + t) (t 0 0 is any time origin) satisfies the equation which can be combined with (23) to get We average over the statistics of the advecting velocity field (9) by using Novikov's theorem [32] 7 e.g. the first term in (25) reduces to Finally equations (9), (23) and (24) yield A spatial Fourier transform allows us to integrate this equation to obtaiñ If we set t = 0, we see thatF θ 2 (k, t 0 , 0) is just the passive-scalar, equal-time structure function; ϕ θ 2 (k, t 0 ) is, therefore, proportional to the Fourier transform of the passive-scalar equal-time structure function S θ 2 (r ). Thus, for a fixed but large value of L, we get, in the Eulerian framework, simple dynamic scaling with an exponent z E 2 = 2. We note that we get an exponent which is not equal to unity (section 1) in this case because of the white-in-time nature of the advecting field.
The analogue of equation (28) for the second-order, quasi-Lagrangian structure function follows from [7] the quasi-Lagrangian form of the AD equation (8).
By substituting for d i j from (12) we obtain, for the isotropic case, A spatial Fourier transformation allows us to integrate this equation to obtain, in the limits κ → 0, η → 0 and L → ∞, ϕˆθ 2 (k, t 0 ) is now proportional to the Fourier transform of the equal-time, quasi-Lagrangian, passive-scalar structure function (which is, of course, the same as the equal-time, Eulerian, passive-scalar structure function). Equation (32) shows that, in the quasi-Lagrangian framework, Fθ 2 (k, t 0 , t) factorizes into a part that depends on t 0 and another which depends only on t.
From the second factor we get simple dynamic scaling with an exponent z 2 = 2 − ξ , which is different from the Eulerian exponent z E 2 = 2 for model A. Such a factorization should also follow for higher-order, time-dependent, structure functions as we show in the next subsection for model B.

Model B
We can use the methods of the previous subsection to obtain analytical expressions for timedependent structure functions for model B. Consider first the second-order structure function whence The angular brackets denote an average over the statistics of u n (t) that are specified by equation (14). By using the complex conjugate of (13) in (34) and Novikov's theorem we get terms of the form Finally by using (14) we obtain where C 2 is a dimensional constant. Integration now yields with A(ξ ) = (2 (2ξ −2) + 2 −(2ξ −2) ) + (2 ξ + 2 −ξ ) + (2 (ξ −2) + 2 −(ξ −2) ). Similarly, we obtain the following exact expression for the fourth-order structure function: φ θ 2 (n, t 0 ) and φ θ 4 (n, t 0 ) are, respectively, the second-and fourth-order, equal-time, quasi-Lagrangian, passive-scalar structure functions. Thus z 2 = z 4 = 2 − ξ , as for the quasi-Lagrangian structure functions of model A. (Recall that we expect quasi-Lagrangian behaviour for shell models since they do not have a direct sweeping effect.) The equality of z 2 and z 4 indicates that we have simple dynamic scaling in model B. We expect all the quasi-Lagrangian exponents z p to be 2 − ξ for this model, but the analytical demonstration of this result becomes more and more complex with increasing p.
Given a factorization of the form shown in (37), it is possible to normalize the timedependent, passive-scalar structure function F θ p (n, t 0 , t) by its value at t = 0 and thus make it independent of t 0 . We cannot prove that such a factorization exists in models C and D, but we present compelling numerical evidence for it in section 6.

Multifractal formalism for models C and D
Equal-time Eulerian and quasi-Lagrangian structure functions are the same for homogeneous and isotropic turbulence [33]. Since quasi-Lagrangian structure functions are required for our study of dynamic multiscaling, we present the multifractal formalism in terms of quasi-Lagrangian variables. Multiscaling in statistically steady fluid turbulence can be rationalized by using the multifractal formalism, which assumes that a turbulent flow has a continuous set of scaling exponents h in the set I ≡ (h min , h max ), instead of a single exponent (e.g. h = 1/3 yields simple K41 scaling) [1]. The scaling exponents h characterize the behaviour of velocity differences δû r (x): for each h ∈ I there exists a set h ⊂ R 3 of fractal dimension Dû(h) and δû r (x)/û L ∼ (r/L) h for separations r in the inertial range if x ∈ h and withû L the velocity at the forcing scale L. Given the measure dµ(h) for the weights of the different values of h, the order-p, equal-time velocity structure function is where the ph term comes from p factors of (r/L) h and the additional factor of (r/L) 3−Dû (h) is the probability of being within a distance ∼ r of the set h , of dimension Dû(h), which is embedded in three dimensions. Dû(h), h min and h max are assumed to be universal. In the limit r/L → 0, of relevance to fully developed turbulence, we get the equal-time scaling exponent ζˆu p = inf h [ph + 3 − Dû(h)] by the method of steepest descents. We now define the order-p, time-dependent, quasi-Lagrangian velocity structure function Fû p (r, {t 1 , . . . , t p }) ≡ [δû (x, r, t 1 ) · · · δû (x, r, t p )] .
For simplicity, we consider t 1 = t and t 2 = · · · = t p = 0, denote the structure function by Fû p (r, t), and suppress the subscript , i.e. we use δû r (x) ≡ δû (x, r, t 1 ). The natural extension of the multifractal formalism to the case of time-dependent structure functions in statistically steady fluid turbulence follows from the Ansatz [6] Fû p (r, t) where the scaling function G p,h ((t/τ p,h )) is assumed to have a characteristic decay time τ p,h ∼ r/δû r (x) ∼ r 1−h and G p,h (0) = 1. For statistically steady passive-scalar turbulence the application of the multifractal formalism is more complicated than it is for fluid turbulence. This is because we must now deal with a joint multifractal distribution of both the velocity and passive-scalar variables. Therefore, the order-( p, p ), equal-time structure function for passive-scalar turbulence has the multifractal representation whereû andθ are assumed to possess a range of universal scaling exponents h ∈ I ≡ (h min , h max ) and g ∈ I ≡ (g min , g max ), respectively. For each pair of h and g in these ranges, there exists a set h,g ⊂ R 3 of fractal dimension Dû ,θ (h, g). The increments in the velocity δû r (x) and the passive-scalar field δθ r (x) scale as δû r (x)/û L ∼ (r/L) h and δθ r (x)/θ L ∼ (r/L) g for separations r in the inertial range if x ∈ h,g .û L andθ L are, respectively, the velocity and the passive-scalar variables at the forcing scale L. For simplicity we will only consider passivescalar structure functions with p = 0 in equation (42), i.e.
which have the multifractal representation as before, the equal-time exponents ζˆθ p can be related to Dû ,θ (h, g) by the method of steepest descents. The corresponding order-p, time-dependent passive-scalar structure function is As in (40), we consider t 1 = t and t 2 = · · · = t p = 0 for simplicity. Given the nature of multiscaling in passive scalars advected by a turbulent velocity field, it is possible to understand passive-scalar turbulence within the framework of the multifractal formalism. However, the analogous expression for time-dependent structure functions in passive-scalar turbulence has to take into account the multifractal nature of the advecting velocity field [7]. We generalize the multifractal representation of time-dependent structure functions in the following way: and we assume that (a) the function G p,h,g (t/τ phg ) has a characteristic decay time τ phg , (b) G p,h,g (0) = 1 and (c) that the dominant contribution to τ phg , in the limits κ → 0, η d → 0 and L → ∞, has the following scaling form Given Fû p (r, t) and Fθ p (r, t), we define the order-p, degree-M, integral-(I ) and derivative-(D) times as follows [6,7] (if the integrals and derivatives in equations (49) and respectively, for r in the inertial range. By substituting the multifractal form (41) in (51), evaluating the time integral first, and then performing the integration over the multifractal measure by the saddle-point method, we obtain the integral bridge relations (see appendix) which were first obtained in [3]. Similarly we get the derivative bridge relations which were first obtained in forced Burgers turbulence [34,35] for the special cases (a) p = 2, M = 1 and (b) p = 2, M = 2, respectively, and for statistically steady fluid turbulence in [6]. We note that within the K41 phenomenology (ζ u,K41 p = p/3) the bridge relations yield the same dynamic exponent z K41 p = 2/3 for both integral-and derivative-time scales. For passive-scalars advected by a turbulent velocity field, the corresponding dynamicmultiscaling exponents are defined via and To obtain bridge relations for dynamic-multiscaling exponents in statistically steady passivescalar turbulence we define the degree-M, order-p integral-time for this case: where Z = 3 + pg − Dû ,θ (h, g) and the argument of the scaling function is suppressed for notational convenience. By using the scaling form of τ pgh (48) and assuming, as in [31], that we get Note that, in contrast to the bridge relations (53) and (54), there is no p-dependence on the right-hand sides of these relations. However, the M-dependence is the signature of nontrivial dynamic multiscaling here. Furthermore, the integral scale bridge relation is meaningful only for those values of M for which ζˆu −M is well defined [7]. We now extend our discussion to the the case of decaying turbulence and define the order-p, time-dependent quasi-Lagrangian velocity structure function (for fluid turbulence) and passive-scalar structure function (for passive-scalar turbulence) as and Fθ p (r, {t 1 , · · · , t p }) ≡ [δθ(x, r, t 1 ) . . . δθ(x, r, t p )] , respectively. For simplicity, we consider the cases t 1 = t 0 + t and t 2 = . . . = t p = t 0 . The derivations of the bridge relations we have given above go through if we assume the following multifractal forms for time-dependent velocity and passive-scalar structure functions, respectively, in decaying turbulence: and i.e. we assume that the multifractal form factorizes into a part Aˆu(r, t 0 ) (or Aˆθ (r, t 0 ) in the case of passive scalars), which depends on the origin of time t 0 , and an integral that is independent of t 0 . This assumption is motivated by the factorization we have seen in section 3 above in (32), (37) and ( In the next section we give compelling numerical evidence in support of our assumptions (64) and (65).

Numerical results
We now present our numerical results for dynamic multiscaling in decaying, homogeneous, isotropic turbulence in models B-D in sections 5.1-5.3, respectively. We show en passant that, for models B and D, equal-time multiscaling exponents are universal in the sense that the exponents obtained from decaying-turbulence runs are equal (within error bars) to the exponents for statistically steady turbulence in these models [6,26,31]. The underlying reason for this universality for passive-scalar turbulence has been discussed earlier [36,37], but we believe that ours is the first numerical demonstration of such universality in these models. Since we use shell models in this section, we replaceû andθ by u and θ , respectively.

Model B
Let us begin with the equal-time exponents ζ θ p for the case of decaying turbulence in model B. As in [26], we use the ESS procedure to obtain exponent ratios from log-log plots ( figure 1(a)) of the equal-time structure functions S θ p (k n ) versus S θ 2 (k n ). In particular, the slopes of the linear regions of the plots in figure 1(a), with 4 n 12, yield the inertial-range exponent ratios ζ θ p /ζ θ 2 . From 50 such statistically independent runs (each run, in turn, is averaged over 5000 independent initial conditions) we calculate the means of the equal-time multiscaling exponent ratios ζ θ p /ζ θ 2 ; these are plotted in figure 1(b) as functions of the order p for 1 p 8. The standard deviations of these exponent ratios (calculated from our 50 runs) provide us with error bars; these are smaller than the symbol sizes used in figure 1(b). Our exponents ratios are in agreement (within the error bars) with earlier results for equal-time exponents for statistically steady passive-scalar turbulence in model B [7,26]. We have also determined the exponent ζ θ 2 directly from log-log plots of θ 2 (k n ) versus k n (equation (22)). We find ζ θ 2 = 1.403 ± 0.003; this agrees well with the analytical prediction [26,38] ζ θ 2 = 2 − ξ , since we use ξ = 0.6 in our simulations. We begin the discussion of our numerical results for time-dependent, passive-scalar structure functions in decaying turbulence by comparing our analytical result (37) with data from our simulations of model B. In figure 2(a), we show the agreement between the two. The error bars have been calculated as in the the case of equal-time multiscaling exponents: the data we present are the means of the values obtained from 50 different statistically independent runs; and the standard deviations of these values yield the error bars. Here and henceforth time is scaled by the large-eddy-turnover time t L . Given the exponential decays of timedependent passive-scalar structure functions in models A and B (equations (32) and (37)), we can extract a unique timescale from plots like the one in figure 2(b). Thus we do not expect time-dependent structure functions to exhibit dynamic multiscaling in models A and B. In particular, equation (37) implies T θ 2 (k n ) = [1/4k 2−ξ n A(ξ )] −1 whence z 2 = 2 − ξ. For our choice of ξ = 0.6 we should, therefore, obtain z 2 = 1.4, and indeed our numerical simulations yield z 2 = 1.398 ± 0.003 ( figure 2(b)). To demonstrate numerically that higher-order time-dependent structure functions also have a dynamic scaling exponent equal to z 2 , as shown in section 3.2 for the fourth-order time-dependent structure function, we show in figures 3(a) and (b) our numerical results for the third-and fourth-orders time-dependent structure functions versus t/t L for 6 n 13 from which we extract T θ 3 (k n ) and T θ 4 (k n ), respectively. The insets show log-log plots of these times versus k n ; the slopes of these plots yield the dynamic-scaling exponents z 3 = 1.398 ± 0.003 and z 4 = 1.402 ± 0.005 in agreement with our expectation z p = 2 − ξ , for all p, since ξ = 0.6 in our calculations.

Model C
The equality of equal-time exponents for decaying and statistically steady fluid turbulence was demonstrated numerically for the Sabra shell model in [23]. We have carried out a  sixth-order scheme. Slopes of log-log plots of T I,u p,1 (n) and T D,u p,2 (n) versus k n (4 n 14) give us z I,u p,1 (figure 5(a)) and z D,u p,2 ( figure 5(b)), respectively. Our results for equal-time and dynamic-multiscaling exponents for the GOY model (for both types of initial conditions) are given in tables 1 and 2, respectively. We compute the multiscaling exponents for equal-time and time-dependent structure functions for 50 different cases. Tables 1 and 2 list the means of these values and their standard deviations yield the error bars. By comparing columns 2 of tables 1 and 2 with column 2 of table II in [6] we confirm, for the GOY model, the weak version of universality [23], i.e. the equal-time exponents ζ u p are the same for both decaying and statistically steady turbulence. Furthermore, our exponents for the GOY shell model agree with those presented in [23] for the Sabra shell model.
A comparison of the remaining columns of tables 1 and 2 with their counterparts in table II of [6] shows that this weak universality also applies to dynamic multiscaling exponents. Moreover, our direct numerical results for z I,u p,1 and z D,u p,2 (columns 4 and 6 in tables 1 and 2) agree with the bridge-relation values of these exponents (columns 3 and 5 in tables 1 and 2) that follow from equations (53) and (54) and ζ u p (columns 2 of tables 1 and 2). Note that the agreement between corresponding entries in tables 1 and 2 shows that our results are insensitive to the type of initial conditions we use. Finally, if we compare these tables with table II in [6], we find that our dynamic-multiscaling exponents for decaying, homogeneous, isotropic turbulence agree with their counterparts for the statistically steady case. Pictorial comparisons of the data in tables 1 and 2 and the results of [6] Figure 6. (a) Plots of z I,u p,1 , with error bars, versus p for model C with data from [6], for statistically steady turbulence, and columns 4 of tables 1 and 2, for decaying turbulence; these plots illustrate the agreement between the three sets of exponents. (b) We compare the derivative-time exponents z D,u p,2 from [6] and columns 6 of tables 1 and 2. As in (a), we find that the three sets of exponents agree, within error bars, with each other.

Model D
From a numerical study of a passive-scalar shell model, with advecting velocities from the Sabra shell model, it was shown [36,37] that the equal-time scaling exponents ζ θ p are universal: they are the same for decaying and statistically steady turbulence; and, in the latter case, they do not depend on the type of forcing. We find, not surprisingly, that this universality holds even when the advecting velocity field is a solution of the GOY model, i.e. for model D: table 3 column 2 shows the values of ζ θ p for 1 p 6, which we have obtained for decaying turbulence; these agree with the results of [7,31] for statistically steady turbulence in model D. Our equal-time exponents are also within error bars of their counterparts for the passive-scalar shell model is equal to µ, with 0 µ 1. We use µ = 0.6, but we have checked in representative cases that our results remain unchanged if we use the range of values 0.3 µ 0.8. Slopes of log-log plots of T I,θ p,1 (n) versus k n yield z I,θ p,1 . To extract the derivative-time scale, we use a centred, finite-difference, sixth-order scheme to obtain T D,θ p,2 (n) from which we get z D,θ p,2 . Integral-and derivative-time multiscaling exponents are extracted from linear fits in the inertial range 4 n 10 as shown in the insets of the representative plots in figures 9(a) and (b) for, respectively, Q θ 5 (n, t) and Q θ 3 (n, t) versus t/t L . Finally, from the bridge relations (61) and our GOY-model, equal-time exponents ζ u −1 = −0.44 ± 0.04 and ζ u 2 = 0.709 ± 0.003, we obtain z I,θ p,1 = 0.56 ± 0.04 and z D,θ p,2 = 0.645 ± 0.003 in agreement with the values from our simulations listed in columns 3 and 4 of table 3. By comparing these columns with their counterparts in table II of [7], we find agreement, within our error bars, between the dynamic-multiscaling exponents for model D for both statistically steady and decaying turbulence.

Conclusions
We have systematized the study of the dynamic multiscaling of time-dependent structure functions in four models (A-D) for passive-scalar (A, B and D) and fluid (C) turbulence. By a suitable normalization of these structure functions, we eliminate their dependence on the origin of time t 0 at which we initiate our measurements. We have shown analytically that for the Kraichnan model of passive-scalar turbulence (models A and B) the two-point time-dependent structure function can be factorized into two parts, one depending on the origin of time t 0 and another that does not. This suggests a suitable normalization of the structure functions by which we can eliminate their dependence on t 0 . Surprisingly the same normalization works for other models of fluid turbulence (model C) and passive-scalar turbulence (model D) as we have shown by extensive numerical simulations. Once this dependence on t 0 has been factored out, the methods developed earlier [6,7] for statistically steady turbulence yield linear bridge relations that connect dynamic-multiscaling and equal-time exponents. We show analytically, for models A and B, and numerically, for models B, C and D, that these exponents and bridge relations are the same for statistically steady and decaying turbulence. Thus we have generalized the universality of equal-time exponents [23] and provided strong evidence for dynamic universality, i.e. dynamic-multiscaling exponents do not depend on whether the turbulence decays or is statistically steady for each of the models A-D. We have shown elsewhere that these exponents do not depend on the dissipation mechanism by studying the GOY model with a hyperviscous term [39]. It is useful to distinguish our results from other studies of temporal dependences of quantities in turbulence. Two such results are described below.
The first example is the result E tot (t) ∼ t −2 , which holds for the total energy E tot at large times t in decaying turbulence [1]. This result holds for times that are longer than the time at which the integral length scale becomes comparable to the size of the system. (In decaying turbulence, the integral scale L int (t) ≡ ( dk E(k, t)/k)/( dk E(k, t)) increases with time since the large-k part of the energy spectrum E(k, t), at time t, gets depleted as t increases.) In all the results we report here, our shell-model analogue of L int is well below the linear size of the system, i.e. L int k −1 0 . Thus our results are not modified significantly by finite-size corrections, nor does the trivial decay of E tot , mentioned above, set in and mask the dynamic multiscaling we have elucidated.
The second example is the intermittency of velocity time increments studied in [40]. These studies calculate structure functions along a single Lagrangian trajectory. Spatial separations, of the sort used in defining equal-time Eulerian structure functions, are replaced by temporal separations along a Lagrangian trajectory. This is distinct from the spatiotemporal structure functions we consider and which are required for a full elucidation of dynamic multiscaling here (and analogous dynamic scaling in critical phenomena). We have described the way in which we have obtained error bars for the equal-time and dynamic-multiscaling exponents. These error bars only account for statistical errors but not the systematic errors associated with the values of n over which we fit inertial-range exponents. One can try to estimate such systematic errors by obtaining local slopes of log-log plots that yield these exponents. However, local slopes can be deceptive in shell models since the values of k n are separated by factors of 2. Instead, we can try to estimate these systematic errors by comparing the exponents we obtain by fitting over the ranges 3 n 15, 4 n 14 and 5 n 13. We have carried out such checks in representative cases for the dynamicmultiscaling exponents we report. The error bars we then obtain are about a factor of 4 larger than those shown in tables 1-3.
We hope our work will stimulate experimental studies of dynamic multiscaling in turbulence. Recent advances in the experimental techniques of particle tracking in turbulent flows [41,42] have made it possible to obtain accurate measurements of Lagrangian properties. To obtain the types of structure functions that we have described it will be necessary at least to track two Lagrangian trajectories of particles that are separated initially by a distance r . At the level of second-order structure functions this has been attempted in the direct numerical simulations of [18]. We now turn our attention to the time-dependent structure function which has the following multifractal form: where G p,h (t/τ p,h ) has a characteristic decay time τ p,h ∼ r/δû r (x) ∼ r 1−h and G p,h (0) = 1. If