Probability&incompressible Navier-Stokes equations: An overview of some recent developments

This is largely an attempt to provide probabilists some orientation to an important class of non-linear partial differential equations in applied mathematics, the incompressible Navier-Stokes equations. Particular focus is given to the probabilistic framework introduced by LeJan and Sznitman [Probab. Theory Related Fields 109 (1997) 343-366] and extended by Bhattacharya et al. [Trans. Amer. Math. Soc. 355 (2003) 5003-5040; IMA Vol. Math. Appl., vol. 140, 2004, in press]. In particular this is an effort to provide some foundational facts about these equations and an overview of some recent results with an indication of some new directions for probabilistic consideration.


Navier-Stokes Equations: Background
The present paper is an attempt to provide a brief orientation to Navier-Stokes equations from a probabilistic perspective developed in the course of working with a focussed research group in this area during the past few years at Oregon State University (OSU). An effort is made to furnish a bit more background for the uninitiated on some of the basics of these equations along with a summary description of the probabilistic methods and results. The approach is in somewhat "broad strokes". The reader should be able to follow and/or supply most basic calculations, but the more technical proofs of some of the main results are left to original references. More detailed accounts of some of the mathematical intricacies, conjectures and open problems identified by the OSU group can be found in Bhattacharya et al. [6,7], Chen et al. [12], and in a recently completed PhD thesis by Orum [45].
A purely mathematical description of the Navier-Stokes equations naturally depends on the framework (function space) in which initial data and external forcings are provided, and in which solutions are sought, e.g. see see Ladyzhenskaya [32,33], Temam [52]. However, rather than begin with the technical nomenclature, let us first take a look at the underlying physical picture. The 3-dimensional incompressible Navier-Stokes (NS) equations governing fluid velocities u(x, t) = (u 1 (x, t), u 2 (x, t), u 3 (x, t)), x ∈ R 3 , t ≥ 0, with initial data u(x, 0 + ) = u 0 (x), are a manifestation of Newton's law of motion and mass conservation for a homogeneous (constant density) fluid: where ∇ = (∂/∂x 1 , ∂/∂x 2 , ∂/∂x 3 ), u · ∇ = 3 j=1 u j ∂ ∂xj , and ∆ = ∇ · ∇ = 3 j=1 ∂ 2 /∂x 2 j ; the latter two are applied component-wise to vector-valued functions. The term p(x, t) is the (scalar) pressure, g(x, t) represents external forcing, and ν > 0 is the kinematic viscosity. Probabilists might want to write the term ν∆ as 1 2 2ν∆ in anticipation of the heat kernel in the form k(x, t) = (4πνt) − 3 2 exp{− |x| 2 4νt } more familiar in applied mathematics. The incompressibility condition refers to the requirement that ∇ · u = 0. ( The (in-viscid) equations corresponding to ν = 0 are referred to as the Euler equations, and are generally considered apart from Navier-Stokes. The equations (1) comprise a (non-linear) system of four partial differential equations governing four unknown scalar quantities; three coordinates of velocity u 1 , u 2 , u 3 and a scalar pressure p. The equations model the motion of a viscous incompressible homogeneous fluid in an inertial frame in free space. They are derived from conservation of mass and momenta of "fluid elements" (in the sense of continuum mechanics) under the additional assumption that there is a local linear relation between the stress and strain rates experienced by fluid parcels. In (1), the left-hand side ∂u ∂t + (u · ∇)u represents the acceleration in an Eulerian frame of reference; i.e. one views the motion as that of a "parade of fluid elements" and u(x, t) is the velocity of the fluid parcel at that location and at that time. Thus if a parcel has moved to location x(t) = (x(t), y(t), z(t)) at time t, then the velocity of the parcel is given by u(x, t) ≡ (u 1 (x, y, z, t), u 2 (x, y, z, t), u 3 (x, y, z, t)) = ( dx dt , dy dt , dz dt ) and, using the chain rule, the acceleration is therefore given by the so-called material or convective derivative The right side of (1) resolves the forces acting on the fluid parcel. The first term ν∆u represents the viscous friction forces and is often referred to as the dissipative term. This term arises in the fluid mechanical derivation as the result of the aforementioned ad-hoc model of linearized shear stresses due to interactions between molecules comprising the fluid. The term g ≡ g(x, y, z, t) represents body forces, such as gravity and/or other externally applied force fields. The pressure gradient term −∇p defines the local pressure force on a fluid parcel (depending on p up to an arbitrary constant). Physically the pressure gradient is a force which acts to maintain a homogeneous fluid density under incompressibility. Thus one may reasonably expect that pressure can be computed from the velocity and boundary conditions (at infinity in this case). We will see that the following orthogonality opens a way to mathematically first remove the pressure for focus on the incompressible velocity. Integrating the following "product formula" over a ball B 0 (R) of radius R ∇ · (pu) = u · ∇p + (∇ · u)p (4) and then applying Gauss' divergence theorem yields the following version of an "integration by parts" type formula Thus, if u is incompressible then the second integral on the right side of (5) vanishes and the first integral is O( pu1(B 0 (R)) ∞ R 2 ); here 1(B 0 (R)) denotes the indicator function. Assuming pu decays sufficiently fast at infinity, one therefore obtains the following basic L 2 −orthogonality R 3 u · ∇p = 0.
Remark. Decay assumptions are implicit to most of the theory to be described throughout this article, from orthogonality conditions (6) to Fourier transform formulations. Such boundary conditions at infinity on Navier-Stokes are the subject of extensive investigation for both free-space and exterior domain problems since Finn's conception of the notion of "physically reasonable" solution; e.g. see Finn [18], Heywood [28], Galdi [23]. A perhaps more familiar illustration of the need to impose decay conditions occurs for example in the simple case of the heat equation, where one observes that since the Gaussian heat kernel is a tempered function, its convolution with a tempered distribution is well-defined as a function. In particular, it is sufficient to assume that the initial data is a tempered distribution for the well-known solution as a convolution with the heat kernel to apply. In summary, two distinguishing features of the equations (1) are: (a) the nonlinear term, referred to as the convective term, (u · ∇)u on the left-hand side, whose role in representing the acceleration is intrinsic and cannot be dropped by changing the model of the forces on the right-hand side, and (b) the incompressibility condition ∇ · u = 0. This latter condition is an additional equation to momentum balance obtained by consideration of mass conservation for a homogeneous fluid of constant fluid density.
Systematic derivations of the Navier-Stokes equations can be found in any number of standard texts on fluids, e.g. Batchelor [3], Landau and Lifschitz [34], Chorin and Marsden [15]. Thumbnail historical sketches of Claude Louis Marie Henri Navier (1785-1836) and George Gabriel Stokes (1819-1903) are available on the web. A brief textual history as it pertains to fluid flow as well as a fundamentally significant role in engineering design can be found in Anderson [1].
Let us briefly turn to the mathematical framework in which one might view these equations. Cannone [9] contains an excellent survey of a variety of natural function spaces in which to analyze the equations. However, rather than attempt to develop the precise function space requirements at the outset for the results to be discussed, we will for the most part continue to confine ourselves to L 2 (R 3 ) 3 and/or spaces of solenoidal, i.e. ∇ · u = 0, tempered distributions u, where the reader will be able to check basic statements and calculations using standard Fourier methods as covered, for example, in Folland [19]. The actual function spaces will evolve naturally in the course of the development of the probabilistic approach.
An in-depth analysis of these equations given by Jean Leray [36] 1 has become a lasting benchmark for contemporary research. Under the assumptions that u 0 , g ∈ L 2 (R 3 ) 3 , ∇ · u 0 = 0, and an energy inequality of the form applies (in the distributional sense), Leray [36] proves the existence of a global (weak) solution to (1) satisfying the inequality (7) . The corresponding uniqueness in L 2 (R 3 ) 3 of Leray's solution remains open to this day for the case of n = 3 space dimensions. From physical considerations, namely conservation of energy, the energy inequality (7) is expected to be an equality. However it has not been proven that such a solution exists in three-dimensions. A formal derivation of conservation of energy is obtained by taking the scalar product of the terms in (1) with u and then integrating over R 3 . Assuming sufficiently rapid decay of u, p at infinity, integration by parts and incompressibility yield that the convective term and pressure both vanish. For the scalar product with the dissipative term one has similarly by performing integration by parts Noting that 1 2 d dt R 3 |u(·, t)| 2 ≡ ∂u ∂t · u one has the equation for conservation of energy 1 2 At another extreme for solution spaces, in defining the Millenium Prize Problem for the Clay Institute, Fefferman [17] stipulates that a physically reasonable solution must be infinitely smooth and decrease more rapidly than any power of the spatial variable x. However this choice of solution space is not without controversy in the international mathematical community. In a critique of this requirement, Ladyszhenskya [33] remarks that "...one must leave both the choice of the phase space and the class of generalized solutions to the researcher without prescribing to him infinite smoothness or some other smoothness of solutions. The only requirement needed is indeed that the uniqueness theorem must hold in the chosen class of generalized solutions." Here she is addressing her basic point that the first fundamental question is whether the equations (and boundary conditions) do uniquely determine the dynamics of an incompressible fluid. In particular, questions regarding the qualitative behavior and properties of the solution naturally follow. In the context of this article boundary conditions refer to decay properties at infinity, however the methods extend to periodic boundary as well. 2 Apart from definitions of classical and various forms of weak (in distributional sense) solutions, some of the mathematical intricacies are further illustrated with the work of Fujita and Kato [21], Kato [30] in which a notion of "mild solution" is advanced as follows: The so-called Leray projection operator P is most easily defined by the orthogonal projection of the Hilbert space the subspace of divergence-free (incompressible) vector fields corresponds to the following basic orthogonality in Fourier frequency space Now, since orthogonality of the projection requires thatû(ξ)− Pu(ξ) be orthogonal to Pu(ξ), this vector difference is in the direction of ξ, i.e.û(ξ) − Pu(ξ) = ( ξ |ξ| ·û(ξ)) ξ |ξ| . Thus, where I = ((δ ij )) 1≤i,j≤3 is the identity matrix and the tensor product ⊗ of two vectors a, b ∈ R 3 is defined by the matrix a ⊗ b = ((a i b j )) 1≤i,j≤3 . The projection extends to other function spaces for which Fourier transforms can be defined. In view of the incompressibility condition and (6) one may apply the projection operator to (1) to obtain for incompressible initial data u 0 belonging to a suitable function space, e.g. L 2 (R 3 ) 3 , the "projected equation" If one views (13) as a (vector) heat equation forced by the (unknown) term w = P(g − u · ∇u) then a standard Duhamel argument, e.g. see Bhattacharya and Waymire [4], may be applied to deduce the following integral equation where {T t ≡ e tν∆ : t ≥ 0} is the (self-adjoint) semi-group associated with 3dimensional Brownian motion Z = {Z(t) : t ≥ 0} with zero drift and diffusion coefficient 2ν, i.e. T t g(x) = E x g(Z(t)) = R 3 g(y) 1 (4πνt) 3 2 exp{− 1 4πνt |y − x| 2 }dy, and applied component-wise. Solutions to (14) satisfying ∇ · u = 0 are referred to as mild solutions. This is the starting point of Kato's approach in which he proved existence of mild solutions for "small" initial data in the special function space L 3 (R 3 ) 3 . Specifically, ignoring external forcing for simplicity, Kato [30] showed that there is a number A > 0 such that for u 0 L 3 ≤ A, one has a solution such that t → u(·, t) ∈ L 3 (R 3 ) 3 is continuous and satisfies (14) on every compact time interval [0, T ], T > 0. However it wasn't until very recently that uniqueness was proven by Furioli, Lemarié-Rieusett and Terraneo [22].
The Picard iteration scheme provides a useful technique naturally associated with the (projected) Navier-Stokes equation. Observe that using incompressibility the non-linear term may be expressed, analogously to uu x = 1 2 (u 2 ) x in one-dimension, in the equivalent form Then the iteration is given by where The convergence of the iterates follows from showing that that these are iterates of a contraction on a ball in a suitable function space; c.f. Meyer (2004) and references therein.
It may be noted that the function space L 3 (R 3 ), or more generally L n (R n ), enjoys the very special norm scale invariance of the form Such Banach spaces are referred to as limit spaces. Excellent surveys of contemporary literature which grew out of Kato's approach are given in Meyer [42], Cannone [10,9]. A more general notion of limit space arises in the context of majorizing kernels and will be given in Section 3; see (64). The Picard iteration approach may also be viewed as a starting point for the probabilistic methods initiated by LeJan and Sznitman [37], though perhaps not in a way one might expect on first glance. Their's is based on a branching representation in Fourier space in which the non-linearity in the Fourier-transformed version of (13) is associated with a "binary branching product" of certain random walks in Fourier frequency space. The essential difference in the two approaches is in the nature of bounds sought for convergence of two distinct but closely related limits. In the former one seeks contraction bounds for convergence of the iteration scheme, while in the latter one seeks integrability bounds on stochastic variables, defined by a.s. limits of a stochastic iteration, to ensure finiteness of their expected values. To highlight the close connections between these two approaches it may be noted that the aforementioned uniqueness proof by Furioli, Lemarié-Rieusett and Terraneo [22] for Kato's solution in L 3 (R 3 ) was obtained by arguments, specifically estimates on the associated bilinear form b(u, v), which the authors explicitly acknowledge to have been inspired by the probabilistic approach to uniqueness for a certain other space obtained by LeJan and Sznitman [37]; namely the space F = {u ∈ S ′ (R 3 ) : sup ξ |ξ| 2 |û(ξ)| < ∞} or, more accurately, the solenoidal distributions belonging to this space F . A distribution whose Fourier transform is essentially bounded is referred to as a pseudo-measure in harmonic analysis. Thus an alternative description of this space is one of locally integrable tempered distributions for which ∆u is a pseudo-measure; now often denoted PM 2 where the superscript 2 refers to the |ξ| 2 . pre-factor, c.f. Cannone [9]. In fact this may also be shown to be a Besov space; see Cannone and Planchon [11]. For the unitiated, a brief historical orientation to the various such function spaces which arise in this context may be found in Peetre [48]. Since, loosely speaking, algebraic products transform to convolution products of Fourier transforms and since derivatives transform by simple algebraic factors, the transformation of the non-linearity under Fourier transform in the (quasilinear) Navier-Stokes equation is deceptively similar to the convolution term one obtains in Fourier space for (semi-linear ) reaction-diffusion equations, e.g. KPP. But one must not forget about the incompressibility condition, which is completely absent for the case of reaction-diffusion. Nonetheless, it is natural to search for branching Brownian motions in physical space; e.g. see McKean [40] in the case of KPP. In the case of the scalar Burgers equation, again devoid of a notion of incompressibility, one may identify such a representation; see Waymire [56], Ramirez [49]. Moreover, recently a probabilistic representation involving semi-Markovian branching Brownian motions in physical space was obtained by Ossiander [47] in which some basic notions from probabilistic potential theory make a natural appearance. In the next section we will further indicate natural probabilistic ties to harmonic analysis in the way of the "background radiation" process in physical space introduced in Gundy and Varapoulos [27] and Gundy and Silverstein [26].
The previously noted basic results capture the essence of the paradigm furnished by iterative methods in a wide range of function spaces. Namely, one either obtains for sufficiently small initial data u 0 existence and uniqueness of global solutions, i.e. for all times T > 0, or one proves existence and uniqueness for arbitrarily large initial data u 0 for short durations T ≡ T (u 0 ) < ∞. A concisely formulated point for caution concerning semi-group approaches was provided in a paper by Montgomery-Smith [44] which suggests a much deeper role for incompressibility than has thus far been captured by these methods. However there are many significant questions beyond existence and uniqueness, e.g. numerical computation and Monte-Carlo simulation algorithms, boundary value problems, and steady-states, which stand to benefit by continued development of probabilistic methods. It is manifestly clear that the introduction of probabilistic reasoning has certainly lead to new insights into the structure of these equations and their solutions. On the other hand, it is also clear that incompressibility is being undervalued in these analyses to date.

Vorticity, Pressure, Incompressibility and Background Radiation
One need do no more than point to the classic papers on diffusions by W. Feller, to name just one example, to highlight and exemplify the impact that the analytic theory of linear (parabolic) partial differential equations has had on the development of probability theory. On the other hand there is no question to the significant counterpoint that the subsequent development of martingale theory, the Itô calculus and the Malliavin calculus has provided in the study of such pde's. However, unlike the fundamental equations governing diffusion, the historic origins of the Navier-Stokes equations are devoid of explicit probabilistic reasoning for their formulation -making the appearance of probabilistic methods here all the more mysterious from the outset.
There is that Laplacian term! One of the earliest applications of probabilistic ideas in connection with incompressible Navier-Stokes was introduced by Chorin [13], based on the associated vorticity equation. The vector cross product ω(x, t) = ∇ ∧ u(x, t), i.e. the curl of the velocity field, defines the vorticity. Although the physical equations are derived in n = 3 dimensions the formulation extends to n ≥ 2 dimensions. In n = 1 dimension the notion of vorticity is lost; the corresponding scalar equation for velocity ∂u ∂t = ν ∂u ∂x + u ∂u ∂x + g being referred to as Burgers equation. However in the lowest dimensional case (n = 2) of velocity u = (u 1 , u 2 ) in (1), one may also identify the vector vorticity ω ≡ (0, 0, ω) with the scalar component ω.
Gradients are annihilated by taking the curl since for any sufficiently smooth scalar function a In particular taking the curl provides another approach to remove the (unknown) pressure gradient in (1). Applying some simple vector identities, e.g. see Chorin and Marsden [15, pp. 163-164], namely (u · ∇)u = 1 2 ∇(|u| 2 ) + ω ∧ u and ∇ ∧ (ω ∧ u) = ω(∇ · u) − u(∇ · ω) + (u · ∇)ω − (ω · ∇)u, and using incompressibility and (18), one may express the curl of u satisfying (1) as A solution ω of (19) may then be used to obtain a solution u to (1) and, in turn, the pressure p as follows: According to the Helmholtz-Hodge decomposition any smooth vector field u which decays sufficiently fast at infinity may be uniquely represented as the superposition of a gradient and a (necessarily incompressible) curl where for a scalar potential Φ and vector potential Ψ obtained by solving a (vector) Poisson equation In particular one may recover u from vorticity ω via the following so-called Biot-Savart law: where L is a linear (convolution) kernel furnished by the divergence free solution Ψ to the Poisson equation, i.e. by inverting the Laplacian via (for n = 3) and hence (for n = 3) Similarly for n = 3 one needs only to apply the appropriate potential kernel to invert the Laplacian. In particular for n = 2 dimensions one has where For numerical computations Chorin [13] viewed (19) by an operator splitting, i.e. along the lines of a Trotter-Kato semi-group product, into successive incremental advective and diffusive flow processes using distinct numerical schemes.
In this connection, assuming n = 2 observe that the term ω · ∇ in (19) vanishes for n = 2; here one should recall the convention of using the scalar ω and vector (0, 0, ω) interchangeably in n = 2 dimensions. Thus, assuming also g ≡ 0, the two-dimensional vorticity equation describes the evolution of a scalar Now, if u were given (computed) over an increment of time, then the scalar component of vorticity would appear to be governed by a diffusion equation.
In this regard, denote two-dimensional standard Brownian motion by : t ≥ 0} be the diffusion defined by the associated stochastic differential equation Then one is naturally led to consider This is the starting point for an approach which has been exploited heavily as a numerical Monte-Carlo method, and has been studied closely by analysts and probabilists alike; e.g. Chorin [13], Marchioro and Pulverenti [39] Goodman [24], Long [38], Szumbarski and Wald [51], Meleard [41]. In this regard, Busnello, Flandoli, and Romito [8] also obtain a natural probabilistic version of the Biot-Savart law in three-dimensions. While we will not pursue these particular approaches any further, we will next take the idea of considering linearization and Duhamels principle in an apparently new direction which exploits incompressibility in an interesting way.
To recover the pressure from its projection onto a divergence free vector field let us recall (12). The matrix elements ξi |ξ| ξj |ξ| may be be viewed as Fourier symbols of iterated Riesz transforms R i R j , i.e. R j f (ξ) = − ξj |ξ|f (ξ); note the signs and constants depend on the version (10) of Fourier transform used. In particular, on appropriate function spaces, e.g. L 2 (R 3 ) 3 , one has where R = ((R i R j )) 1≤i,j≤3 is the matrix of iterated Riesz transforms. To exploit incompressibility in (1) note that the terms linear in u are divergence free. Thus take the divergence to annihilate the linear terms ∂∇·u ∂t = 0 = ∆(∇ · u) = ∇ · (∆u), and assume no forcing g = 0, so that one has Taking Fourier transforms with continued use of incompressibility one arrives at the following determination of the pressure Apart from what may be gleaned from its potential theoretic role in making projections onto divergence free vector fields, a basic understanding of incompressibility remains illusive from a probabilistic perspective. However, even in the case of linearized Navier-Stokes equations, the incompressibility condition plays a fascinating role. Two common linearizations are (i.)Stokes equations and (ii.)Oseen equations, respectively given by linearization around zero velocity or by linearization around a non-zero constant U, respectively. Specifically one has, respectively, for Stokes and for Oseen For simplicity and without any significant loss of generality let us consider the Stokes equation (34). Observe that for given twice continuously differentiable and incompressible initial data u 0 = (u 0 ), one has that in the unforced case (g ≡ 0) the solution (u, p) = ((u 1 , u 2 , u 3 ), p) is given by the familiar semi-group representation is three-dimensional standard Brownian motion, and c is an arbitrary constant. In particular the linearity in u and the apparent affine form of (34) due to −∇p are reconciled by the incompressibility. In other words, the incompressibility of u (at times t > 0) is directly inherited from that of u 0 , causing −∇p to instantaneously "drop out" at t = 0. Accordingly one may view the matrix-valued Gaussian (heat) kernel furnished by the transition probabilities of the Brownian motion, as the fundamental solution to (34) for incompressible initial data. Let us consider, on the other hand, the evolution of smooth but compressible initial data u 0 . Also let us include a possibly compressible forcing term g. We seek an incompressible vector field u(x, t) ≡ Pu(x, t), t > 0, such that by orthogonality (6) of ∇p to incompressible u for t > 0, The compressible initial data u 0 will necessitate that u(x, t) ≡ Pu(x, t) for t > 0, be discontinuous at t = 0. Applying Duhamel's principle to (38) one where and R = ((R j R k )) 1≤j,k≤3 is the matrix of iterated Riesz transforms. The simplest way to check the commutativity required for the second equation in the display (39) is by Fourier transforms along the lines indicated earlier for (31).
To compute the pressure one simply takes the divergence in (34) and use incompressibility to get For t = 0 one may apply the Helmholtz-Hodge decomposition to recover u 0 from Pu 0 defining the pressure at t = 0 by In particular, therefore, One may note that (39) actually embodies two Duhamel principles; one applied and one derived! Namely, we apply a Duhamel principle for Brownian motion to (38) with P(g) as the forcing, and derive a Duhamel principle for (34) with arbitrary initial data using the fundamental solution (40) of the unforced linear incompressible pressurized equation and adding g as the forcing. The representation (40) of the (signed) fundamental solution Γ of Stokes equation for (u, p) is the well-known classic formula of Oseen [46]; also see Solonnikov [50], Thomann and Guenther [53]. One may check a semi-group property of the form Moreover, Γ has the following equivalent representation where Hess denotes the Hessian matrix. A simple way to check this is via Fourier transforms since ξj ξ k |ξ| 2 is the Fourier symbol for the operator ∂ 2 ∂xj x k ∆ −1 . Thomann and Guenther [53] recently made an explicit computation of the full fundamental solution in terms of special functions which has had major implications for analysis of Stokes equation in unrestricted domains. Their precise formulae may be expressed as where Erf(s) = 2 √ π s 0 e −y 2 dy denotes the error function and , is a confluent hypergeometric function. This is a hugely significant formula in view of the rich mathematical literature on parametric dependencies, asymptotic growth and/or decay properties etc of the special functions which appear in the formula (46).
The matrix Γ(x, 0) is the orthogonal projection of the (vector) Dirac distribution δ x onto divergence-free vector fields. Thus, the effect of incompressibility is instantaneous in that the solution is obtained by projecting the initial data onto divergence-free vector fields and then solving (34) with incompressible initial data as above. In fact one may observe that The transformation Γ(x 0 , 0)u 0 is the (instantaneous) projection of u 0 onto its divergence free (incompressible) component. The sign changes and discontinuity at t = 0 in the fundamental solution make the issue perplexing from a probabilistic standpoint. However, there is some hope that a probabilistic unraveling of this semi-group might lead to ways to deal with boundaries. Thomann and Guenther (personal communication) have pointed out that if one poses a Neumann boundary condition on the halfspace, the usual "reflection method" applied to the fundamental solution Γ for the unrestricted problem fails.
In view of the prominent role of Riesz transforms it seems natural to explore a role for Gundy-Varapoulos-Silverstein's background radiation process in the present framework. Thus we close this section with elements of such a new connection. The background radiation process Z = {Z t ≡ (Y t , X t ) : −∞ < t ≤ 0} belongs to the class of "approximate Markov processes" defined, as a collection of measurable functions, on a sigma-finite measure spaces (Ω, F , P ) taking values in the half-space H 3 = [0, ∞) × R 3 . Foregoing its construction, the process is defined by the following properties (Gundy [25]) : is distributed as a one-dimensional Brownian motion on (0, ∞) started at (arbitrary) y 0 > 0 and killed upon entry to the state 0. 2 lim t↑0 Y t = 0 a.s.
is an independent Brownian motion with initial distribution (at time "t = −∞") given by Lebesgue measure on R 3 in the sense that for each y > 0, defining τ y = inf{t ≤ 0 : Y t = y}, the process {(y, X τy+t ) : 0 ≤ t ≤ −τ y } is distributed as a Brownian motion started with Lebesgue measure on the hyper-plane {y} × R 3 at t = 0. 4 X 0 = lim t↑0 X t exists a.s. and is distributed as Lebesgue measure on R 3 , i.e. the induced measure P • X −1 0 is Lebsesgue measure. These properties were used by Gundy and Varapoulos [27] in Gundy and Silverstein [26] to prove the following formula for iterated Riesz transforms of a function f belonging to Schwartz space S where f 0 denotes the harmonic extension of the function f on the R 3 to the half-space H 3 , the matrix A jk = e j ⊗ e k , i.e. having 1 in location (j, k) and 0's elsewhere, defines the indicated martingale transform. With a little more work one may in fact calculate b = − 1 2 using their indicated methods. In view of (40) and viewing conditional expectations as disintegration formulae, one may express the solution to Stokes linearization (34) in terms of background radiation as follows, writing u = (u (1) , u (2) , u (3) ): where {X * t : t ≥ 0} is the reversed-time process (Brownian motion) introduced in Gundy and Silverstein [26], and θ (jk) is the martingale transform defined by and for fixed t, x, K 0 is the harmonic extension of the the heat kernel z → K(z− x, t) to the half-space H 3 . In particular, defining the matrix-valued function one has the representation Form here one may regard the non-linear term in (1) as a forcing and exploit a Duhamel principle to obtain an integral equation involving background radiation. While this connection has not been pursued beyond this survey article, an alternative real-space semi-Markov cascade representation of Navier-Stokes is fully developed in Ossiander [47] which also contains intriguing connections with the fundamental solution obtained by Thomann and Guenther [53] for linearized flow.

Multiplicative Stochastic Cascade
One may give purely fluid mechanical derivations of diffusion equations based on mass conservation and linear flux laws such as Ficks and Fourier laws, however the most compelling physical arguments are no doubt Einstein's implicit use of the central limit theorem and law of large numbers. While the modern approach of hydrodynamic scaling limits of interacting particle systems seeks to arrive at the Navier-Stokes equations from a more intrinsic probabilistic or statistical mechanical foundation , such approaches are fundamentally much more complex than the connections that are possible between Brownian motion and the heat equation; see e.g. Esposito, Marra and Yau [16], Yau [55]. That said, owing to a remarkable insight of LeJan and Sznitman [37] it is nevertheless possible to approach certain classes of solutions to these equations by probabilistic constructions of branching random walks and multiplicative cascades in a way which preserves some basic Markovian semi-group structure. This is the topic to be briefly described in this section.
As remarked earlier the multiplicative cascade approach appears as a mathematical device, closely related to Kato's contraction map techniques for obtaining a mild solution on appropriate function spaces. While a pure analyst may take solace in this realization, as a probabilistic technique it has additional virtues. Chief among these is the "natural" emergence of solution spaces (and their attendant norms) for which one may obtain unique global solutions under a "small" initial data constraint as a consequence of the probabilistic representation. The calculation of steady state solutions are also quick to follow in these function spaces; see Bhattacharya et al. [7]. New questions and ideas with regard to Monte Carlo methods emerge naturally in more pragmatic considerations of numerical solutions; see Ramirez [49]. However, whether this approach can lead to truly significant new results on the existence and uniqueness of global solutions under smooth initial data without regard to "size" remains to be seen. It is clear that a more substantial use of incompressibility is lacking in the existing theory; c.f. Montgomery-Smith [44]. Example (Warm-up). Some of the basic ideas under multiplicative cascade representations are more generally applicable to diverse classes of evolution equations, including certain linear parabolic and fractional diffusion equations, semi-linear reaction-diffusions and some quasi-linear equations such as twodimensional incompressible Navier-Stokes equations and one-dimensional Burgers' equation; see Bhattacharya et al. [7]. Although interesting from a number of other points of view, e.g. physical space representations, Monte-Carlo simulations, etc., these examples are not capable of furthering understanding of incompressibility. That said,, the following extremely simple example serves to illustrate some of the most basic graph theoretic and probabilistic ideas involved in this approach. Consider in n ≥ 1 dimensions, where a > 0, and b ∈ R n are constants. To quickly get the flavor of the method, consider the spatial Fourier transform of (53) Now consider the random linear tree τ θ (t) rooted at a vertex θ of type ξ θ = ξ which, after an exponential length of time is replaced by a single vertex 1 of the same type ξ 1 = ξ. Proceeding in this manner one may calculate that the very familiar solution exp(−a|ξ| 2 + ib · ξ)û 0 (ξ) is the expectation of the random product (θ, t) initialized by ξ θ = ξ and consisting of factors m(ξ) = ib · ξ/a|ξ| 2 at each vertex until termination where one attaches the end factorû 0 (ξ). The indicated stochastic recursion then yields and the solution is represented aŝ where N (t) is the Poisson process with parameter λ(ξ) := a|ξ| 2 which counts the number of times the exponential clocks ring before time t. In particular the Poisson process formally suggests a Fourier dual role to that played by the standard Brownian motion in the real space expectation formula. Similarly one may obtain a formal "dual" Feynman-Kac formula (in the sense of this Fourier transform representation) under the complex measure condition on coefficients given by Itô [29]; see Kolokoltsov [31], and Chen et al. [12]. In particular this approach makes Itô's complex measure condition completely natural from a probabilistic point of view. One may also obtain a Fourier version of McKean's [40] branching Brownian motion formula for KPP, and other interesting equations; see Orum [45]. Also included are, for example, the generalized fractional Burgers equation of the type considered by Woyczynski, Biler, and Funaki [54]. The essential idea introduced by LeJan and Sznitman [37] and presented here in the generalized form developed in Bhattacharya et al. [6,7] is as follows. Without loss of generality assume that the forcing g is incompressible; else replace g by Pg below. Either by directly taking Fourier transforms in the projected equation (13) and then introducing an exponential integrating factor, or by taking Fourier transforms in (14), one may deduce the following integral equation in Fourier frequency space: where for complex vectors w, z is the projection of w onto the plane orthogonal to ξ, and ν > 0 is the viscosity parameter. For ξ = 0, LeJan and Sznitman [37] rescale the equation (FNS) to normalize the integrating factor e −ν|ξ| 2 s to the exponential probability density ν|ξ| 2 e −ν|ξ| 2 s . They then observe that the resulting equation is precisely the form for a branching random walk recursion for χ(ξ, t) := ν|ξ| 2û (ξ, t), for which the kernel |ξ−η| −2 |η| −2 is naturally constrained by integrability to dimensions n ≥ 3 for normalization to a transition probability. Alternatively, Bhattacharya et al. [6] introduce a Fourier multiplier 1/h, where h(ξ), ξ ∈ W h := {ξ ∈ R 3 : ξ = 0} is a positive function such that 0 < h * h(ξ) < ∞, which is used to re-scale the Fourier transformed equation (FNS) by factors 1/h(ξ) in compensation for the temporal normalization of the exponential factor to a probability. Namely, we consider the equation (FNS h ) obtained from (FNS) as and H(ξ, dη 1 × dη 2 ) is for ξ ∈ W h the transition probability kernel, with support contained in the set In view of the conservation of wave-numbers η 1 + η 2 = ξ reflected in the support of the transition probability kernel, the definition of the supporting set of frequencies W h may be relaxed to a semi-group W h ⊆ {ξ ∈ R 3 : ξ = 0}. However in general we include the following additional exterior condition in defining (FNS h ): One may also note that with the exception of ∂û/∂t andĝ, the Fourier transform of each term in (1) vanishes at zero frequency. Thus the value ofû(0, t) ≡ R 3 u(x, t)dx is explicitly determined from the initial data and forcing by The probabilistic approach is based upon an interpretation of the re-scaled integral equation. This is achieved in terms of expectation values of multiplicative cascade solutions to stochastic recursions generated by a multi-type branching random walks in Fourier space. The process may be initiated at the root vertex θ of a binary tree graph having wave-number type ξ = 0. After an exponential holding time S θ with rate ν|ξ| 2 an independent coin flip is made. If the outcome is tail, denoted [κ θ = 0], then the tree terminates as a single edge of length S θ connecting θ and < 1 >. If the outcome is head, denoted [κ θ = 1], then the edge of length S θ branches into two vertices < 1 >, < 2 > of respective types ξ 1 , ξ 2 , ξ 1 + ξ 2 = ξ, selected according to the probability kernel H(ξ, ·). Given ξ two independent exponential clocks with respective rates ν|ξ 1 | 2 and ν|ξ 2 | 2 are reset at the vertices < 1 >, < 2 > and the process is repeated with each of these vertices as "roots". The root θ has an associated birth time B θ ≡ 0 and an end-time S θ at which it is terminated or it branches. Similarly each de- ..,vj > . and an end-time B v + S v . At a given time t > 0 one assigns multiplicative weights to vertices v such that either In the first case the re-scaled forcing ϕ(ξ v ) is assigned, while in the latter case the re-scaled initial data χ(ξ v ) is assigned, where ξ v denotes the frequency type of v; see Figure 1.
Since the underlying discrete parameter branching process defined by the equally likely Bernoulli offspring distribution is a critical binary Galton-Watson process, the recursion will terminate in a finite number of iterations with probability one. In particular, (θ, t) is a.s. a finite product of values of χ 0 and/or ϕ evaluated at randomly selected Fourier frequencies selected by the branching random walk. For example, the functional evaluation of the sample tree in the figure is given by In general (·) is a random multiplicative functional of scalar values m(·) and Fourier transformed initial data and/or forcing (vector) values over the vertices of a multi-type branching random walk tree τ θ (t) initiated from a single progenitor of type ξ θ = ξ. In general the scalar and vector valued factors are evaluated at the wave-number (type) of the respective vertices appearing in the tree τ θ (t), with the initial and forcing terms appearing at the end-nodes. This generalizes branching random walks in the sense of LeJan and Sznitman [37] for n = 3 dimensions where h(ξ) = |ξ| −2 , ξ ∈ W h = R 3 \{0}. The essential requirement for this approach is that the following expected values exist The simplest proof of this uses the strong Markov property to check that if the indicated expected values are finite then (FNS) is satisfied. Existence of such expected values has been obtained by Bhattacharya et al. [6] through constructions of a particular class of the Fourier multipliers, referred to as majorizing kernels, and defined below, from which one also obtains uniqueness in a naturally associated function space. Uniqueness is proven by arguments which exploit discrete parameter martingale structure in successive truncations by generation of the branching random walk product. The point is that uniqueness proofs are largely insensitive to the holding time distribution. In the original development of the theory the holding times between branchings were exponentially distributed with holding time parameters determined from the principal part (Laplacian) of the equation as above. However, recent extensions by Bhattacharya et al. [7], Orum [45], and Ossiander [47] have been obtained based on semi-Markov cascades, which in turn have lead to new representations of both local and global multiplicative cascade solutions to (1). A first order approach to obtain finite expected values of the branching random walk cascade will be seen to result from the observation that the product ⊗ ξ satisfies |w ⊗ ξ z| ≤ |w||z|, w, z ∈ C n , and the coefficients m(ξ) may be controlled by selecting Fourier multipliers such that m(ξ) ≤ 1 or, equivalently, for suitable B > 0. Such h is referred to as a majorizing kernel with constant B. The choice of B is linked to the size of |û 0 | and may be standardized to B = 1 in this context by consideration of suitable initial data. In any case, h is a standardized majorizing kernel if and only if for some, and hence any, B > 0 Bh is a majorizing kernel with constant B.
The following slightly more general definition is suitable for extensions to generalized Navier-Stokes equations with fractional Laplacian, semi-linear equations and for considerations of local solutions. Definition. A positive locally integrable function h on W h ⊆ R n \{0} whose closure W h is a semi-group and such that i h is continuous on W h ii h * h > 0 a.e. on W h iii h * h(ξ) ≤ B|ξ| θ h(ξ), for ξ ∈ W h and some real exponent θ and some B > 0, will be referred to as an FNS-admissible majorizing kernel with constant B and exponent θ. We define h = 0 on W c h and refer to W h as the support of h. Those majorizing kernels h(ξ) which are defined and positive for all ξ = 0 are said to be fully supported.
Formulated in these terms, the results of LeJan and Sznitman [37] may be interpreted in terms of two exponent one, standardized fully supported majorizing kernels, π 3 /|ξ| 2 and αe −α|ξ| /2π|ξ|. These kernels are respectively non-integrable and integrable, with equality in (iii) of the above Definition. The former is the Fourier transform (in the sense of distribution) of the Riesz potential (convolution kernel) defining the integral operator (−∆) −1 , while the latter is the Fourier transform of the Bessel potential defining the integral operator (I − ∆) −1 ; c.f. Folland [19, pp. 149-154].
Apart from their role in existence, uniqueness and expected value representations valid in particular function spaces, the majorizing also play a role in tracking such structure as regularity, analyticity, support size, complexification, etc. of initial data in the evolution of solutions through subspaces corresponding to such more detailed specifications as can be reflected in their Fourier transforms. More specifically, let us define a Banach space F h,γ,T with a norm that depends on a Fourier multiplier 1/h as the completion of the set under the indicated norm, where γ ∈ {0, 1} serves to conveniently index two different norms we wish to consider. Here S ′ is the space of tempered distributions on R n . The definition of majorizing kernel and of the function spaces F h,γ,T imply that the product of distributions in F h,γ,T is itself a distribution. Also, implicit to the definition of the Banach space F h,γ,T is the requirement that tempered distributions belonging to this space have Fourier transforms which are functions. In passing we wish to point out that the notion of limit space referred to in Section 1 has a natural generalization in this context. In particular if h is a majorizing kernel of exponent θ = 1 then h λ (ξ) = λ −2 h(ξ/λ) is also a majorizing kernel of the same exponent. Moreover if u ∈ F h,γ,T then u λ ∈ F h λ ,γ,T and | |u| | F h,γ,T = | |u λ | | F h λ ,γ,T . Thus, while an exponent one majorizing kernel h such that h = h λ defines a limit space F h,γ,T in the usual sense, the norm relation is more general. We will conclude this section with statements of general existence/uniqueness theorems possible within majorizing spaces. The main results use majorizing kernels of different exponents to establish existence, uniqueness and regularity properties of the solutions of the (FNS). Moreover these solutions have the expected value representation in terms of the multiplicative stochastic functional (θ, t) of a multitype branching random walk in Fourier wavenumber space. This will be followed with some basic structure theorems which are useful for the construction of majorizing kernels.
The first result illustrates the use of a majorizing kernel h of exponent 1, to obtain existence and uniqueness of the solution given by the stochastic cascade representation for small enough initial data and forcing on a time interval that is solely constrained by the length of time for which the forcing remains small. Specifically one has the following theorem.
In the statements of these results, (−∆) −β , 0 < β < 3/2 denotes the negative power of the Laplacian defined as the integral operator with symbol |ξ| −2β ; i.e. the convolution with a Riesz potential.
As already noted the existence question is made trivial by the stochastic iteration so long as the indicated expected values are defined. A proof of uniqueness can be based on a discrete parameter martingale construction by truncated geneology adapted from LeJan and Sznitman [37] which is used to obtain a bound on the difference between two solutions in this space by 2P (T > n) for all n = 0, 1, 2, . . ., where T denotes the time to extinction of the underlying discrete parameter critical binary branching process.
It should be remarked that some regularity properties of the solutions can be inferred from the particular majorizing kernel being used via Payley-Wiener theory, Sobolev embedding, etc. For example, note that the majorizing kernel h 0 (ξ) = π 3 /|ξ| 2 gives existence and uniqueness, but no control over tracking regularity properties of the solutions from those of initial data. However, solutions obtained using the majorizing kernels h (α) β = |ξ| β−2 e −α|ξ| β , 0 < β ≤ 1, α > 0, maintain the C ∞ − regularity of the initial data. as can be seen from the bound on the Fourier transform of the solution Moreover one may construct, for example, smooth compactly supported initial data whose evolution is more efficientlly dominated by a kernel h (α) β with β < 1, whose decay at infinity is consistent with the compactness of support property; see Bhattacharya and Rao [5, pp. 83-89].
On the other hand, as will be elaborated below, working in the function spaces F h,1,T it is possible to use majorizing kernels to obtain spatial analyticity of the solution. However, it should be noted that the size constraints imposed on the initial data and forcing to obtain analyticity are substantially more restrictive than those required in Theorem 3.1; i.e. existence and uniqueness of such solutions is within a smaller ball of the function space. Specifically one has Theorem 3.2. Let h(ξ) be a standard majorizing kernel with exponent θ = 1.
ρν 2 e −1/2ν for some 0 ≤ ρ < 1. Then there is a unique solution u in the ball B 1 (0, R) centered at 0 of radius Under the conditions of Theorem 3.2 the asserted solution satisfies the following decay condition Thus Theorem 3.2 provides another approach generalizing that of Lemarié-Rieusset [35], Foias and Temam [20] to obtain conditions for regularity in the stronger form of spatial analyticity. In particular, by means of the following key inequality along the lines of Foias and Temam [20] and Lemarié-Rieusset [35], the majorization ofû(ξ, t) can be shifted to e √ t| |ξ| |û (ξ, t); see Bhattacharya et al. [6,Proposition 4.1,p. 5031]. As a consequence, for example, if e −δ|ξ| h(ξ) ∈ L 1 for some δ then, for times √ t > δ, u(x + iy, t) is complex analytic in the region is a complex analytic extension since for |y| ≤ √ t − δ. In particular, in the case e −δ|ξ| h(ξ) ∈ L 1 for all δ > 0, analyticity holds for all times t > 0. This is the case for h(ξ) = 1/|ξ| 2 treated by Lemarié-Rieusset [35], as well as for such majorizing kernels larger at infinity as indicated in examples below; c.f. Example (viii).
Two additional corollaries to the theory for global solutions follow immediately from the stochastic cascade representation: (a) time-asymptotic solutions, c.f. Bhattacharya et al. [7], Cannone [9] and references therein, and (b) self-similar solutions, c.f. LeJan and Sznitman [37], Cannone [9] and references therein. Specifically, for the case of (a), under the majorizing conditions of Theorem 3.1 with T = ∞, suppose further that lim t→∞ĝ (ξ, t) =ĝ ∞ (ξ) exists for each ξ = 0. Then, since the underlying discrete parameter binary branching is critical, lim t→∞ (θ, t) exists a.s. as a finite random product. Moreover, under the majorizing conditions one has for each t ≥ 0, with probability one | (θ, t)| ≤ 1. Thus, by Lebesgue's Dominated Convergence Theorem exists for each ξ. Now again apply Dominated Convergence to (F N S) h to obtain (69) Multiplication through by h(ξ) establishes that exists and satisfies the steady state Navier-Stokeŝ With regard to (b), mild self-similar solutions to Navier-Stokes of the form were originally considered by Leray [36] in connection with backward-in-time evolution of Navier-Stokes as an unrealized possible approach to determine singular flow. Cannone [9] provides a brief historical overview and ultimate demise of backward self-similarity, but the existence of forward self-similarity. In this regard the simplicity of the branching random walk cascade is overwhelming. Specifically, one observes that for majorizing kernels h defined by equality in (62) satisfying the limit-space scale invariance it follows that (i.) m(ξ) is constant and therefore trivially invariant under scale change ξ → ξ/λ, λ > 0; (ii)a ⊗ ξ b depends on ξ only through e ξ ≡ ξ |ξ| ; under the scale change ξ v → ξ v /λ, λ > 0, the distribution of the holding time S v changes to that of λ 2 S v ; and, finally, the spatial transition probability h(η)h(ξ−η) h * h(ξ) dη is invariant under re-scaling of Fourier frequencies. Thus taking initial dataû(ξ)/h(ξ) and forcingĝ(ξ)/h(ξ) (constant in time) to both be homogeneous of degree 0 one obtains thatû (ξ, t) = λ −2û (λ −1 ξ, λ 2 t).
Another variation on the general approach presented here involving semi-Markov branching, in the sense that the holding time distributions are not required to be exponentially distributed, leads to conditions for a local existence and uniqueness theory in all dimensions.; see Bhattacharya et al. [7], Orum [45], Ossiander [47].
Some sense of the class of majorizing kernels may be derived by noting from Hölder's inequality that the set of fully supported majorizing kernels with a given exponent is a log-convex set. Also if h(ξ) is a majorizing kernel then so is ce a·ξ h(ξ) for arbitrary fixed vector a and positive scalar c. One may also show that the only fully supported homogeneous majorizing kernels in n ≥ 3 dimensions are those of degree n − 1. Along these lines the Theorem 3.4 below may be viewed as something of a "toolkit" for constructing examples.
6. If h ∈ H n,θ , then e a·ξ h(ξ) ∈ H n,θ for any fixed a ∈ R n , 7. If h ∈ H n,θ then for any pseudo-metric ρ on a subset of R 3 containing W h , e −aρ(ξ0,ξ) h(ξ) ∈ H n,θ for any a > 0 and ξ 0 fixed. In particular e −a|ξ| β h(ξ) ∈ H n,θ for any a > 0 and 0 < β ≤ 1 In general, existence and classification of majorizing kernels is non-trivial. For example, it can be shown that any piecewise continuous h ∈ H 1,1 must have W h = (0, ∞) or W h = (−∞, 0). This illustrates the tradeoff between n and θ; if exponent θ > 0, the existence of majorizing kernels with support R n \{0} is problematic for smaller values of n. There are however fully supported majorizing kernels of exponent θ = 0 for all n ≥ 1. Using the toolkit Theorem 3.4 one may identify a number of majorizing kernels; see Bhattacharya et al. [6].
The kernels in Example (vii) are closely related to the Bessel kernels of Aronszajn and Smith [2]. They can also be combined with other kernels of the Example (vii) along the lines of the previous proposition to construct kernels in H n,θ for 0 < θ < 1. One may apply the Laplace method for estimating integrals to show that the Bessel type kernels h = h n,β,γ are also regularizing kernels in the sense that the distributions in the corresponding function space F h,0,T are C ∞ − functions.
In view of the role of majorizing kernels in providing bounds on the Fourier transformed forcing and/or initial data, the theory contains a dual problem which is to identify classes of divergence free vector fields in physical space which are so dominated. This section is concluded with some examples in this category.

Some Concluding Remarks
We wish to conclude with pointers to two major problem areas pertaining to incompressible Navier-Stokes equations which also fall within the subject area of probability theory.
Probabilistic methods in the study of Navier-Stokes equations are not new, howeverr new ideas have emerged recently which are based on objects of broad probabilistic interest such as multi-type branching random walks. The theory to date is in a form closely tied to the contraction map methods of Kato [30], so much so that most can be obtained by either route. Bounding the stochastic cascade product by one is a bit simpler than the contraction argument only in the case that the contraction is not strict. Similarly, as illustrated in this survey, owing to the a.s. finiteness of critical branching, uniqueness, time-asymptotic limits and certain self-similarity structure are extremely simple observations within the stochastic cascade framework.
The iterates of the Picard scheme can be probabilistically identified as expectations over truncations ( by generation) of the branching random walk. In particular, the n−th iterate is represented by an expected value restricted to averages over trees whose evolution by time t consists of only vertices v such that |v| ≤ n and such that the n-th generation vertices |v| = n survive beyond time t, i.e. n j=0 S v|j > t; see Proposition 4.3 of Bhattacharya et al. [6]. This and the uniqueness proof are two instances in which minimal but significant consideration is given to this underlying branching evolution.
However suitable exploitation of the genealogy of the branching random walk in ways which take fuller advantage of incompressibility has the potential to yield significant new results. Problems of this latter type are indeed at the next frontier of this research.
The Fourier representations are naturally constrained to free space and/or periodic boundary conditions. A second major problem for probabilists is the continued quest for representations in physical space which are amenable to the imposition of boundary conditions of various types.