On the decay of Burgers turbulence

arXiv:physics/9709002v1 [physics.flu-dyn] 2 Sep 1997 On the decay of Burgers turbulence S.N. Gurbatov1,2 , S.I. Simdyankin1,3 , E. Aurell3,4, U. Frisch2 & G. Tóth5,6 1 Radiophysics Dept., University of Nizhny Novgorod, 23, Gagarin Ave., Nizhny Novgorod 603600, Russia. (Permanent address) 2 Observatoire de la Côte d’Azur, Lab. G.D. Cassini, B.P. 4229, F-06304 Nice Cedex 4, France. 3 Center for Parallel Computers, Royal Institute of Technology, S–100 44 Stockholm, Sweden. 4 Mathematics Dept., Stockholm University, S–106 91 Stockholm, Sweden. 5 Dept. Atomic Physics, Loránd Eötvös University, Budapest 1088, Puskin u. 5–7, Hungary. 6 Sterrekundig Instituut, Postbus 80000, 3508 TA Utrecht, The Netherlands. J. Fluid Mech. 1997, 344, 339–374 (Received 9 December 1996 and in revised form 16 April 1997) 1 Abstract This work is devoted to the decay of random solutions of the unforced Burgers equation in one dimension in the limit of vanishing viscosity. The initial velocity is homogeneous and Gaussian with a spectrum proportional to kn at small wavenumbers k and falling off quickly at large wavenumbers. In physical space, at sufficiently large distances, there is an “outer region”, where the velocity correlation function preserves exactly its initial form (a power law) when n is not an even integer. When 1 < n < 2 the spectrum, at long times, has three scaling regions : first, a |k|n region at very small ks with a time-independent constant, stemming from this outer region, in which the initial conditions are essentially frozen; second, a k2 region at intermediate wavenumbers, related to a self-similarly evolving “inner region” in physical space and, finally, the usual k−2 region, associated to the shocks. The switching from the |k|n to the k2 region occurs around a wavenumber ks (t) ∝ t−1/[2(2−n)] , while the switching from k2 to k−2 occurs around kL (t) ∝ t−1/2 (ignoring logarithmic corrections in both instances). When −1 < n < 1 there is no inner region and the long-time evolution of the spectrum is self-similar. When n = 2, 4, 6, . . . the outer region disappears altogether and the long-time evolution is again self-similar. For other values of n > 2, the outer region gives only subdominant contributions to the small-k spectrum and the leading-order long-time evolution is also self-similar. The key element in the derivation of the results is an extension of the Kida (1979) log-corrected 1/t law for the energy decay when n = 2 to the case of arbitrary integer or non-integer n > 1. A systematic derivation is given in which both the leading term and estimates of higher order corrections can be obtained. It makes use of the Balian– Schaeffer (1989) formula for the distribution in space of correlated particles, which gives the probability of having a given domain free of particles in terms of a cumulant expansion. The particles are, here, the intersections of the initial velocity potential with suitable parabolas. The leading term is the Poisson approximation used in previous work, which ignores correlations. High-resolution numerical simulations are presented which support our findings. 2 1 Introduction The Burgers equation ∂t v + v∂x v = ν∂x2 v, (1) describes a variety of nonlinear wave phenomena arising in the theory of wave propagation, acoustics, plasma physics and so on (see, e.g., Whitham 1974; Gurbatov, Malakhov & Saichev 1991). It was originally introduced by Jan M. Burgers (1939) as a toy model for turbulence. On the one hand, it shares a number of properties with the Navier–Stokes equation : the same type of nonlinearity, of invariance groups and of energy-dissipation relation, the existence of a multidimensional version, etc. On the other hand, it is known to be integrable (Hopf 1950; Cole 1951) and therefore lacks the property of sensitive dependence on the initial conditions (chaos). When this was realized, interest in the Burgers equation as a toy model for Navier–Stokes turbulence abated somewhat. More recently, it became clear that a number of questions connected to the Burgers equation are far from trivial. One class of problems, arising e.g. in surface growth (Barabási & Stanley 1995), leads to a Burgers equation with a random driving force (usually white-noise in time) added to the right hand side (Kardar, Parisi & Zhang 1986; Bouchaud, Mézard & Parisi 1995; Cheklov & Yakhot 1995; Polyakov 1995; Gurarie & Migdal 1996; E et al. 1997). Even in the absence of a driving force, the explicit solution provided by the Hopf–Cole method leads to non-trivial problems in the limit of vanishing viscosity, when random initial conditions are assumed. An example, motivated by work on large-scale structures in the Universe (Zeldovich 1970; Gurbatov & Saichev 1984; Gurbatov, Saichev & Shandarin 1989; Shandarin & Zeldovich 1989), is when the initial velocity field is a Brownian (or fractional Brownian) process in the space variable (She, Aurell & Frisch 1992; Sinai 1992; Vergassola et al. 1994; Aurell et al. 1997). Other examples arise in the sudy of disordered systems (Bouchaud & Mézard 1997). Another instance, which is at the core of the present study, is the longtime behavior of freely decaying solutions in the limit of zero viscosity and, in particular the law of decay of the energy. Before reviewing this issue it is useful to put it into a broader perspective, namely the law of decay for three-dimensional Navier–Stokes turbulence. Kármán & Howarth (1938; see also Lin & Reid 1963), investigating the 3 decay of high-Reynolds number, homogeneous isotropic three-dimensional turbulence, proposed a self-preservation (self-similarity) ansatz for the spatial correlation function of the velocity : the correlation function keeps a fixed functional shape; the integral scale L(t) and the root-mean-square velocity u(t) have a long-time power-law behavior; there are two exponents which can be related by the condition that the energy dissipation per unit mass should be proportional to u3 /L. But an additional relation is needed to actually determine the exponents. Kolmogorov (1941) realized this and proposed using the so-called “Loitsyansky (1939) invariant” as the additional relation and, thereby, derived a law of energy decay u2 (t) ∝ t−10/7 . Proudman & Reid (1954) and Batchelor & Proudman (1956) showed that the integral considered by Loitsyansky is actually not invariant. Under the assumption that the energy spectrum E(k, t) ∝ k 4 at small wavenumbers k, the Loitsyansky integral is, within a numerical factor, just the coefficient of k 4 . For isotropic turbulence, the contribution of the nonlinear term to the time derivative of the energy spectrum, which is called the transfer function, should be proportional to k 4 at small k, thereby preventing the constancy of the Loitsyansky integral. However, if initially E(k) ∝ |k|n with a “spectral exponent” n < 4, then the coefficient of |k|n will be invariant. This is called here the principle of “permanence of large eddies” (PLE). Kolmogorov’s argument is easily adapted to that case and gives then u2 (t) ∝ t−2(n+1)/(n+3) (see, e.g., Frisch 1995, Section 7.7). Almost exactly the same arguments may be applied to the one-dimensional Burgers turbulence. In one dimension the transfer function should be proportional to k 2 at small k (see Section 3). Hence, if the initial energy spectrum E0 (k) ∝ |k|n with −1 < n < 2, Kolmogorov’s argument predicts again u2 (t) ∝ t−2(n+1)/(n+3) . We shall see that this is correct only when −1 < n < 1. For Navier–Stokes turbulence, a rigorous derivation of the law of energy decay seems ruled out at present.1 The situation is of course more favorable for Burgers-type turbulence and it will be our main goal in this paper to critically examine the issues of self-similarity and energy decay at long times. In this paper we shall exclusively consider the behavior of the solutions in the limit of vanishing viscosity for Gaussian initial data. For the effect of a finite 1 It is not even clear, in three dimensions, if for the initial value problem there is any energy decay at all in the limit of vanishing viscosity, since the issue of breakdown of smoothness for the Euler flow is moot (see, e.g., Frisch 1995, Section 7.8). 4 (small) viscosity and/or non-Gaussian initial data (e.g., piecewise constant or shot noise velocities), the reader is referred to recent papers by Funaki, Surgailis & Woyczynski (1995), Surgailis (1996, 1997) and Hu & Woyczynski (1996). The paper is organized as follows. In Section 2, we formulate our problem and list some elementary results about the Burgers equation to be used in what follows. Section 3 is devoted to phenomenological derivations of the main results about the decay of random solutions to the Burgers equation, the nature of which depends crucially on the spectral exponent n. The most delicate results, obtained when n > 1, are rederived by systematic methods in Sections 5 and 6. In Section 4 we discuss previous work on this issue. In Section 5 we derive the law of decay of energy as the leading term in a systematic expansion. Our basic tool is the cumulant expansion (Balian & Schaeffer, 1989) of the probability of finding regions in space without particles. In Section 6 we calculate the spatial correlation function, and prove the permanence of large eddies. In Section 7 we summarize our main findings and derive consequences, such as the absence of globally self-similar decay for the energy spectrum when 1 < n < 2. Section 8 presents highresolution numerical experiments which confirm such findings. Section 9 presents concluding remarks. The Appendix contains a rederivation of the Balian–Schaeffer formula in a form appropriate for our problem. 2 Formulation and elementary results We shall be concerned with the initial-value problem for the unforced Burgers equation in one dimension : ∂t v + v∂x v = ν∂x2 v, v0 (x) ≡ v(x, 0). (2) Use will also be made of the conservation form of the Burgers equation : ∂t v = ∂x v2 − + ν∂x v 2 ! (3) and of the formulation in terms of the (velocity) potential : 1 ∂t ψ = (∂x ψ)2 + ν∂xx ψ, 2 5 (4) where v = −∂x ψ. (5) As is well-known, the solution to the Burgers equation with positive viscosity ν, has an explicit integral representation obtained by Hopf (1950) and Cole (1951). We are mainly interested here in the solutions in the limit of vanishing viscosity. Use of Laplace’s method then leads to the following “maximum representation” for the potential (Hopf 1950) in the limit of vanishing viscosity (henceforth always understood) : ψ(x, t) = max Φ(x, a, t) a " Φ(x, a, t) ≡ ψ0 (a) − (6) (x − a)2 , 2t # (7) where ψ0 (a) is the initial potential. Expansion of the square shows that (6) is basically a Legendre transformation. From (5) and (6) it follows that v(x, t) = x − a(x, t) , t (8) where a(x, t) is the coordinate at which Φ(x, a, t) achieves its (global) maximum for given x and t. It is easily checked that a is the Lagrangian coordinate from which emanates the fluid particle which will be at x at time t (see, e.g., Gurbatov et al. 1991). In what follows we shall work with homogeneous random velocity fields. We list the main definitions. The correlation function is Bv (x, t) ≡ hv(x, t)v(0, t)i, (9) where the angular brackets denote ensemble averages. Its Fourier transform, the energy spectrum is 1 E(k, t) ≡ 2π Z ∞ Bv (x, t) exp(ikx) dx. (10) −∞ The initial values of the correlation function and of the energy spectrum are denoted B0v (x, t) and E0 (k), respectively. The (mean) energy is E(t) ≡ hv 2 (x, t)i = Z ∞ −∞ 6 E(k, t) dk. (11) The initial energy is denoted σv2 ≡ hv02 (x)i = Z ∞ E0 (k) dk. (12) −∞ We now introduce the basic assumptions of this paper. The initial velocity field is homogeneous and Gaussian with an energy spectrum of the form E0 (k) = α2 |k|n b0 (k), (13) where n is the spectral exponent. The even and non-negative function b0 (k) is assumed to satisfy b0 (0) = 1. It is smooth and decreases faster than any inverse power of k at infinity (this implies, of course, that the initial velocity correlation function is very smooth). b0 (k) can be characterized by a wavenumber k0 around which lies most of the initial energy and which is, in order of magnitude, the inverse of the initial integral scale L0 . The usual definition of the integral scale, which involves the integral of the correlation function, is not appropriate since this integral vanishes for any n > 0. Instead, we use σψ L0 ≡ , (14) σv where σψ2 is the variance of the initial potential E0 (k) dk. (15) k2 −∞ This definition of L0 will be convenient in deriving various long-time laws in Sections 5 and 6. As we shall see, the value of the spectral exponent n determines the universality class of the decaying solution at long times. In order for the velocity to be homogeneous its variance must be finite; hence, we must have n > −1. When −3 < n < −1 the velocity has infinite variance and only increments are homogeneous; the velocity itself behaves, at large distances, as a Brownian or fractional Brownian motion process, a case considered by She et al. (1992), Sinai (1992), Vergassola et al. (1994), which is not within the scope of the present paper. When −1 < n < 1 the velocity is homogeneous whereas the potential behaves, at large distances, as a Brownian or fractional Brownian motion process. This case has already been investigated by Burgers (1974), by She et al. (1992) and by Avellaneda & E (1995) and will be considered here only briefly. We shall concentrate on the case n > 1, when both the velocity and the potential are homogeneous. σψ2 ≡ hψ02 (x)i = 7 Z ∞ 2.1 Elementary results We observe that the smoothness of the energy spectrum at the origin depends on the spectral exponent n. When n is an even integer, our assumption implies that the initial correlation function B0v (x), being the nth derivative of a smooth function of x decreasing rapidly at infinity, also has these properties. When n is odd or non-integer the spectrum is not smooth. It follows that the correlation function decreases as a power law (Pitman 1968; Bleistein & Handelsman 1975), viz 2 B0v (x) ≃ α2 Cn |x|−n−1 , |x| → ∞, (16) where πn . (17) 2 Next, we observe that, when n > 1, the initial potential ψ0 (x) is a homogeneous random function while, for n < 1, it only has homogeneous increments. In the former case, the potential has a translation-invariant correlation function, denoted by B0ψ (x), the Fourier transform of which is E0 (k)/k 2 . In the latter case, it is simpler to characterize the initial potential by its structure function S0ψ (x) ≡ h(ψ0 (x) − ψ0 (0))2 i, (18) Cn = −2Γ(n + 1) sin which is still translation-invariant and exists as long as n > −1. It is related to the energy spectrum by S0ψ (x) = 4α2 Z 0 ∞ (1 − cos(kx)) k n−2 b0 (k) dk. (19) When n > 1 the structure function goes to a finite limit at large separations while, when n < 1, it grows without bound : S0ψ (x) ∼ ( σψ2 n > 1 α2 |x|1−n n < 1, |x| → ∞. (20) An important consequence of (4), (5) and (11), valid whenever the velocity is homogeneous, irrespective of the potential being or not being homogeneous, is E(t) = hv 2i = h(∂x ψ)2 i = 2∂t hψ(x, t)i. (21) 2 In this paper the symbol ≃ is used for the leading-order term in asymptotic expansions, ∼ is used to relate two quantities which may differ by an order unity dimensionless constant and ∝ when the quantities are just proportional. 8 It relates the energy to the time-derivative of the mean potential, which does not dependent on the position x. We now recall some results about the Burgers equations which are relatively standard and will be used in what follows. When n ≥ 0, the correlation invariant J(t) = Z ∞ hv(x, t)v(0, t)i dx = E(0, t) (22) −∞ is time independent. Indeed, using the conservation form of the Burgers equation (3), it is checked that the time derivative of the integral in (22), limited from −X to +X is expressed in terms of correlations between quantities evaluated at points a distance X apart. Such correlations vanish as X → ∞ because, with our assumptions, the initial conditions have decreasing correlations at large separation and this property will be preserved.3 Note that J is just the value of the energy spectrum at k = 0. A straightforward consequence of the Burgers equation with a homogeneous velocity field is 1 ∂t Bv (x, t) = ∂x S3 (x, t) + 2ν∂x2 Bv (x, t), 6 (23) which relates the correlation function Bv (x, t) and the third order structure function S3 (x) ≡ h[v(x, t) − v(x, 0)]3 i. This equation is the analogue of the Kármán–Howarth equation for Navier–Stokes turbulence. Just as the latter implies Kolmogorov’s four-fifths law (see, e.g., Frisch 1995), eq. (23) implies a “twelve law” : in the limit of vanishing viscosity and at scales much smaller than the integral scale, we have S3 (x) = −12εx, (24) where ε ≡ ∂t hv 2 /2i is the mean energy dissipation. Eq. (23) or (24), in the limit of zero viscosity, puts a strong constraint on self-similar decaying solutions both for Navier–Stokes and Burgers turbulence. By self-similar decaying solutions, we understand the following : when the velocity v and the space coordinate x are rescaled by suitable 3 This is expected because influence cannot propagate faster than the initial velocities. 2 2 So, correlations over a distance X should not exceed a value ∝ e−(1/2)X /(tσv ) ) . It should be possible to prove this rigorously for the kind of initial conditions assumed here. 9 time-dependent scaling factors u(t) and L(t), the (single-time) statistical properties of the solution become independent of the time. This definition implies in particular the following scalings for the correlation function and the third-order structure function of the velocity : 2 Bv (x, t) = u (t)B̃v S3 (x, t) = u3 (t)S̃3 ! x L(t) ! x , L(t) (25) (26) with dimensionless time-independent functions B̃v and S̃3 . Substituting (25)(26) into (23) for ν → 0 and demanding that B̃v and S̃3 be time independent, we obtain two first-order differential equations relating u(t) and L(t), namely, L̇ u̇ = α1 , L̇ = α2 u, (27) u L where the dot denotes time differentiation and α1 and α2 are two dimensionless undetermined constants. Note that exactly the same equations are obtained with a self-similarity ansatz for Navier–Stokes turbulence. Eq. (27) has power-law solutions4 but, as already pointed out by Kolmogorov (1941), the exponents in the power laws cannot be determined solely from the selfsimilarity ansatz. Additional ingredients are needed. For what follows it will be convenient to give a precise definition of the scaling factors. We choose u(t) to be the r.m.s. velocity u(t) ≡ hv 2 (x, t)i1/2 , (28) and we define the “integral scale” L(t) by taking α2 = 1 in (27), that is L̇(t) = u(t). (29) Let us finally note that in Fourier space the self-similarity ansatz (25), together with (29), translates into E(k, t) = L3 (t) Ẽ(kL(t)), t2 (30) 4 In what follows, we shall be concerned with decaying solutions which satisfy the selfsimilarity conditions (25)-(26) and (27) only in an asymptotic sense, for long times and to leading order. We shall then see that u(t) and L(t) may take the form of power laws with logarithmic corrections. 10 Ẽ(k̃) = 1 2π Z ∞ B̃v (x̃) exp(ikx̃) dx̃, (31) −∞ where our choice of normalization of the energy as u2 (t) imposes that Ẽ(k̃) dk̃ = 1. R 3 3.1 Phenomenology of the decay Permanence of large eddies and self-similar decay An important role in the decay at long times is played by the principle of permanence of large eddies (PLE). We shall give two formulations almost – but not quite – equivalent, one in physical space and the other one in Fourier space. PLE : physical space formulation. We assume that the Gaussian initial velocity has a correlation function with a power-law behavior at large distances : Bv (x) ≃ α2 Cn |x|−n−1, |x| → ∞, (32) with an exponent n > −1. Then, (32) remains true at any later time, with the same exponent n and the same constant α. PLE : Fourier space formulation. We assume that the Gaussian initial velocity has a spectrum with a power-law behavior at small wavenumbers : E(k) ≃ α2 |k|n , |k| → 0, (33) with an exponent 2 > n > −1. Then, (33) remains true at any later time, with the same exponent and constant. In its physical-space formulation PLE can be supported by the following argument. From (32), it follows that the characteristic velocity associated with a large separation x is u0 (x) ∼ α|x|− n+1 2 , (34) so that the corresponding “turnover time”, the characteristic time for significant changes through nonlinearity, is t0 (x) ∼ n+3 x ∝ |x| 2 , u0 (x) 11 (35) which grows without bound at large separations. Hence, at any finite time, eddies of sufficiently large size should remain essentially as they were initially. This simple argument could be upset if two fluid particles, initially strongly correlated (say, separated by not more than one initial integral scale), were to propagate this correlation to a large separation x in a finite time.5 However, with Gaussian initial conditions this occurrence has a very small probability, which can be bounded by ∝ exp[−(1/2)x2 /(σv t)2 ], as explained in the footnote following the derivation of the constancy of J (eq. (22)) in the previous section. Thus, the initial correlation function at large separations should remain unaffected to leading order. In its Fourier space formulation, PLE has already been proven in the previous section, for the case of the spectral index n = 0. When n < 2, PLE is a consequence of the observation that the time-rate of change of the energy spectrum, the so-called transfer function, is typically proportional to k 2 at small wavenumbers, because of the beating interactions of pairs of Fourier components with nearly opposite wavevectors. We may also derive the constancy of the coefficient of |k|n just from the constancy of the coefficient of |x|−n−1 at large x. There are reasons why the two formulations are not quite equivalent. First, the n = 0 result in Fourier space is not an instance of PLE in physical space. Indeed, when n is an even integer, (32) implies a non-analytic spectrum at small k with logarithmic corrections. Second, when n > 2 and the initial correlation function has the behavior (32) at large separations, the energy spectrum at any later time will begin with an analytic (time-dependent) k 2 term with a subleading non-analytic correction proportional to |k|n (with a time-independent coefficient). We now show how PLE, used as a “boundary condition” at large x or small k in connection with the self-similarity ansatz (25) and (29), allows us to derive the law of temporal variation of the integral scale and the energy. We first give the argument in physical space. Let us assume that the distance beyond which PLE applies is within the scope of the self-similarity formula (25). It then immediately follows from (25), (29) and (32) that B̃(x̃) ∝ 5 For the case of incompressible Navier–Stokes turbulence, pressure effects can instantaneously affect the large-scale correlations and thereby invalidate PLE, at least in the anisotropic case (Batchelor & Proudman 1956). 12 |x̃|−n−1 and 2 (36) L(t) ∼ (αt) 3+n and thus, by (29), 4 E(t) = u2 (t) ∼ α 3+n t− 2(n+1) 3+n . (37) In Fourier space the self-similarity ansatz (30), together with (33) gives Ẽ(k̃) ∝ |k̃|n , at small k and the same relations for integral scale and the energy as written above. Clearly, this argument cannot be applied with initial data such that the spectral index n ≥ 2, since the later spectrum has now a k 2 dependence at small k with a time-dependent coefficient. Actually, we shall now show that the laws (36) and (37) apply only when n < 1. We begin with the case −1 < n < 1 when the initial potential has homogeneous increments. Many aspects of this case are well understood, thanks in particular to Burgers’ own work (1974; see also Gurbatov et al. 1991). The phenomenology is quite simple. Increments ∆ψ0 (L) of the initial potential over a distance L = x − a can be estimated from the square root of the structure function S0ψ (L) of the potential (18). For a given position x, the maximum in (6)-(7) will come from those as such that the change in potential is comparable to the change in the parabolic term q S0ψ (L) ∼ L2 . t (38) At large times L is also large and we can use (20) to obtain precisely (36) and, hence, the energy law (37). Alternatively, we could argue that when t is so large that the parabolas appearing in (7) have a radius of curvature much larger than the typical radius of curvature of features in the initial potential, we can plausibly replace that initial potential by fractional Brownian motion of exponent h = (1 − n)/2 (in order to be consistent with (20)). Without loss of generality we may assume that this fractional Brownian motion starts at the origin for x = 0. This function is then statistically invariant under the transformation x → λx and ψ0 → λh ψ0 . It is then elementary, using (6)-(7), to prove that a rescaling of the time is (statistically) equivalent to a suitable rescaling of x distances and of ψ(x, t). This implies (36) and (37).6 6 For details, see Vergassola et al. (1994, Section 4) and Avellaneda E (1995). 13 3.2 Elementary derivation of the Kida law of decay The case n > 1 when the initial potential is homogeneous is more difficult. Because of this homogeneity, the structure function becomes constant at large distances. If we then formally use (38) we obtain a t−1 law for the energy and a t1/2 law for the integral scale. This is almost, but not quite, the right result. There are logarithmic corrections, first discovered by Kida (1979), which can actually be captured by phenomenology, as now explained. We shall estimate the energy in terms of the mean potential at x = 0, using (21) and the maximum representation (6)-(7). We assume that, at long times t, the distribution of values of the potential is peaked around its (non-vanishing) mean value hψi ≡ hψ(0, t)i. This quantity is thus the typical maximum value of ψ0 (a) − a2 /(2t). Equivalently, we ask : what is the value hψi by which we should vertically shift the parabola a2 /(2t) so that the graph of ψ0 (a) has a substantial probability of being everywhere below the shifted parabola ? To estimate this probability we divide the a-axis into equal intervals, I1 , I2 , I3 , . . . , with a common length of about one initial integral scale L0 . The initial potential takes approximately independent values in the different intervals. Thus, the probability of no intersection with the shifted parabola is (approximately) the product of the probabilities q1 , q2 , . . . of no intersections in the various intervals : Pno inters. ∼ q1 q2 . . . . (39) Each of the qi s is of the form 1 − pi , where pi is the probability of at least one intersection in the interval Ii . This can be estimated as a2 pi ∼ W hψi + i . 2t ! (40) Here, W (ψ) is the single-point cumulative probability of ψ0 , integrated from ψ to infinity and ai is any point in the interval Ii , since the parabolic function changes very slowly with a when t is large. All the pi s are small and thus Pno inters. ∼ 1 − p1 − p2 . . . . (41) The sum of the pi s can be approximated by an integral from −∞ to +∞. When hψi is taken very large this integral is very small. As we decrease 14 hψi, the probability of having an intersection somewhere becomes significant when hψi is chosen such that a2 W hψi + 2t −∞ Z ∞ ! da = O(1). L0 (42) Since we assumed Gaussian initial conditions with r.m.s. value σψ ,the cumulative probability W (ψ) is approximately exp [−(1/2)(ψ/σψ )2 ]. The leading order of the integral in (42) is obtained by keeping only the first two terms in the expansion of the square of hψi + a2 /(2t). Thereby, we obtain hψi ∼ σψ ln1/2 where tnl ≡  t , tnl  (43) L20 σψ (44) is the characteristic nonlinear time. Using (21), we then obtain the Kida (1979) log-corrected 1/t law for the energy decay : E(t) ∼ t−1 σψ ln−1/2  t . tnl  (45) Finally, assuming (29), we obtain Kida’s (1979) result for the integral scale : L(t) ∼ t1/2 (σψ )1/2 ln−1/4  t . tnl  (46) We must, of course, stress that this derivation is only heuristic. We shall come back to the issue of a more systematic derivation (Sections 5 and 6), including estimates of the subleading corrections. It is straightforward to extend this phenomenology to initial potentials which are homogeneous and non-Gaussian, as the only ingredient used is the p.d.f. of ψ0 . For example, when the p.d.f. has a tail ∝ exp(−|ψ|β ), the logarithmic correction to the 1/t law has the exponent (1 − β)/β. Hence, when n > 1 the law of decay could be very sensitive to the functional form of the p.d.f. of the initial potential. Such heuristic considerations for nonGaussian potentials have not yet been tested by more systematic theory and simulations. 15 Figure 1: Sketch of the energy spectrum at long times when 1 < n < 2. 3.3 The Kida law and the issue of self-similarity when 1<n<2 One of the most questionable aspects of our phenomenological derivation is its validity when n is not an even integer. Indeed, this gives rise to longrange algebraically decreasing correlations, which are neglected in this phenomenology and which Kida (1979) explicitly ruled out. Actually, we shall see in Section 5 that the Kida law remains valid nevertheless as long as the initial potential is homogeneous, that is for any n > 1. We now address for the first time one of the central issues of this paper. The predictions (36) and (37) for E(t) and L(t), obtained by assuming selfsimilarity and using PLE as a “boundary condition” at large x (or small k) are not consistent with the Kida law (45)-(46). Because of the time-dependence of the coefficient of k 2 in the spectrum at small k, the applicability of the former is obviously restricted to n < 2. Since the latter applies when n > 1, there may be a contradiction in the region of overlap 1 < n < 2. Actually, we shall see that n > 1 leads to a self-similar decay of the solution in an “inner” region of not too large distances (or, in Fourier space, of not too small 16 wavenumbers), controlled by the integral scale L(t) (46). In this inner region correlations of ψ decrease like a Gaussian function and the Fourier signature is a range of wavenumbers at which the energy spectrum is ∝ k 2 . In addition, if n > 1 and is not an even integer, there exists an “outer” frozen region at which PLE applies. Thus, there is, at very large distances, an algebraic tail in the correlation function and the spectrum has a time-independent |k|n contribution at small k. When 1 < n < 2 this is the leading term at small k so that the complete spectrum has three ranges : proportional to |k|n at very small k, then proportional to k 2 at intermediate wavenumbers (up to approximately the inverse of the integral scale L(t)) and, finally, proportional to k −2 at high k (the usual signature of shocks), as shown in Fig. 1. We shall return to the issue of the global picture in Section 7. 4 Homogeneous initial potential : previous work In this section we discuss previous work dealing with the limit of vanishing viscosity for initial potentials which are homogeneous and Gaussian. The main papers under discussion are Kida (1979), Gurbatov & Saichev (1981; see also Gurbatov et al. 1991), Fournier & Frisch (1983) and Molchanov, Surgailis & Woyczynski (1995). For initial conditions such that the velocity field is homogeneous but the potential is not, or such that the initial conditions are not Gaussian, the reader is referred to Burgers (1974), Albeverio, Molchanov & Surgailis et al. (1994), and Esipov & Newman (1993), Esipov (1994), Hu & Woyczynski (1996) and references therein. Let us first state the essence of what is at presently known to hold with certainty. In the limit of vanishing viscosity, as the time t tends to infinity, the statistical solution becomes self-similar in the sense of Section 2. The integral scale L(t) and the energy E(t) are given, to leading order, by t , 2πtnl   t −1/2 −1 E(t) ≃ t σψ ln , 2πtnl L(t) ≃ (tσψ )1/2 ln−1/4 17   (47) (48) where tnl ≡ L0 L20 = , σψ σv L0 ≡ σψ . σv (49) The non-dimensionalized self-similar correlation function B̃v (x̃), defined by (25), which is a function of x̃ = x/L(t), is given by d (x̃P (x̃)) , dx̃ B̃v (x̃) = (50) where P (x̃) = and 1 2 Z dz ∞ −∞ g  x̃+z 2  exp g(z) ≡ h Z (x̃+z)2 8 z i +g  x̃−z 2  exp h (x̃−z)2 8 s2 e− 2 ds. i (51) (52) −∞ Note that the properties of the self-similar state are universal in so far as they are expressed solely in terms of two integral characteristics of the initial spectrum, namely the initial r.m.s. potential σψ and r.m.s. velocity σv .7 It may be shown that the function P (x̃) is the probability of having no shock within an Eulerian interval of length x̃L(t) (Gurbatov, Malakhov & Saichev 1991). Important consequences of (50)-(51) are the asymptotic behavior of the non-dimensionalized correlation function B̃v (x̃) at large x̃ and of its Fourier transform, the non-dimensionalized spectrum Ẽ(k̃) (31) at small k̃ = kL(t) (Gurbatov, Malakhov & Saichev 1991) : π x̃2 B̃v (x̃) ≃ − , x̃ exp − 32 8 ! r Ẽ(k̃) ≃ µ2 k̃ 2 , |k̃| → 0, Z ∞ 1 µ2 = x̃2 P (x̃) dx̃ ≈ 1.08. π 0 |x̃| → ∞, (53) (54) (55) Most of these results have already been obtained by Kida (1979). He was interested in the long-time behavior when the initial potential is homogeneous 7 Observe that the spectral exponent does not directly enter, in contrast to what happens when n < 1 (cf. (36)-(37)). 18 and the initial correlations decrease rapidly at large separations, without particularly emphasizing the Gaussian assumption. He introduced a model of discrete independent maxima whose p.d.f. had three free parameters (an exponent and two constants); their relation to the properties of the initial conditions (say, the spectrum) were left unspecified. For the case of a p.d.f. with a Gaussian tail, he obtained the functional form of the results given above. At long times, the contacts of the parabolas, in the maximum representation (6)-(7) of the solution, are near large-value local maxima of the initial potential. Kida’s assumption of independence of the successive maxima implicitly assumes a Poisson process. Gurbatov & Saichev (1981) who conjectured the asymptotic existence of a Poisson process8 and Molchanov, Surgailis & Woyczynski (1995) who proved it, showed that, in the x-ψ plane, the density of the points is uniform in the x-direction and exponential in the ψ-direction. This permits the calculation of the one- and two-point p.d.f.’s of the velocity (Gurbatov & Saichev 1981). Molchanov, Surgailis & Woyczynski (1995) also calculated the full N-point multiple time distributions. Their proof made crucial use of some fine properties of extremes of Gaussian processes (Leadbetter, Lindgren & Rootzen 1983).9 From the point of view of the present paper it is important to stress that Molchanov, Surgailis & Woyczynski’s (1995) result is valid for any n > 1 and not just for n ≥ 2. 10 Assuming a Poisson distribution, Gurbatov & Saichev (1981) and Fournier & Frisch (1983) showed that the statistical properties of the points of contact between the parabolas and the initial potential can be obtained from the statistical properties of their intersections, whose mean number can be calculated using the formula of Rice (1954; see also Leadbetter, Lindgren & Rootzen 1983; Papoulis 1991). Thus, they could express the parameters in the asymptotic formulas in terms of the r.m.s. initial potential and velocity. None of the above theories give a control over how fast in time the asymptotic results are reached. Furthermore, the rescaling involved in these theories can only capture what we have called in the previous section the “inner re8 See also Fournier & Frisch (1983). What they proved is not exactly in the form of the long-time asymptotic relations given above but an equivalent result in which the time and space variables and the amplitude of the initial conditions are rescaled. 10 They assume that the initial correlations of ψ decrease faster than the inverse of ln |x| at large |x| and can be Taylor-expanded to fourth order at small |x|. 9 19 gion”, where PLE is irrelevant. We shall show in the next section that, for finite long times, the Poisson distribution is just an approximation which can be systematically refined, so as to address such issues. 5 Systematic derivation of the Kida law for the energy decay using the Balian–Schaeffer formula The (mean) energy is given by (21) in terms of the time derivative of the mean potential hψi. Here, we use the method of Fournier & Frisch (1983) and the Balian–Schaeffer formula (A.18), established in Appendix A, to derive the asymptotic behavior of the mean potential for the case of a homogeneous initial potential, that is, n > 1. Let P (A, t) denote the probability density of ψ at time t and let Q(A, t) denote the (cumulative) probability to have ψ < A, given by Q(A, t) = ZA ′ ′ P (A , t) dA . (56) −∞ The mean value of the potential at time t is then expressible as hψi = Z∞ AP (A, t) dA. (57) −∞ An integration by parts gives hψi = Z∞ [1 − Q(A, t)] dA − 0 Z0 Q(A, t) dA. (58) −∞ We shall see that the second integral in (58) may be bounded by a timeindependent constant, while the first grows; so it is only the first integral that will be of importance in what follows. From (6)-(7) and the homogeneity of the initial potential, we have y2 Q(A, t) = Prob ψ0 (y) − − A < 0, 2t " 20 # for all y . (59) The argument of the probability in (59) can be written as the condition of non-intersection of the following two graphs, shown in Fig. 2 : G1 : y 7→ ψ0 (y), y2 G2 : y 7→ A + . 2t (60) We can immediately use this geometric picture to eliminate the second integral in (58) from further discussion. Indeed, the probability that G1 and G2 never intersect is obviously less than the probability that ψ0 (y) is less than A at the single point y equal to zero. The integral over A from minus infinity to zero of this quantity is finite and bounded by a quantity of the order of σψ . This follows from the observation that, for negative A, the probability Q(A, t) given by (59) is less than the same expression evaluated at the single point y = 0. Hence, the contribution of the second integral in (58) is bounded by a time-independent constant and will not contribute to the logarithmically growing leading-order mean potential. We now treat the intersections of the graphs G1 and G2 as “particles” and apply the Balian–Schaeffer formula (A.18) of Appendix A to calculate the probability of having zero particles (no intersections at all). We shall limit ourselves to the case where G1 intersects G2 going upward (see the footnote on p. 30 for what happens when all intersections are counted.). Such intersections are depicted by circles in Fig. 2. From (59), (60) and (A.18), we have (−)n Q(A, t) = p0 = exp cn , n=1 n! ∞ X ! (61) where the cumulants cn are constructed from the n-particle distribution functions Fn as explained in Appendix A. R In particular, c1 = hNi = F1 (y)dy, the mean number of particles on the whole real line. We now claim that we have the Rice (1954) formula : hNi = Z F1 (y)dy = Z ′  δ (Φ (y, t)) {Φy }+ dy , 2 (62) where Φ(y, t) = ψ0 (y) − y2t − A and the shorthand {x}+ denotes the positive part of x, i.e. equal to x if x is greater than zero and zero otherwise and δ is the delta-function. 21 Figure 2: Intersections of the initial potential ψ0 (y) (curve G1 ) with the parabola G2 . Indeed, the integral in (62) only gets a contribution when the function Φ vanishes. This happens precisely when the two graphs G1 and G2 intersect. Let y ∗ be a point of intersection. An integral over a small interval around the intersection then contributes to (62) a quantity |Φ′y (y ∗, t)|−1 {Φ′y (y ∗, t)}+ . If the derivative of Φ at y is positive, that is, if G1 intersects G2 going upward, then the contribution is one. If, on the other hand, G1 intersects G2 going downward, then the contribution is zero, because {Φ′y (y ∗ , t)}+ vanishes. The integral in (62) therefore counts the number of upward intersections (upcrossings) between the two graphs. Note that the assumed smoothness of the process ψ0 (y) (which follows from the fast decrease at infinity of its spectrum) renders the applicability of the Rice formula robust. R The F1 (y)dy is thus determined by the joint probability density  integral  ′ ′ WΦ,Φ′ Φ, Φ ; y of Φ and its derivative Φ at the point y : hNi = Z F1 (y)dy = Z dy Z ′ ′  ′  dΦ {Φ }+ WΦ,Φ′ Φ, Φ ; y . (63) We have similar expressions for the functions Fn (y1 , . . . , yn ). Let us introduce the n-point joint probability density Wψ,v (ψ1 , . . . , ψn ; v1 , . . . , vn ) of 22 the initial potential ψ and the initial velocity v, the latter being equal to minus the space derivative of the initial potential. Then we can write, analogously to (63) : y2 y2 Fn (y1 , . . . , yn ) = dv1 . . . dvn Wψ,v 1 + A, . . . , n + A; v1 , . . . , vn · 2t 2t ·{v1 − y1 /t}+ {v2 − y2 /t}+ · · · {vn − yn /t}+ . (64) Since the initial potential is Gaussian, all the joint probability densities are completely determined by the initial two-point correlation function : ! Z B0ψ (x) = hψ0 (y + x) ψ0 (y)i . (65) In our case we need the joint probability densities for (initial) fields ψ at various points and their (negative) space derivatives v at those same points. They are given by 1 Wψ,v (ψ1 , . . . , ψn ; v1 , . . . , vn ) = h i1/2 · (2π)2n det M   1 T −1 · exp − (ψ1 , . . . , ψn , v1 , . . . , vn ) M (ψ1 , . . . , ψn , v1 , . . . , vn ) , 2 (66) where M is the enlarged correlation matrix Mψℓ ,ψℓ′ = B0ψ (yℓ − yℓ′ ) ′ Mψℓ ,vℓ′ = −Mvℓ′ ,ψℓ = B0ψ (yℓ − yℓ′ ) (67) ′′ Mvℓ ,vℓ′ = −B0ψ (yℓ − yℓ′ ) = B0v (yℓ − yℓ′ ) . In the last line we have introduced the initial velocity correlation B0v . The diagonal elements of this matrix are Mψℓ ,ψℓ = σψ2 , Mvℓ ,vℓ = σv2 . (68) Here, σψ2 and σv2 are the variances of the initial potential and the initial velocity, respectively, as previously defined. We now return to the Balian–Schaeffer expression (61) for the cumulative probability of the potential. Inspection of (58) shows that the mean potential 23 comes predominantly from those arguments A such that Q(A, t) is not too close to unity or, by (61), such that the alternating sum of cumulants is not much smaller than unity. As A grows, the changeover of the alternating sum to a value of order unity happens around a value which we shall denote A∗ (t). 5.1 The Poisson approximation We are interested in the energy decay at large times and it will be shown that in this case the first cumulant c1 = f1 is the leading term in (61) and that, to leading order, we have the Poisson approximation : Q (A, t) ≃ exp {−c1 } . (69) It is so called for the natural reason that if the upcrossings obeyed a Poisson distribution, then (69) would be true exactly. The rest of this subsection contains a computation of the first cumulant c1 . The function F1 is determined by the single-point probability density of the initial potential ψ and the velocity v, the subscript 0 being dropped for conciseness. For a homogeneous random function, these two variables are uncorrelated. We then obtain from (63) and (64) F1 (y1 ) =  1 2πσv2 σψ2  Z 2 v2 1 (y12 /2t + A) + exp − 2 σψ2 σv2 " #! · {v − y1 + } dv. (70) t It is now useful to introduce the auxiliary function Z (x) =  1 2π 1/2 Z v2 {v − x}+ e− 2 dv, (71) which can obviously be expressed in terms of the error function. We only need to use that, for x → 0, the function tends to Z(0) = (2π)−1/2 and that, for x → ∞, we have Z(x) ≃ x. In terms of Z(x) the mean number of upcrossings is given by (  2) 2 (y /2t + A) y1 q 1 1/2πσψ2 exp − c1 = F1 (y1 ) dy1 = dy1 σv Z . tσv 2σψ2 (72) The cumulant c1 is a dimensionless number. It is therefore preferable to rewrite (72) in a manifestly dimensionless form. By (14), σψ /σv is equal to Z Z  24 L0 , the initial integral scale. In (72) we introduce the new dimensionless variables z1 = y1 /L0 and a = A/σψ and the dimensionless nonlinear time tnl = L0 /σv already defined (49). Eq. (72) can then be written c1 = Z z1 t/tnl dz1 Z    !q 1/2π exp −    z12 2(t/tnl ) 2 2   +a  . (73)   It is seen that (73) involves two spatial scales : when z1 is small compared to the “prefactor scale” t/tnl we can simply substitute Z(0)qfor Z(z1 /(t/tnl )); when z1 is large compared to the “Gaussian scale” l(a, t) ≃ (t/tnl )/a there is a large damping by the term in the exponent quadratic in z1 . If we assume, for the moment, that the typical value of a∗ q = A∗ /σψ where c1 becomes order unity does not grow in time faster than t/tnl ,11 it follows that it is self-consistent to assume that the prefactor scale is much larger than the Gaussian scale, so that we may safely use Z(0) in (73). The remaining integral in (73) cannot be carried out explicitly, but it is nevertheless easy to evaluate it to relevant leading order. We have two terms in the exponent that depend on z1 : one quadratic and one quartic. The quartic term is larger for z1 much beyond the Gaussian scale, but then the damping is strong anyway; so, the correction from the quartic term is small, of relative order a−2 and thus negligible. It is then easily shown that, for those large values of a ∼ a∗ which give the dominant contribution to the mean potential at large-time, we can use the following leading-order approximation for the cumulant c1 : t c1 (A, t) ≃ tnl 2πa  1/2 2 /2 e−a , A = σψ a. (74) We now determine a∗ (t). For this, we assume that c1 is the dominant term in the Balian–Schaeffer formula (61) for the probability of no intersections, as will be shown in Section 5.2. a∗ (t) is then determined by the condition that c1 (σψ a∗ , t) be order unity. To avoid undetermined constants, we define a∗ (t) by the condition c1 (σψ a∗ , t) = 1, that is through the transcendental equation  11 t tnl 2πa∗ 1/2 ∗ )2 /2 e−(a = 1, (75) Actually, it will be seen to grow only as a square root of the logarithm of the time. 25 from which it is easily shown that, to leading order, at large times 1/2 ∗ a ≃ ln  t . tnl 2π  (76) We now return to the cumulative probability of Q(A, t). Using the Poisson approximation (69) and the leading-order approximation (74) for c1 , we obtain the following doubly exponential distribution :12 Q(A, t) ≃ exp{− exp{−ζ}}, (77) where ζ is obtained from a by shifting and rescaling : ∗ a=a ! ζ 1+ ∗ 2 . (a ) (78) Since a∗ (t) is logarithmically large at large times, it follows from (77)-(78) that the values of a = A/σψ contributing to the mean potential (58) are = ∆a ≃ (a∗1)2 ≃ concentrated near a∗ (t) with a small relative width ∆A A∗ a∗ ln 2πtt nl . Hence, to leading order, the mean potential is given in the Poisson approximation by : hψ(t)i ≃ σψ a∗ ≃ σψ ln1/2  t , 2πtnl  (79) which implies the energy decay law (48). Notice that, within the Poisson approximation, the energy decay is fully determined by the local properties of the initial potential correlation function ′′ B0ψ (x), namely, σψ2 = B0ψ (0) and σv2 = −B0ψ (0) and does not involve its large-scale behavior. 5.2 Corrections to the Poisson approximation So far we have only reproduced a known result. The Balian–Schaeffer formula allows us to find the applicability condition for the Poisson approximation by estimating contributions of higher-order cumulants to (61). 12 This doubly exponential distribution is equivalent to the distribution of maxima being Poisson with uniform density in space and exponentially decaying in a (Molchanov, Surgailis & Woyczynski 1995). 26 We thus turn to the calculation of the second cumulant c2 appearing in (61), given in terms of the one- and two-particle densities by (A.17), repeated hereafter for convenience : c2 = Z Z [F2 (y1 , y2 ) − F1 (y1 )F1 (y2 )] dy1 dy2 . (80) The function F2 is now given by (64). Our purpose in this section will only be to show that the contribution of c2 – and of the higher-order cn s – are subdominant compared to that of c1 . For this purpose, it will be enough to estimate the order of magnitude of c2 (for values of a ∼ a∗ (t)). Thus, there is no need to derive the full leading-order expression for c2 with all the numerical constants. From (64)-(67) it follows that, unlike c1 , the second cumulant c2 involves the entire initial correlation function B0ψ (x). Let us denote by ∆corr a typical correlation length for the initial potential.13 We shall now distinguish those contributions to (80) where the two arguments y1 and y2 are separated by not more than ∆corr and those where the separation is larger or, possibly, much larger, a case of interest when the initial correlation function has a slow algebraic decrease. The former will be called “local” (loc.) and the latter “non-local” (nloc.). Using (64) and (80) we can estimate cloc. ∼ 2 Z Z F2 (y1 , y2 ) dy1dy2 ∼ ∆corr Z (F1 (y1 ))2 dy1 . (81) Introducing as in (73) new (dimensionless) variables z1 and a, we get cloc. ∼ 2 ∆corr L0 Z dz1 Z z1 t/tnl !2   z12 !2   1 exp − +a  2π 2(t/tnl ) . (82)  The integral in (82) is practically the same as for the first cumulant, so we have   ∆corr t 1/2 −a2 loc. c2 (A, t) ∼ e , A = σψ a. (83) L0 tnl a The main differences with the expression (74) for c1 is (i) the presence of a factor ∆corr /L0 and (ii) the fact that the argument of the exponential is now 13 We could have taken ∆corr = L0 but keeping the distinction has some interest, as we shall see. 27 −a2 instead of −a2 /2. From (75) and (83) we find that, around the value a∗ (t) which give the dominant contribution to the mean potential, the second order local cumulant is roughly cloc. ∼ 2  t ∆corr · L0 tnl a∗ (t)  !−1/2 . (84) It easily follows that the local second order cumulant may be neglected compared to c1 as soon as   ∆corr 2 t ≫ tnl . (85) L0 Here, we make a short digression. There are instances where (∆corr /L0 )2 can be large. Consider for instance a quasimonochromatic signal with a center wavenumber k0 ∼ L−1 and a width ∆k ∼ [∆corr ]−1 ≪ k0 . In this 0 case, the condition t/tnl ≫ 1 is not enough for the Poisson approximation (69) to hold. The asymptotic regime described by Kida’s law is then obtained much later. As long as tnl ≪ t ≪ tnl (∆corr /L0 )2 , the energy decays much faster : E(t) ∼ L20 /t2 (Gurbatov & Malakhov 1977). The physical reason for this is a strong correlation of the shocks in the early stage of the evolution, which prevents the rapid merging of shocks.14 Let us now turn to the non-local contribution to c2 which is important only when n is not an even integer, so that there is an algebraic tail in the correlation function : B0ψ (x) ∼ σψ2  ∆corr x n−1 , x ≫ ∆corr . (86) Owing to this, when 1 < n < 2 the main contribution to c2 comes from the non-local region where |y1 − y2 | ≫ L0 while |y1 |, |y2 | ≤ L(t). Here, L(t) is the integral scale at time t, given by (46). The main technical difficulty now is that we have to express the inverse matrix of correlations that appears in the exponent in the Gaussian joint probability densities (64)-(67). Actually, there is no need to fully invert the matrix. It is sufficient to note that when the points are far apart, correlations become weak and the off-diagonal elements will be small compared to the 14 Similar remarks can be made when the initial energy spectrum has a very long plateau at intermediate wavenumbers. 28 diagonal ones. We may therefore invert perturbatively to first order in the off-diagonal elements. Furthermore, the velocity-velocity correlations and the velocity-potential correlations decay much faster with separation than the potential-potential correlations. Therefore, it is consistent to keep in the inverted matrix only off-diagonal elements proportional to B0ψ (yℓ − yℓ′ ). As 1 before we can substitute for the function Z(x) its value at the origin, (2π)− 2 . After all these simplifications we can write the following expression for the non-local part of the second cumulant c2 : 1 cnloc. = 2 (2πL0 )2    1 dy1 dy2 · exp − 2   2σ NL ψ " ! ( ( 2 y1 B0ψ (y1 − y2 ) +A · exp σψ4 2 Z !2 !2    ·  y12 y2 +A + 2 +A 2 2 !#) ) y22 +A −1 . 2 (87) Here, the subscript NL (non-local) is a reminder that we only integrate over points y1 and y2 sufficiently far apart. The next step is to expand the second exponential, which contains the correlation function, to linear order. As in the evaluation of c1 , we change the variables y1 , y2 and A into the dimensionless variables z1 , z2 and a, respectively and obtain, using (86) : cnloc. ∼ 2  · ∆corr L0 Z NL n−1 · 1 · (2π)2   !2 z12 +a 2(t/tnl ) 1 dz1 dz2 · exp −  2 ·|z1 − z2 | 1−n z12 +a 2(t/tnl ) ! !2   z12 +a + 2(t/tnl ) z22 +a , 2(t/tnl ) ! ·  (88) the integration being now over the region such that |z1 − z2 | > ∆corr /L0 . The integral (88) is estimated in the following way : the exponentials effectively delimit the region of integration to a circle in the (z1 , z2 )-plane with radius l(t) ∼ (t/(atnl ))1/2 . Inside this circle, most of the factors in (88) are of order unity. Therefore, we can estimate the non-local part of c2 as cnloc. ∼ 2  ∆corr L0 n−1 1 · (2π)2 Z 2 NL dz1 dz2 · e−a a2 |z1 − z2 |1−n , 29 (89) where NL is now the region |z1 | <  t tnl a 1 2 , |z2 | <  t tnl a 1 2 , |z1 − z2 | > ∆corr . L0 (90) When n > 2 the non-local component of c2 is much smaller than the local one and we may ignore it. When 1 < n < 2 the integral in (89) is readily evaluated, and we obtain the estimate cnloc. ∼ 2  ∆corr L0 n−1 −a2 2 ·e a  t tnl a 1 2 ·  t tnl a  2−n 2 . (91) Finally, estimating cnloc. for A ≃ A∗ (t), we obtain 2 cnloc. ∼ 2  ∆corr L0 n−1  tnl t (n−1)/2 . (92) Comparison with (84) indicates that the non-local part of c2 is much larger than the local one, but still small compared to unity, so that we may indeed neglect c2 compared to c1 .15 Let us observe that when n is just slightly larger than one, an exceedingly long time may be needed before the asymptotic Kida law applies since c2 is then just marginally smaller than c1 . Using similar procedures we can show that cloc. and cnloc. for m > 2 m m are even smaller than the c2 -corrections. This establishes the validity of the Kida law (48) at long times. 6 The correlation function at long times and the preservation of long-range correlations In this section we calculate the (spatial) correlation function of the potential at long times for the solution to the Burgers equation, using a Balian– Schaeffer formula, adapted now to the evaluation of two-point quantities. 15 This statement holds provided only upcrossings or only downcrossings between the graphs G1 and G2 are counted. Otherwise, two things happen : (i) the mean number of crossings and hence c1 is doubled; (ii) strong correlations are present between successive up and down large-ordinate crossings; these invalidate the Poisson approximation; actually, the coefficient c2 becomes large and exactly cancels the aforementioned factor two. 30 As we shall see, in evaluating the far tail of the correlation function, the dominant contribution will come from the second rather than from the first cumulant. For reasons of symmetry we take the following definition for the correlation function of the potential : Bψ (x, t) = hψ(x/2, t)ψ(−x/2, t)i − hψ(x, t)i2 , (93) from which the the velocity correlation function is obtained as Bv (x, t) = − ∂ 2 Bψ (x, t) . ∂x2 (94) To evaluate (93) we need, in addition to the one point probability density P1 (A, t), the joint probability density P2 (A, B; x, t) of the potential at the same time t for two points with spatial separation x. From (93), we obtain Bψ = Z Z AB(P2 (A, B; x, t) − P1 (A, t)P1 (B, t)) dAdB. (95) We also introduce the cumulative probability Q2 (A, B; x, t), related to P2 (A, B; x, t) by ∂2 Q2 (A, B; x, t). ∂A∂B P2 (A, B; x, t) = (96) Using this in (95) and integrating by parts twice, we have Bψ (x, t) = Z Z ∞ [Q2 (A, B; x, t) − Q1 (A, t)Q1 (B, t)] dAdB. (97) −∞ By the maximum representation solution (6) we can express Q2 in the form Q2 (A, B; x, t) = Prob (G1 never intersects G2 ) , (98) where G1 and G2 are now the following graphs, shown in Fig. 3 : G1 : y 7→ ψ0 (y),  x  y− 2 G2 : y → 7 g(y) = min  2t 31 2 + A,  y+ 2t x 2 2  + B . (99) Figure 3: Intersections of the initial potential ψ0 with two parabolas. Thus, the procedure of calculation of Q2 is similar to that used in the preceding section for calculating the mean potential, but with a more complex function g(y). Hence, we may again use the Balian–Schaeffer formula and obtain ! ∞ X (−)n cn , (100) Q2 (A, B; x, t) = exp n=1 n! where the cn s depend now on A, B, x and t. The first cumulant c1 (A, B) in the expansion (100) is now c1 (A, B) = Z ( 1 1 g(y1) Z (g (y)/σv ) exp − 1/2 (2πσψ ) 2 2σψ2 ′ ) dy1 . (101) From (99) and Fig. 3 we see that the integral in (101) can be divided into two integrals over pieces of the two parabolas intersecting at y ∗ : 1 c1 (A, B) = 2πL0         Z  y∗   exp −  −∞        32 " (y+ x2 ) 2t 2 #2     +B    2σψ2       dy + y∗ = Z −∞ y∗        exp − A−B t. 2x "       (y+ x2 ) 2t 2 #2     +A    2σψ2    , dy          (102) (103) Here, we have assumed from the start that t ≫ tnl , to replace Z(x) by Z(0) = (2π)−1/2 . It may now be checked that (102) behaves in different ways depending on how the widths of the Gaussians compare to the separation x. As long as x is smaller than or comparable to the integral scale L(t) given by (47), the leading-order term in (100) is found to be c1 , that is, the Poisson approximation remains valid. One then derives the expressions (50)-(51) for the correlation function previously obtained by Kida (1979), Gurbatov & Saichev (1981) and Molchanov, Surgailis & Woyczynski (1995). Actually, to derive the self-similarly evolving correlation function given by (50), the only possible spatial scaling factor is (up to a multiplicative constant) the integral scale L(t) given by (47). It may also be checked that L(t) is (up to a constant factor) equal to hψ 2 (t)i1/2 /hv 2(t)i1/2 When |x| ≫ L(t) the situation is quite different. Now, the two parabolas intersect at a very large ordinate so that, to leading order, we have two independent integrations over two different almost complete parabolas. Thus, to leading order, we have c1 (A, B; x, t) ≃ c1 (A, t) + c1 (B, t), x ≫ L(t). (104) This, however, implies Q2 (A, B; x, t) − Q1 (A, t)Q1 (B, t) ≃ exp{−c1 (A, t) − c1 (B, t)} − exp{−c1 (A, t)} exp{−c1 (B, t)} = 0, (105) which, in turn, implies the vanishing of Bψ (x, t). It follows that the leadingorder contribution cannot come from the Poisson approximation. This leading-order contribution actually comes from the second order cumulant Z Z c2 = [F2 (y1 , y2) − F1 (y1 )F1 (y2 )] dy1dy2 , (106) 33 Figure 4: The different regions of integration appearing in (106). Near the diagonal y1 = y2 we have regions I and I’ of (107). Away from the diagonal we have regions II and II’. where F1 and F2 are now the one- and two-particle densities constructed from the crossings of the G1 and G2 graphs defined in (99). It may be checked that when |x| ≫ L(t) the significant regions in (106) are : x x , y2 ∼ ), 2 2 x x II : (y1 ∼ , y2 ∼ − ), 2 2 I : (y1 ∼ x x I’ : (y1 ∼ − , y2 ∼ − ), 2 2 x x II’ : (y1 ∼ − , y2 ∼ ), 2 2 (107) also illustrated in Fig. 4. In the regions I and I’ we have the same estimates (84)-(92) for the value of c2 as in the calculation of the mean potential. The corresponding contribution 34 to the correlation function goes to zero with separation even faster than the c1 -contribution and is not relevant. In the regions II and II’, for t ≫ tnl and |x| ≫ L(t) we obtain after integration over v1 , v2 , the following expression 1 1 c2 = 2 (2πL0 )2     Z Z   dy1 dy2 exp − −∞ · exp B0ψ (x)     ∞  y12 2t  +A σψ4   y22 2t  y12 2t 2σψ2   +B  ) ( 2 +A   ! −  y22 2t  − 1 ABB0ψ (x) −1 . ≃ c1 (A)c1 (B) exp σψ4 +B 2σψ2 2    ·    (108) Just as in Section 5.1, the joint probability density of A, B is localized in a narrow region near A ≃ B ≃ hψi ≃ σψ (ln t/2πtnl )1/2 . We introduce in (108) the new variables ζ and ξ through A = hψi + σψ2 σ2 ζ, B = hψi + ψ ξ , hψi hψi (109) and we observe that, when ζ and ξ are O(1), we have, by (74) and (79), c1 (A) ≃ e−ζ , c1 (B) ≃ e−ξ . (110) We now evaluate the correlation of the potential (97) for |x| ≫ L(t). For this we use (100), truncated to its first two terms (108) and (110), and the fact that B0ψ (x) is small, so that we can expand exponentials with B0ψ (x) in their arguments. We finally obtain : Bψ (x, t) ≃ B0ψ (x), t ≫ tnl , |x| ≫ L(t). (111) This establishes the preservation of long-range correlations, an instance of the permanence of large eddies (PLE), as defined in Section 3. 7 The global picture at long times We are now in a position to derive the main results of this paper, some of which have already been presented in a phenomenological way in Section 3. 35 Our derivation of the Kida law (48) for the energy decay in Section 5, using the Balian–Schaeffer formula, proves the validity of this law for any n > 1 and not just n > 2 as thought previously. Similarly, for any n > 1 the integral scale L(t) is given by (47). More novel results emerge when considering the global structure, at long times, for the correlation function (in physical space) and for its Fourier transform, the energy spectrum. We shall consider only the situation at scales larger than the integral scale L(t). Indeed, we have already stated in Section 6 that, at scales small or comparable to the integral scale, we recover known results (see Kida 1979; Gurbatov, Malakhov & Saichev 1991). We observe that for any n > 1 which is not an even integer the correlation function Bv (x) has two different asymptotic regions. In an “outer asymptotic region” at extremely large x, permanence of large eddies (in the form (111) established in Section 6) implies that Bv (x, t) ≃ Bv (x, 0) ≃ Cn α2 |x|−n−1 , |x| → ∞. (112) There is also an “inner asymptotic region” where the correlation function has the (approximately) Gaussian decrease given by (53) in dimensionless form with x̃ = x/L(t). Using (25) we can rewrite this in dimensioned form as Bv x L(t) ! x x2 π 2 ≃− , u (t) exp − 2 32 L(t) 8L (t) ! r (113) which can be reexpressed in terms solely of L(t), using u(t) = L̇(t). The separation ls (t) at which the correlation function switches from its inner to its outer form, called here the switching distance, is obtained by equating (112) and (113), yielding16  ls (t) ∼ L(t) (n − 1) ln  t tnl 1/2 . (114) Hence, the switching distance is only logarithmically larger than the integral scale.17 Still, the fact that L(t) and ls (t) have (slightly) different scaling laws 16 To truly characterize the transition from the inner to outer asymptotic expansion it is necessary to perform matched asymptotic expansions which require more than the leading order. Equating the leading-order inner and outer terms correctly predicts the order of magnitude of the scale at which this transition takes place. 17 When n is an even integer the algebraic outer region disappears altogether. 36 in time implies that, globally, the correlation function does not evolve in a self-similar fashion. This phenomenon becomes much more prominent if we consider the energy spectrum, the Fourier transform of the correlation function. In dimensionless inner-region variables, the energy spectrum has the following form (Gurbatov, Malakhov & Saichev 1991) : Ẽ(k̃) ≃    0.36k̃ −2, k̃ ≫ 1 1.08k̃ 2 , k̃ ≪ 1 , k̃ ≡ kL(t). (115) The k̃ −2 region is the signature of shocks and will not interest us further, while the k̃ 2 region comes about because the (inner) correlation function of the potential has a non-vanishing integral from −∞ to +∞. In dimensioned variables the small-k inner-region behavior of the spectrum is thus E(k, t) ≃ L5 (t) 2 k ≃ A(t) k 2 , t2 where kL(t) ≪ 1, 1 inner region, (116) t t 2 −5/4 . (117) ln A(t) = tnl 2πtnl Instead of a Gaussian region (in physical space), we have a spectrum with an algebraic k 2 region and a time-dependent coefficient A(t). We must now distinguish two cases. When n > 2, the k 2 contribution (116) dominates everywhere over the |k|n contribution, the latter being then only a subdominant term in the small-k expansion.18 The more interesting case is 1 < n < 2, already sketched in Fig. 1. The permanence of large eddies implies now that, at extremely small k, σv2 L30  E(k, t) ≃ α2 |k|n ,  for k → 0,  outer region. (118) This relation holds only in an outer region |k| ≪ ks (t) where (118) dominates over (116). The switching wavenumber ks (t), obtained by equating (116) and (118), is given by ks (t) ≃ α2 t2 L5 (t) ! 1 2−n ≃ L−1 0  t tnl − 1 2(2−n) ln 5 4(2−n)  t . 2πtnl  (119) 18 This subdominant term is easily seen in numerical simulations at very short times and moderately small wavenumbers, for example when n = 3 (A. Noullez, private communication). 37 Some comments are now made. First, it may seem paradoxical that the switching wavenumber ks (t) is not the inverse of the switching distance ls (t). The reason is that in physical space we have a Gaussian competing with a power law, whereas in Fourier space we have two power laws competing; this renders inapplicable naı̈ve Fourier phenomenology in which distances are roughly inverse wavenumbers. Second, let us define an energy wavenumber kL (t) = L−1 (t), which is roughly the wavenumber around which most of the kinetic energy resides. From (47) kL (t) ∼ (tσψ )−1/2 (ignoring logarithmic corrections). We then have from (119), still ignoring logarithmic corrections : n−1 ks (t) t − 2(2−n) ∼ . (120) kL (t) tnl Hence, the switching wavenumber goes to zero much faster than the energy wavenumber : in Fourier space the separation of the two regions is much more manifest. This is why the results of the numerical simulations in Section 8 will be presented in Fourier space. Let us also observe that the ratio of the energy in the outer region to the total energy, a measure of how well the Kida law is satisfied, is equal to 3(n−1) (t/tnl )− n−2 (up to logarithms) and thus becomes very small when t ≫ tnl , unless n is very close to unity. Thus, we have established the central claim of this paper, namely that for 1 < n < 2 there is no globally self-similar evolution of the energy spectrum. Of course, as n → 2 the inner k 2 region overwhelms the outer |k|n region and as n → 1 the converse happens, so that in both instances global self-similarity tends to be reestablished.  8  Numerical experiments The purpose of this section is to obtain numerical confirmation of the rather unusual long-time regime predicted by the theory when the spectral exponent n is between one and two. There should then be three different power-law regions in the energy spectrum : E(k) ∝ |k|−2 for |k| ≫ L−1 (t); E(k) ∝ |k|2 for ks (t) ≪ |k| ≪ L−1 (t); E(k) ∝ |k|n for |k| ≪ ks (t), 38 (121) where ks (t) and L(t) are given by (119) and (47), respectively. Clearly, such simulations require extremely high spatial resolution since three power-law ranges are present. Fortunately, for the Burgers equation there is no need to numerically integrate (1) : advantage can be taken of the maximum representation (7) to directly construct, from the initial data, the solution in the limit of vanishing viscosity at any given time t. The principles of the method have been presented in detail elsewhere (She, Aurell & Frisch 1992; Noullez & Vergassola 1994; Vergassola et al. 1994; Aurell, Gurbatov & Simdyankin 1996) and will only be briefly recalled here. The method exploits two observations. First, the maximum representation (7) becomes a Legendre transformation after expansion of the quadratic term. Second, when the problem is discretized on a grid of N points, the search for the maximizing Lagrangian point a(x, t) can be done by a dichotomic procedure which uses the non-decreasing property of the inverse Lagrangian map. The number of operations can thus be reduced from a brute-force estimate O(N 2 ) to O(N ln2 N). In the simulations it is of course necessary to introduce both large-scale and small-scale cutoffs. The small-scale cutoff must be chosen sufficiently small compared to the integral scale to leave room for a k −2 inertial range. Since the integral scale grows with the time, this constraint is not very stringent. In practice, we take the initial integral scale about a factor ten larger than the mesh. For the large-scale cutoff, we assume periodic boundary conditions which facilitate Fourier transformations and make no significant difference as long as we analyse scales much smaller than the spatial periodicity. In our simulations, the maximum wavenumber is taken to be 1/2, so that wavevectors run from ±1/N to ±1/2, where N is the resolution, that is the number of grid points in physical space. In all our simulations we take N = 221 ≈ 2 × 106 . The energy spectrum for the initial Gaussian velocity is E0 (k) = αn2 |k|n − e k2 2k2 0 , k = 0, ± 2 1 1 ,± ,...,± , N N 2 (122) with k0 ≈ 0.07 and αn chosen such that hψ02 (x)i = 1. The initial integral scale and nonlinear time, as defined in (14) and (44), respectively, are then 39 given by L0 =  1/2 n−1 1  2   n+1 k0 Γ 2  Γ  ,   n−1 1 Γ 2  . tnl = 2 k0 Γ n+1 (123) 2 With the chosen value of k0 , the integral scale is around ten and the nonlinear time around one hundred. The initial potential is generated by fast Fourier transform (FFT) from its Fourier coefficients ψ̂0 (k) which are complex Gaussian variables of variance h|ψ̂0 (k)|2 i = E0 (k)/(Nk 2 ). After the solutions at suitable output times have been constructed with the procedure described above, they are Fourier transformed and their energy spectra are calculated by averaging typically over 103 realizations.19 The precise value of the spectral exponent 1 < n < 2 in our runs depends, of course, on the particular tradeoffs desired : the closer n is to 1, the more conspicuous is the persistence of the initial |k|n range; and the closer it is to 2, the more conspicuous is the (nonlinearly generated) k 2 range. A natural unit for the output times is t∗ = N 2 . Indeed, we know from (47) that the integral scale grows roughly as t1/2 . Thus t∗ is roughly the time needed for the integral scale to grow from an initial value O(1) to a value comparable to the largest scale available in the simulation, which is O(N). In practice, we take output times ti = 10i−6 N 2 , i = 1, 2, . . . , 8. (124) We describe now the results of our numerical experiments. Fig. 5 shows the evolution of the spectrum for n = 1.7 and is mostly meant to bring out qualitative aspects. At large wavenumbers a k −2 inertial range develops very quickly.20 At short times the spectrum preserves its initial value at all but the largest wavenumbers, in accordance with the permanence of large eddies. The wavenumber where the spectrum achieves its maximum, which we use as an operational measure of the inverse integral scale at time t,21 moves to smaller and smaller wavenumbers. At large times the spectrum appears to change in a self-similar way. 19 This Monte Carlo procedure is easily parallelized in the implementation on an IBM SP2 machine. 20 The time scale for shocks to develop is tnl . 21 It is actually almost exactly equal to this (see Gurbatov, Malakhov & Saichev 1991, Fig. 5.13). 40 −8 10 (1) (2) (3) energy spectrum −10 10 (4) (5) −12 (6) 10 (7) (8) −14 10 −16 10 −7 10 −6 10 −5 10 −4 −3 10 10 wavenumber −2 10 −1 10 0 10 Figure 5: Evolution of the energy spectrum with an initial spectrum (dashed line) proportional to k n (n = 1.7) at small wavenumbers k. Resolution N = 221 . Spectra averaged over 103 realizations. The labels correspond to output times ti = 10i−6 N 2 , (i = 1, 2, . . . 8). We turn to a more quantitative description. Fig. 6 shows the energy spectrum at the third output time t3 = 10−3 N 2 . We also plot the initial spectrum with its |k|n tail which should be preserved at small wavenumbers (leading term of outer expansion) and the leading term of the inner expansion, obtained by Fourier transformation of (50), which has a k 2 tail at small wavenumbers. The intersection of the latter two is the switching wavenumber ks (t) as defined by (119). It is seen that, except in the immediate neighborhood of the switching wavenumber, the leading terms of the inner and outer expansions are in very good agreement with our simulations. Fig. 7, communicated by A. Noullez, gives additional evidence for the permanence of large eddies. 41 −7 10 −8 10 −9 10 −10 energy spectrum 10 −11 10 −12 10 −13 10 −14 10 −15 10 k = ks −16 10 −7 10 −6 10 −5 10 −4 −3 10 10 wavenumber −2 10 −1 10 0 10 Figure 6: Energy spectrum at t = t3 = 10−3 N 2 . Same conditions as Fig. 5. Initial spectrum : dashed line. Leading-order asymptotic theory in inner region : dashed dotted line. The vertical line shows the switching wavenumber. We examine now the temporal evolution of various quantities for which we have theoretical predictions. Fig. 8 shows the evolution of the integral scale compared to the theoretical prediction(47). The agreement is excellent provided that the logarithmic correction is taken into account. By (28) and (29) this also implies the validity of Kida’s log-corrected law for the energy decay. Fig. 9 shows the corresponding information for the coefficient A(t). Fig. 10 shows the evolution of the switching wavenumber ks (t), compared with the theory (119). Finally, Fig. 11 shows the evolution of the switching wavenumber, but now for a spectral index n = 1.5. The agreement with the leading-order theory is slightly better than for n = 1.7. 42 6 -12 t = 10 -10 t = 10 -8 t = 10 -6 t = 10 -4 t = 10 compensated spectrum E(k)⋅k -1.5 5 4 3 2 1 0 0 10 10 2 10 4 10 6 wavenumber Figure 7: Evolution of the energy spectrum compensated by a k −n factor in order to reveal the region of “permanence of large eddies” (PLE) as a plateau with unit height. Spectral exponent n = 1.5. Resolution N = 220 . Spectra averaged over 30,000 realizations. The unit of time is different from the other figures, the spatial period being here unity. Observe the progressive shrinking of the PLE zone. (Courtesy A. Noullez.) 9 Concluding remarks The main conclusions have already been presented in Section 7 (see also the Abstract). Here, we make additional comments. We have given an explicit formula only for the leading-order term of the law of energy decay. The higher-order terms have just been bounded : they appear to be smaller (in relative terms) by negative powers of the time t. It would be useful to have explicit expressions for the subdominant terms. Our derivation of the law of energy decay uses an asymptotic expression for the probability of non-crossing of the initial potential and suitable parabo43 5 10 4 integral scale 10 3 10 2 10 1 10 2 10 4 10 6 10 time 8 10 10 10 Figure 8: Evolution of the computed integral scale L(t) (asterisks) compared to the theoretical leading-order prediction (solid line) and the same without the logarithmic correction (dashed line). Same conditions as Fig. 5. las. Kida’s (1979) original derivation obtained this from the distribution of shocks, which he related to the distribution of (high local) maxima in the initial potential. Actually, it should be possible to reformulate our theory in terms of maxima rather than of upcrossings. Roughly, we expect that high-level upcrossings should occur very close to high-level maxima. This statement cannot be exactly true since two successive up- and downcrossings of the (initial) potential with a parabola need not be separated by a maximum of the potential. However, crossings of a parabola having a high-lying apex are very likely to happen in the flat region near the apex of the parabola where up- and downcrossings are separated by a maximum. Quantifying such remarks would be of special interest for extension of the theory to the multidimensional Burgers equation. Indeed, in d > 1 space dimensions, the equivalent of the crossings are (d−1)-dimensional manifolds for which the extension of the method presented here is far from straightforward. But (local) maxima are still point processes and can be handled via the d-dimensional version of the Balian–Schaeffer formula (Balian & Schaeffer 1989). 44 −1 10 −2 10 −3 A(t) 10 −4 10 −5 10 −6 10 2 10 4 10 6 10 time 8 10 10 10 Figure 9: Evolution of the coefficient A(t) in front of the k 2 term at intermediate wavenumbers (asterisks), compared to the leading-order theoretical prediction (117) (solid line) and the same without its logarithmic correction (dashed line). Same conditions as Fig. 5. There are a number of other interesting extensions of the present work, such as the effect of viscosity and of non-Gaussian initial potentials. Instances of these have already been studied by Funaki, Surgailis & Woyczynski (1995), Surgailis (1996), Hu & Woyczynski (1996) and Surgailis (1997). Finally, a completely open question is the extension of the present ideas to three-dimensional incompressible turbulence. Does the mechanism which prevents globally self-similar decay for 1 < n < 2 have a counterpart for the Navier–Stokes problem ? Acknowledgements. We have benefited from discussions with A. Noullez. W.A. Woyczynski and S. Lehr have kindly sent us various papers indicated by a referee. This work was supported by RFBR-INTAS through grant 95IN-RU-0723 (S.N.G., S.I.S., U.F. and E.A.), by RFBR through grant 96-0219303 (S.N.G. and S.I.S.), by the Swedish Institute (S.I.S.), by the Swedish Natural Science Research Council through grant S-FO-1778-302 (E.A.), by a 45 2 10 0 10 −2 switching wavenumber 10 −4 10 −6 10 −8 10 −10 10 −12 10 2 10 4 10 6 10 time 8 10 10 10 Figure 10: Evolution of the switching wavenumber ks (t), compared to the leading-order theoretical prediction (119) (solid line) and the same without its logarithmic correction (dashed line). Same conditions as Fig. 5. grant from the Fondation des Treilles (E.A., U.F., S.N.G. and G.T.), by the French Ministry of Higher Education (S.N.G. and G.T.) and by the Hungarian Science Foundation grant OTKA F-4491 (G.T.). S. Gurbatov and G. Tóth thank the Observatoire de la Côte d’Azur, S. Simdyankin and E. Aurell the Center for Parallel Computers, and E. Aurell the Radiophysical Department of Nizhny Novgorod University, for hospitality and creating nice and pleasant working conditions. 46 0 10 −1 10 −2 switching wavenumber 10 −3 10 −4 10 −5 10 −6 10 −7 10 −8 10 2 10 4 6 10 10 8 10 10 10 time Figure 11: Evolution of the switching wavenumber ks (t) for n = 1.5. Otherwise, as in Fig. 10. Appendix A The Balian–Schaeffer formula for “particle” distributions The original work of Balian & Schaeffer (1989) was concerned with the following type of problem. One has a distribution of “particles” (galaxies, in the applications considered by them) in the three-dimensional space which is homogeneous, i.e. translation-invariant and one wants to determine the probability that there is no particle at all within some finite domain. For the application to the Burgers equation, the role of the particles will be played by the upcrossings between the initial potential and a given parabola, which constitute a point process in the sense of Leadbetter, Lindgren & Rootzen (1983). Although the potential is a homogeneous random function, the presence of the parabola makes our problem inhomogeneous. Furthermore, our 47 domain is the whole real line and is thus not bounded.22 Nevertheless, our approach will remain rather close in spirit to the original Balian–Schaeffer approach. In the remainder of this appendix we shall use the word “particle” rather than “intersection”, since the formalism has obviously applications beyond the Burgers equation. Let there be given a random set of particles on the real line. The total number of particles is assumed to be random but almost surely finite. We are actually interested in finding the probability of having no particle at all. We shall characterize this random set by the so-called n-particle densities. F1 (y1) dy1 is the probability of having a particle within the interval [y1 , y1 + dy1 ]. Similarly, Fn (y1 , y2 , . . . , yn ) dy1dy2 . . . dyn is the probability of having (at least) n different particles within, respectively, the intervals [y1 , y1 + dy1 ], [y2 , y2 + dy2], . . . and, finally, [yn , yn + dyn ]. Note that when the particles are independent (an assumption not made here), Fn is just a product of F1 R factors. Note also that the Fn s are not normalized. For example F1 (y1 )dy1 is the mean number of particles.23 Balian and Schaeffer’s work leads to expressions for pn , the probability of having exactly n particles in the domain, in terms of the Fn s. We shall here concentrate on the evaluation of p0 , the probability of having no particle at all. We begin by introducing the exclusive n-particle densities corresponding to the case where there are exactly n particles. Pn (y1 , y2 , . . . , yn ) dy1dy2 . . . dyn is the probability of having a particle within the interval [y1 , y1 +dy1 ], another particle within the interval [y2 , y2 + dy2], . . . , a particle within the interval [yn , yn + dyn ] and no other particle. Observe that Z Pn (y1 , y2 , . . . , yn ) dy1dy2 . . . dyn = n! pn . (A.1) The factor n! is present because the Pn s are here defined for unlabeled particles. We shall use (A.1) to calculate pn for n = 1, 2, .. and then obtain p0 from the obvious relation p0 = 1 − p1 − p2 − . . . pn − . . . , (A.2) 22 The formalism developed hereafter applies identically to bounded and unbounded domains. The integral sign will always denote a simple or multiple integration over the whole domain. 23 This follows from the observation that the mean number of particles in the interval [y1 , y1 + dy1 ] differs by O(dy12 ) from the probability of having one particle within the interval [y1 , y1 + dy1 ]. 48 expressing that the total probability is unity. We show now that the pn s may be expressed in terms of the Fn s. First we relate the Fn s to the Pn s. To have a particle in [y1 , y1 + dy1] we must have one of the following mutually exclusive situations : (i) a total of exactly one particle in [y1 , y1 + dy1 ], (ii) a total of exactly two particles, one of them in [y1 , y1 + dy1 ] and the other one at any place, . . . . Hence, F1 (y1 ) = P1 (y1 ) + + 1 2! Z Z P2 (y1 , y2 ) dy2 P3 (y1 , y2 , y3 ) dy2dy3 + . . . . (A.3) The presence of the factorials in the denominators comes again from the unlabeled character of the particles. Similarly, we have Fn (y1 , . . . yn ) = ∞ X l=0 1 l! Z Pn+l (y1 , . . . yn , yn+1, . . . , yn+l )dyn+1 . . . dyn+l . (A.4) Next we introduce the quantities fn = Z Fn (y1 , y2 , . . . , yn )dy1dy2 . . . dyn . (A.5) Integrating (A.3) over dy1 , we obtain f1 = p1 + 2!p2 + 3! m! p3 + . . . pm + . . . . 2! (m − 1)! (A.6) More generally, integrating (A.4) over dy1 , . . . , dyn , we have fn = n!pn + (n + 1)!pn+1 + . . . m! pm + . . . . (m − n)! (A.7) Multiplying (A.7) by (−)n /n! and summing over n, we obtain (−)n fn . p0 = n=0 n! ∞ X (A.8) Let us briefly digress and consider the case of independent particles, when Fn (y1 , . . . , yn ) = F1 (y1 )F1 (y2) · · · F1 (yn ), and thus fn = f1n . From (A.8) we then obtain p0 = exp (−f1 ) . (A.9) 49 R Since f1 = F1 (y1 )dy1 is the mean number of particles, (A.9) may be viewed as just a consequence of the Poisson distribution. It is referred to in this paper as the “Poisson approximation”. The general case, of correlated particles, is handled by introducing correlation functions through a cluster expansion : F1 (y1 ) = C1 (y1 ) F2 (y1 , y2 ) = C1 (y1 )C1 (y2 ) + C2 (y1 , y2 ) .. . X Fn (y1 , y2, . . .) = · · · Cl (yr1 , . . . , yrl ) · · · all partitions (A.10) (A.11) (A.12) We also need the cumulants cn = Z Cn (y1 , y2 , . . . , yn )dy1dy2 . . . dyn . (A.13) Note that the cn s for n > 1 vanish in the case of independent particles. The name “cumulant” for the cn s is justified by the observation that the cn s and the fn s defined in (A.5) are related like moments and cumulants : if we introduce the generating function Ω(z) = zn fn , n=0 n! ∞ X (A.14) we easily obtain from (A.10)-(A.12) zm Ω(z) = exp cm . m=1 m! ∞ X ! (A.15) Expansion of the r.h.s. of (A.15) in powers of z and identification with the expansion (A.14) allows the successive determination of the set of fn s in terms of the set of cn s and, conversely, identification of the expansion of the logarithm of the r.h.s. of (A.14) with the logarithm of (A.15) gives the cn s in terms of the fn s. In particular, we have c1 = f1 = Z F1 (y1 ) dy1, c2 = f2 − f12 = Z (A.16) (F2 (y1 , y2 ) − F1 (y1 )F1 (y2 )) dy1 dy2 . 50 (A.17) Comparing (A.8) and (A.14), we have p0 = Ω(−1). Then, using (A.15), we obtain the required expression for the probability of having zero particles : (−)n p0 = exp cn . n=1 n! ∞ X ! (A.18) Of course, if we truncate the series to its first term, we recover the Poisson approximation. References [0] Albeverio S., Molchanov A.A. & Surgailis D. 1994. Stratified structure of the Universe and Burgers’ equation – a probabilistic approach, Prob. Theory Relat. Fields 100, 457–484. [0] Aurell E., Gurbatov S.N. & Simdyankin S.I., 1996. Numerical proof of self-similarity in Burgers’ turbulence, TRITA-PDC Report 1996:2, patt-sol/9602005. [0] Aurell E., Frisch U., Noullez A. & Blank M. 1997. Bifractality of the Devil’s staircase appearing the Burgers equation with Brownian initial velocity, J. Stat. Phys., to appear. [0] Avellaneda M. & E W. 1995. Statistical properties of shocks in Burgers turbulence, Comm. Math. Phys. 172, 13–38. [0] Balian R. & Schaeffer R. 1989. Scale-invariant matter distribution in the universe, Astron. Astrophys. 220, 1–29. [0] Barabási A.-L. & Stanley H.E. 1995. Fractal Concepts in Surface Growth. Cambridge University Press, Cambridge. [0] Batchelor G.K. & Proudman I. 1956. The large-scale structure of homogeneous turbulence, Philos. Trans. Roy. Soc. 248, 369–405. [0] Bleistein N. & Handelsman A. 1975. Asymptotic Expansions of Integrals, Holt, Rinehart and Winston. Reprinted by Dover (1986). [0] Bouchaud J.–P. & Mézard M. 1997. Turbulent renormalisation group in disordered systems. Preprint. [0] Bouchaud J.–P., Mézard M. & Parisi G. 1995. Scaling and intermittency in Burgers turbulence. Phys. Rev. E 52, 3656–3674. 51 [0] Burgers, J.M. 1939. Mathematical examples illustrating relations occurring in the theory of turbulent fluid motion, Kon. Ned. Akad. Wet. Verh. 17, 1–53 (also in ‘Selected Papers of J.M. Burgers’, eds. F.T.M. Nieuwstadt & J.A. Steketee, pp. 281–334, Kluwer, 1995). [0] Burgers J.M. 1974. The Nonlinear Diffusion Equation. D. Reidel, Dordrecht. [0] Cheklov A. & Yakhot V. 1995. Kolmogorov turbulence in a randomforce-driven Burgers equation: anomalous scaling and probability functions, Phys. Rev. E 52, 5681–5684. [0] Cole J.D., 1951. On a quasi-linear parabolic equation occurring in aerodynamics, Quart. Appl. Math. 9, 225–236. [0] E W., Khanin K., Mazel A. & Sinai Ya.. 1997 Probability distribution functions for the random forced Burgers equation, Phys. Rev. Lett. 78, 1904–1907. [0] Esipov S.E. & Newman T.J. 1993. Interface growth and Burgers turbulence : the problem of random initial conditions, Phys. Rev. E 48, 1046–1050. [0] Esipov S.E. 1994. Energy decay in Burgers turbulence and interface growth : the problem of random initial conditions. II, Phys. Rev. E 49, 2070–2081. [0] Fournier J.D. & Frisch U. 1983b. L’équation de Burgers déterministe et statistique, J. Méc. Théor. Appl. (Paris) 2, 699–750. [0] Frisch U. 1995. Turbulence: the Legacy of A.N. Kolmogorov. Cambridge University Press. [0] Funaki, T., Surgailis, D. & Woyczynski, W.A. 1995. Gibbs–Cox random fields and Burgers turbulence, Ann. Appl. Probab. 5, 461–492. [0] Gurarie V. & Migdal A. 1996. Instantons in Burgers equation, Phys. Rev. E 54, 4908–4914. [0] Gurbatov S.N. & Malakhov A.N. 1977. Statistical characteristics of random quasi monochromatic waves in nonlinear media, Sov. Phys. Acoust. 23, 325–329. [0] Gurbatov S.N., Malakhov A.N. & Saichev A.I. 1991. Nonlinear Random Waves and Turbulence in Nondispersive media : Waves, Rays, Particles. Manchester University Press, Manchester. 52 [0] Gurbatov S.N. & Saichev A.I. 1981. Degeneracy of one-dimensional acoustic turbulence for large Reynolds numbers, Sov. Phys. JETP 80, 589–595. [0] Gurbatov S.N. & Saichev A.I. 1984. Probability distribution and spectra of potential hydrodynamic turbulence, Radiophys. Quant. Electr. 27, 303–313; Izv. VUZ Radiofiz. (USSR) 27, 456–468. [0] Gurbatov S.N., Saichev A.I. & Shandarin S.F. 1989. The large-scale structure of the Universe in the frame of the model equation of non-linear diffusion, Month. Not. R. astr. Soc. 236, 385–402. [0] Hopf E. 1950. The partial differential equation ut +uux = uxx , Comm. Pure Appl. Mech. 3, 201–230. [0] Hu, Y. & Woyczynski, W.A. 1996. Shock density in Burgers’ turbulence, in Nonlinear Stochastic PDEs, T. Funaki & W.A. Woyczynski, eds. 157–165. Springer, Berlin. [0] Kardar M., Parisi G. & Zhang Y.C. 1986. Dynamical scaling of growing interfaces, Phys. Rev. Lett. 56, 889–892. [0] Kármán T. von & Howarth L. 1938. On the statistical theory of isotropic turbulence, Proc. R. Soc. Lond. A 164, 192–215. [0] Kida S. 1979. Asymptotic properties of Burgers turbulence, J. Fluid Mech. 93, 337–377. [0] Kolmogorov A.N. 1941. On degeneration (decay) of isotropic turbulence in an incompressible viscous liquid, Dokl. Akad. Nauk SSSR 31, 538–540. [0] Leadbetter M.R., Lindgren G. & Rootzen H. 1983. Extremes and Related Properties of Random Sequences and Processes. Springer, Berlin. [0] Lin C.C. & Reid W.H. 1963. Turbulent flow, theoretical aspects, in Handbuch der Physik, Fluid Dynamics II, 438–523, eds. S. Flügge & C. Truesdell. Springer, Berlin. [0] Loitsyansky L.G. 1939. Some basic laws for isotropic turbulent flow, Trudy Tsentr. Aero.-Gidrodin. Inst, 3–23. [0] Molchanov S.A., Surgailis D. & Woyczynski W.A. 1995. Hyperbolic asymptotics in Burgers’ turbulence and extremal processes, Comm. Math. Phys. 168, 209–226. 53 [0] Monin A.S. & Yaglom A.M. 1975. Statistical Fluid Mechanics. vol. 2, MIT Press, Cambrige, Mass. [0] Noullez A. & Vergassola M. 1994. A fast algorithm for discrete Legendre transforms, J. Sci. Comp., 9, 259–281. [0] Papoulis A. 1991. Probability, Random Variables, and Stochastic Processes, 3rd edition. McGraw-Hill, New York. [0] Pitman, E.J.G. 1968. On the behaviour of the characteristic function of a probability distribution in the neighborhood of the origin, J. Austral. Math. Soc. 8, 422–443. [0] Polyakov, A.M. 1995. Turbulence without pressure, Phys. Rev. E 52, 6183–6188. [0] Proudman, I. & Reid, W.H. 1954. On the decay of a normally distributed and homogeneous turbulent velocity field, Phil. Trans. R. Soc. Lond. A 247, 163–189. [0] Rice, S.O. 1954. Mathematical analysis of random noise, in Selected Papers on Noise and Stochastic Processes, 133–294, ed. N. Wax. Dover. [0] Shandarin S.F. & Zeldovich Ya.B. 1989. The large-scale structure of the Universe: turbulence, intermittency, structures in a selfgravitating medium, Rev. Mod. Phys. 61, 185–220. [0] She Z.S., Aurell E. & Frisch U. 1992. The inviscid Burgers equation with initial data of Brownian type, Commun. Math. Phys. 148, 623–641. [0] Sinai Ya. 1992. Statistics of shocks in solutions of inviscid Burgers equation, Commun. Math. Phys. 148, 601–622. [0] Surgailis, D. 1996. Intermediate asymptotics of statistical solutions of Burgers’ equation, in Nonlinear Stochastic PDEs, T. Funaki & W.A. Woyczynski, eds. 137–145. Springer, Berlin. [0] Surgailis, D. 1997. Asymptotics of solutions of Burgers’ equation with random piecewise constant data, in Stochastic Models in Geosystems, S.A Molchanov & W.A. Woyczynski, eds., vol. 85 of “The IMA Volumes in Mathematics and its Applications”, 427–441. Springer, Berlin. [0] Tatsumi T. & Kida S. 1972. Statistical mechanics of the Burgers model of turbulence, J. Fluid Mech. 55, 659–675. 54 [0] Vergassola M., Dubrulle B., Frisch U. & Noullez A. 1994. Burgers’ equation, Devil’s staircases and the mass distribution for large-scale structures, Astron. Astrophys. 289, 325–356. [0] Whitham G.B. 1974. Linear and Nonlinear Waves. Wiley, New York. [0] Zeldovich Ya.B. 1970. Gravitational instability: an approximate theory for large density perturbations, Astron. & Astrophys. 5, 84–89. 55 FIGURE CAPTIONS Figure 1: Sketch of the energy spectrum at long times when 1 < n < 2. Figure 2: Intersections of the initial potential ψ0 (y) (curve G1 ) with the parabola G2 . Figure 3: Intersections of the initial potential ψ0 with two parabolas. Figure 4: The different regions of integration appearing in (106). Near the diagonal y1 = y2 we have regions I and I’ of (107). Away from the diagonal we have regions II and II’. Figure 5: Evolution of the energy spectrum with an initial spectrum (dashed line) proportional to k n (n = 1.7) at small wavenumbers k. Resolution N = 221 . Spectra averaged over 103 realizations. The labels correspond to output times ti = 10i−6 N 2 , (i = 1, 2, . . . 8). Figure 6: Energy spectrum at t = t3 = 10−3 N 2 . Same conditions as Fig. 5. Initial spectrum : dashed line. Leading-order asymptotic theory in inner region : dotted line. The vertical line shows the switching wavenumber. Figure 7: Evolution of the energy spectrum compensated by a k −n factor in order to reveal the region of “permanence of large eddies” (PLE) as a plateau with unit height. Spectral exponent n = 1.5. Resolution N = 220 . Spectra averaged over 30,000 realizations. The unit of time is different from the other figures, the spatial period being here unity. Observe the progressive shrinking of the PLE zone. (Courtesy A. Noullez.) Figure 8: Evolution of the computed integral scale L(t) (asterisks) compared to the theoretical leading-order prediction (solid line) and the same without the logarithmic correction (dashed line). Same conditions as Fig. 5. Figure 9: Evolution of the coefficient A(t) in front of the k 2 term at intermediate wavenumbers (asterisks), compared to the leading-order theoretical prediction (117) (solid line) and the same without its logarithmic correction (dashed line). Same conditions as Fig. 5. Figure 10: Evolution of the switching wavenumber ks (t), compared to the leadingorder theoretical prediction (119) (solid line) and the same without its logarithmic correction (dashed line). Same conditions as Fig. 5. Figure 11: Evolution of the switching wavenumber ks (t) for n = 1.5. Otherwise, everything as in Fig. 10. 56 View publication stats