A Renormalization Group Approach to Spontaneous Stochasticity

We develop a theoretical approach to ``spontaneous stochasticity'' in classical dynamical systems that are nearly singular and weakly perturbed by noise. This phenomenon is associated to a breakdown in uniqueness of solutions for fixed initial data and underlies many fundamental effects of turbulence (unpredictability, anomalous dissipation, enhanced mixing). Based upon analogy with statistical-mechanical critical points at zero temperature, we elaborate a renormalization group (RG) theory that determines the universal statistics obtained for sufficiently long times after the precise initial data are ``forgotten''. We apply our RG method to solve exactly the ``minimal model'' of spontaneous stochasticity given by a 1D singular ODE. Generalizing prior results for the infinite-Reynolds limit of our model, we obtain the RG fixed points that characterize the spontaneous statistics in the near-singular, weak-noise limit, determine the exact domain of attraction of each fixed point, and derive the universal approach to the fixed points as a singular large-deviations scaling, distinct from that obtained by the standard saddle-point approximation to stochastic path-integrals in the zero-noise limit. We present also numerical simulation results that verify our analytical predictions, propose possible experimental realizations of the ``minimal model'', and discuss more generally current empirical evidence for ubiquitous spontaneous stochasticity in Nature. Our RG method can be applied to more complex, realistic systems and some future applications are briefly outlined.


I. INTRODUCTION
Classical dynamical systems (ODE's and PDE's) with smooth vector fields v have unique solutions of the Cauchy initial-value problem for any fixed initial data x 0 . However, understanding has slowly consolidated that unicity is effectively broken when the field v becomes nearly singular and then small random perturbations of the deterministic dynamics can persist in the joint limit of vanishing noise and increasing singularity. In that case, deterministic trajectories become non-unique and solutions remain stochastic. This phenomenon occurs for Lagrangian fluid particles in high Reynolds-number turbulence [1,2], where it has been termed spontaneous stochasticity [3]. Although the effect was implicit in the work of L. F. Richardson on two-particle dispersion [4], it was only realized much later that classical determinism of Lagrangian trajectories breaks down and provides the mechanism for anomalous dissipation of scalars advected (passively or actively) by a turbulent flow [5][6][7][8]. Similar Lagrangian spontaneous stochasticity of magnetic fieldlines in turbulent astrophysical plasmas breaks the fluxfreezing constraint at high electrical conductivity, necessary to explain fast magnetic reconnection [9][10][11]. Entropy cascade in nearly collisionless plasmas likewise requires spontaneous stochasticity of Hamiltonian chargedparticle motions in the self-consistent electromagnetic fields of the limiting Vlasov kinetic equation [12,13]. Even before these Lagrangian manifestations were clearly understood, it was pointed out by E. N. Lorenz in a pioneering work [14] that Eulerian solutions of fluid equations can exhibit similar indeterminacy for multiscale turbulent flow. This Eulerian spontaneous stochasticity must be directly connected with non-uniqueness of dissipative weak solutions to the Cauchy problem for ideal fluid equations obtained in the high Reynolds limit or the "Nash non-rigidity" phenomenon [15][16][17]. There are profound scientific implications for predictability of weather and climate [18,19], but more generally for practical indeterminism of diverse problems in geophysics, astrophysics, and cosmology [20][21][22]. Despite the significant implications, there are few general theoretical tools available to predict when spontaneous stochasticity will occur and to determine what the observed statistics will be. We develop here a novel approach based on the renormalization group method from statistical physics and field theory [23][24][25]. Although the method can be applied much more generally, we focus on the "minimal model" that exhibits the basic phenomenon, the one-dimensional ODĖ with Hölder exponent/roughness exponent h.
Here v * (x) may be regarded as a caricature of a turbulent velocity field in the limit Re → ∞ which, according to Onsager's theorem [26,27], must develop singularities h ≤ 1/3 in order to account for anomalous dissipation of kinetic energy. Although the model velocity v * (x) is singular at only a single point (x = 0) when 0 < h < 1, it provides the standard textbook example of breakdown of uniqueness (e.g. see [28], p.2). The non-unique solutions, labelled by α = (±, τ ), have the form of delayed power-laws, with δ = 1 − h: x α (t) = 0 0 ≤ t ≤ τ, ±(A δ (t − τ )) 1/δ t > τ, (I. 3) where ± is the direction and τ ∈ [0, ∞] is a waiting time at the origin. The set of non-unique solutions is here uncountably infinite and, in fact, this is generally true for any bounded, continuous v by Kneser's theorem ( [28], Theorem 4.1). The above simple one-dimensional ODE and various generalizations have been extensively studied in the probability theory literature when supplemented with a white-noise random perturbation: See [29][30][31][32]. This initial-value problem is stochastically well-posed and has a unique (random) solution. In the limit of vanishing noise D → 0, it has been shown that the solution remains stochastic, with equal probabilities of 1/2 for the two extremal solutions with no initial delay (τ = 0) to occur: x ± (t) := ±(A δ t) 1/δ , (I. 5) and all other solutions of the limiting deterministic ODE have zero probability to occur [29,31]. The precise rate of vanishing of the probability to observe non-extremal solutions has also been derived as a "singular large deviations" [30], distinct from Freidlin-Wentzell predictions for smooth dynamical systems [33]. We wish to reconsider this model from the perspective of statistical physics. To render the model in principle experimentally realizable, it is necessary to impose an IR cut-off length L and also a short-distance (UV) regularization scale , for example, v(x) = A x (x 2 + 2 ) δ/2 , |x| < L (I. 6) so that "rough" power-law scaling v(x) ∼ A|x| h occurs only in the range |x| L. In a physical turbulent fluid flow, the UV cut-off is provided by viscosity. The precise behavior for |x| > L is immaterial but one of our main findings is that the details of the cutoff at |x| < influence the asymptotic behavior in the limit L/ 1, and render generally inapplicable many of the specific results obtained in the probability theory literature for the idealized mathematical problem with = 0. The statistics of the time-histories x(t) with added white-noise, as in (I.4), can be expressed in terms of path-integrals with the Onsager-Machlup action [34,35] ds |ξ(s) − v(ξ(s))| 2 (I. 7) e.g. for the transition probability density function: Dξ δ(x − ξ(t))e −S [ξ] .
(I. 8) The stochastic dynamics thus corresponds to a (0+1)dimensional Euclidean field theory, or a 1-dimensional equilibrium statistical system with the Onsager-Machlup action S[ξ] as "Hamiltonian". The joint limit L/ → ∞ and D → 0 is then quite closely analogous to a zerotemperature phase-transition in an equilibrium system such as the 1D Ising model, which appears in the joint thermodynamic limit L/a → ∞ and low-temperature limit T → 0 [36]. In this analogy, the non-unique solutions (I. 3) of the initial-value problem play the role of non-unique ground states of the Hamiltonian and the physics problem is to determine which "ground-states" appear in the limit and their precise probabilities. We shall apply a version of the renormalization group (RG) method to answer this question. To understand why RG is applicable to this problem, there are two key points. First, the physical mechanism of spontaneous stochasticity is always the forgetting of the initial data at sufficiently long times. This mechanism can be well-illustrated by the singular, deterministic ODE (I.2) whose exact solution for general initial data x 0 = 0 is easily found for h < 1 to be When t |x 0 | δ /Aδ, then the general solution becomes indistinguishable from the extremal solutions (I. 5) and all memory of the initial data is lost except for the overall ± sign. At long times one expects a universal statistical state to be achieved independent of the precise initial data and their statistics. There is thus a close analogy with critical phenomena in equilibrium spin systems, in which a universal statistical behavior is expected at length scales much greater than the lattice spacing a, independent of the microscopic Hamiltonian within some broad universality class, which is determined by an IR stable RG fixed point and its domain of attraction [23][24][25]. The second key fact that makes RG relevant is that the singular dynamics (I.2) possesses an exact scale symmetry according to which the transformation x b (t) := b −1/δ x(bt) (I. 10) for any real number b > 0 maps a solution x(t) to another solution x b (t). This symmetry is broken by the regularization of type (I.6) at length-scales < and also by the Langevin noise proportional to √ D in the stochastically perturbed dynamics (I. 4). However, the symmetry is expected to be restored, at least statistically, in the joint limit → 0 and D → 0. Here we may note that the extremal solutions (I.5) are scale-invariant, x ± b (t) = x ± (t) for all b > 0, and are thus perfectly "ordered" states, analogous to the homogeneous ground states of the 1D Ising model with all spins up or down. The existence of such a scale-symmetry (or set of scale-symmetries) of the limiting singular dynamics is not required for spontaneous stochasticity, but it is a typical feature of the ideal equations that govern high Reynolds-number fluids and plasmas [19,37]. Based upon these two ingredients, we shall construct an RG framework that permits a systematic analysis of SS in the 1D model problem.
The fundamental result we obtain is that there is an RG fixed point that governs the "critical limit" → 0, FIG. 1 Phase diagram for the position variance Υ at the inertial timet = 1 as a function of the state variables X = ln(Re) and Y = ln(P e). The white (black) solid lines plot Y = X + c 1(2) exp(X/2) and Y = c 3(4) for c1 = 0.10, c3 = 6.74 (c2 = 0.94, c4 = 3.97), the theoretical phase boundaries for large X, Y from the RG analysis.
D → 0 of our model. Just as for closely analogous T = 0 critical phenomena such as quantum critical points [38,39] or the order-disorder transition in the 1D Ising model [36], the singular critical point in spontaneous stochasticity can never be strictly achieved but its drastic effects are felt for finite values of and D over a very wide region. This is illustrated by the phase diagram in X − Y -parameter plane, Figure 1, with X = ln(Re) and Y = ln(P e), where Re is the "Reynolds number" and P e the "Péclet number", dimensionless versions of 1/ and 1/D, respectively. (See section III A). As order parameter for spontaneous stochasticity we take the variance of a similarly non-dimensionalized positionx(t) at timet = 1: This variance is suitable as order parameter because deterministic behavior is implied by Υ = 0 while spontaneous stochasticity corresponds to a value Υ = Υ * > 0 determined by the RG fixed point and independent of both X and Y. The plot in Figure 1 clearly distinguishes three phases, a "deterministic phase" D for extremely large Y where Υ 0, a "noise-driven phase" N at small Y where Υ(Y ) strongly increases as Y decreases, and an extensive "spontaneously stochastic phase" (SS) at sufficiently large values of X, Y where Υ = Υ * . The values of Υ are obtained from a numerical simulation of the model (section IV) but the "phase boundaries" shown as solid lines are analytic results from an RG analysis of the approach to the fixed point (section III C). Although the boundaries between the "phases" are somewhat arbitrary, they are qualitatively well-defined and match well our analytical predictions. In any case, it is clear from Figure 1 that the critical behavior at X = Y = ∞ is attained already for a broad range of finite X, Y values.
Before we present the details of our RG analysis, we must emphasize that spontaneous stochasticity is qualitatively distinct from the so-called "butterfly effect" or "deterministic chaos" of differentiable dynamics [40,41]. This familiar behavior corresponds to Lyapunov exponential growth of errors, which is observed in our 1D toy for h = 1. In that case, the solution of equation (I.2) for general initial data is Regarding x 0 = 0 as a small perturbation away from 0, it is obvious that for h = 1 the initial error is exponentially magnified over time but never "forgotten." In particular, for any fixed time t one observes that x(t) → 0 as x 0 → 0, whereas from (I.9) for h < 1 one observes at any fixed time t that x(t) → x ± (t) = 0 as x 0 → 0. This fundamental distinction is what characterizes spontaneous stochasticity. Even though the 1D toy problem has a smooth velocity when regularized as in (I.6), we shall see that the exponential error growth (I.12) either does not occur at all or else is an extremely short-lived transient that rapidly crosses over to power-law growth as in (I.5). As we shall argue below, this is the generic situation in high Reynolds-number turbulent flows and in similar problems with near-singular dynamics in geophysics and astrophysics.

II. REVIEW OF CURRENT EVIDENCE FOR SPONTANEOUS STOCHASTICITY
Spontaneous stochasticity is conjectured or even demonstrated to underlie some of the most basic turbulence phenomena and should be ubiquitous in turbulent flows. A simple Lagrangian manifestation of spontaneous stochasticity is the universal Richardson t 3 -growth in mean-square width of plumes of smoke or of other aerosols advected by a turbulent flow [5], which many observations suggest to be independent of the details of the source and of the mass diffusivity of the aerosol ( [42], section 24.3). More recently, it has been argued [22] that Eulerian spontaneous stochasticity explains how universal statistics are attained in a finite time for turbulent mixing layers during their phase of non-stationary, selfsimilar "equilibrium" spread [43,44]. In addition to many hints and indirect pieces of evidence for spontaneous stochasticity, there are a few controlled, precision studies which provide direct corroboration of the phenomenon. It is useful to review these here in order to provide context and motivation for our subsequent analytical investigations of our model problem.
We first consider Lagrangian fluid particles in turbulent flows, where spontaneous stochasticity is necessary to explain enhanced turbulent mixing. In fact, anomalous scalar dissipation for Re 1, P e 1 strictly independent of molecular transport coefficients D, ν can occur only if Lagrangian particle trajectories are spontaneously stochastic (at least away from solid walls) [5,6,8]. Therefore, the accumulated evidence from experiments and simulations for a scalar dissipative anomaly [45] provides indirect support for spontaneous stochasticity. It is impossible from empirical observations, however, to ever conclusively rule out a very slow decrease of scalar dissipation for increasing Re, P e. It is therefore desirable to have direct evidence for spontaneous stochasticity. The most basic observable effect on displacements r(t) = x (t)−x(t) of pairs of advected particles at positions x(t), x (t) would be the "forgetting" of the initial displacement r(0), analogous to the "forgetting" of x 0 in (I.9) for our model problem in the limit t |x 0 | δ /Aδ. Unfortunately, this phenomenon has proved difficult to verify in laboratory experiments. The papers [46,47] describe an experiment in a turbulent flow at Re 4.4 × 10 4 , seeded with neutrally bouyant polystyrene particles which were illuminated by two crossed laser beams and tracked with high-speed cameras. Those experiments observed a "ballistic" separation of particles, with r.m.s. separation r (t) ∝ v (r(0))t, exhibiting a strong dependence upon r(0) through v (r(0)) = |v(r(0)) − v(0)| 2 . The analogous dependence is seen in our model problem at short times by Taylor-expanding (I.9) for Since the time t 0 = r(0)/v (r(0)) required to "forget" r(0) in a turbulent flow grows ∝ (r 2 (0)/ε) 1/3 (assuming the Kolmogorov "2/3-law" [48] withε the mean energy dissipation per mass) a plausible explanation of the observations of [46,47] is that initial separations r(0) were insufficiently small to observe "forgetting". In fact, difficulties in tracking particles in densely seeded flows prevented [46,47] from accessing separations smaller than about 43 Kolmogorov dissipation lengths.
Numerical simulations are not subject to such limitations on r(0). The study [49] computed two forced isotropic turbulent flows at Reynolds numbers Re 1.4 × 10 4 and 3.5 × 10 4 and investigated deterministic Lagrangian particles for initial separations ranging from 2 − 24 Kolmogorov lengths. Figure 2 of [49] showed a very clear "forgetting" effect, with mean-square displacement r 2 (t) transitioning from a "ballistic" regime ∝ (v (r(0)) 2 t 2 for t t 0 instead to a universal regime ∝εt 3 as predicted by Richardson [4], independent of r(0) for t t 0 A similar effect had been seen even earlier for stochastic Lagrangian particles satisfying a Langevin equation such as (III.1) but with velocity taken from a numerical simulation of forced, isotropic turbulence at Re 1.25 × 10 4 [10]. This study released all realizations of the simulated Brownian particles at the same point and averaged over that release location. For two different values of Schmidt number, Sc = 1 and Sc = 0.1 [50], Figure 1 of [10] showed that r 2 (t) crossed from a diffusive scaling ∝ Dt for t t D ∼ (D 3 /ε) 1/4 over to a universal Richardson law ∝εt 3 for t t D , independent of the magnitude of Sc. The noise level mea-sured by Sc is thus "forgotten" at sufficiently long times. Using the same simulation database, paper [8] extended the previous study by considering stochastic Lagrangian particles without averaging over the release location x 0 and computed the probability distributions P (x, t) for three values Sc = 0.1, 1, 10. It was verified for t equal to about one large-eddy turnover time that these distributions were nearly the same for all three Sc values.
Empirical evidence has emerged also for Eulerian spontaneous stochasticity of the type first proposed by Lorenz [14] for multiscale turbulent systems. This effect must occur constantly in laboratory turbulence, as obvious from the apparent non-reproducibility of the specific details of such flows. It would be very challenging, however, to perform controlled laboratory studies because of the difficulty in creating an ensemble of initial conditions whose members differ in only fine details. Numerical experiments, on the other hand, can be readily performed and the most important scientific implications are in fact for the unpredictability and extreme ill-conditionedness of simulations for such singular turbulent flows. Numerical studies of predictability of threedimensional isotropic turbulence [51,52] have considered the mean-square "pair-separation" v(t) − v (t) 2 in energy-norm · for two Navier-Stokes solutions v(t), v (t) with initial data slightly different at small scales. These works verify the prediction [53] that that meansquare separation should grow ∝εt for long enough times, but unfortunately the crucial "forgetting" of the initial separation v(0) − v (0) was not investigated.
More direct evidence for Eulerian spontaneous stochasticity was found in [18] which studied a solution of the surface quasi-geostrophic (SQG) model for turbulent flow over a fractal topography. It should be noted here that weak solutions of the inviscid SQG equations expected to describe the infinite-Re limit are known to be nonunique for exactly prescribed initial data [54]. Figure 5 of [18] presented data on the r.m.s. separation ψ(t) − ψ (t) 2 of the stream-functions ψ(t) for the base flow and ψ (t) for a flow with initial condition perturbed by three sets of "error fields" of decreasing magnitude and length-scale. This separation was found to be almost unchanged as the initial errors decreased, consistent with Eulerian spontaneous stochasticity, or what [18] termed "the real butterfly effect". Even more compelling evidence comes from the paper [22] which studied the development of a turbulent mixing layer from a Kelvin-Helmholtz unstable singular vortex-sheet, where again solutions of the limiting inviscid Euler equations are known to be non-unique [55,56]. Perturbing the initial vortex-sheet with random "errors" of decreasing magnitude and simultaneously decreasing the viscosity ν, the crucial "forgetting" of initial error v(0) − v (0) was observed in Figure 1 of [22] for the mean-square separation v(t) − v (t) 2 after a sufficient time-interval. Furthermore, after this initial short transient, the spontaneously stochastic ensemble of solutions entered a selfsimilar growth phase with universal statistical properties independent of regularization and noise level, consistent with earlier observations [44].
We may mention finally that quantum spontaneous stochasticity in the semiclassical limit /m → 0 has been derived for a Hamiltonian version of our 1D model: when that is quantized in the canonical manner. See [57,58]. The second paper [58] also verified their analytical results for a regularized potential by numerical simulation of the Schrödinger equation and discussed some possible realizations by laboratory experiments. However, it would be extremely challenging in an experiment to verify quantum fluctuations as the origin of the stochasticity, precisely because the limiting spontaneous statistics are universal and do not depend upon the precise source of the random perturbations. Temperatures near absolute zero would be required to suppress thermal fluctuations and environmental noise sources would have to be very carefully controlled.

III. RENORMALIZATION GROUP THEORY
The experimental and numerical evidence presented in the previous section provide some solid confirmation that spontaneous stochasticity is generally present in high-Reynolds-number turbulence and that it underlies such basic phenomena as anomalous dissipation, enhanced mixing, unpredictability, and universal statistics in non-stationary, evolving flows. Presently, however, there are first-principles analytical derivations of spontaneous stochasticity only in some model problems, such as Lagrangian particle histories in the Kraichnan model [5,6,59,60] and in Burgers shock solutions [61], or semiclassical particle trajectories in some singular 1D quantum models [57,58]. The mathematical techniques employed in these works are special to the models considered and give no general physical intuition. We here develop an RG approach based upon a close analogy with zerotemperature critical phenomena, which provides not only physical insight into the nature of spontaneous stochasticity but also a powerful tool with which to derive the universal statistics, both analytically and numerically.

A. Dimensional Analysis of the Model
We first perform a careful dimensional analysis of the family of models we consider, which may all be expressed by the 1D Langevin equation where η(t) is white-noise with mean zero and covariance but with different choices of the drift velocity v(x). In all cases, v(x) is a "quasi-singular" function satisfying with some roughness exponent 0 < h < 1, while it is regularized for |x| < . In most of our analysis we assume that v (x) ≥ 0 and a simple, concrete regularization satisfying that monotonicity condition is the linear interpolation Any such regularization must have exactly one point |x eq | < where v(x eq ) = 0 and for most natural regularizations, such as the linear interpolation, x eq = 0. Although it is not necessary to our arguments, we assume that below for convenience of presentation. It is also convenient to impose the strict equality (III.3), although it suffices in fact to assume only that uniformly on compact sets, as in the example (I.6). To simplify the mathematics we shall analyze the model with no explicit IR cutoff L and simply restrict attention to times less than T = L δ /Aδ, before the extremal solutions reach |x| = L. In a realizable experiment, L would be set by the size of the apparatus but here we let x range over the real line. The details of the velocity field (or even its meaningfulness) for |x| > L do not matter for the effects we study and the only requirement is that L . These toy models may be interpreted physically in terms of a Brownian tracer particle in a 1D "turbulent" flow, where x represents particle position, t is time, D is molecular diffusivity, v(x) is fluid velocity, is viscous dissipation length, and L is a large-scale outer length. The exponent 0 < h < 1 gives the velocity Hölder exponent in the "inertial range" of the flow for < |x| < L, with h = 1/3 the Kolmogorov-Onsager value. The physical dimensions of the various quantities appearing are In defining dimensionless parameter groups it is useful to adopt the fluid-mechanical interpretation. We thus introduce a "kinematic viscosity" ν with dimensions (length) 2 /time via and an "integral-scale velocity" V = AL h , from which we construct the dimensionless parameters which we refer to as "Reynolds number" Re and "Péclet number" P e. The first number is a dimensionless measure of the singularity of the velocity, while the second number is a dimensionless measure of the weakness of noise. We may also introduce their ratio the "Schmidt number" Sc, which measures the importance of noise-weakness relative to velocity-singularity. In a physical fluid both molecular transport coefficients ν, D are thermodynamic functions only of temperature and pressure, and the second coefficient depends also upon molecular properties of the solute and solvent, such as their molar masses. Thus, in fluid turbulence experiments, Re and P e are usually made large by increasing V or especially L, while Sc remains nearly constant. Of course, in the context of our mathematical toy model, an arbitrary dependence of Sc upon Re may be assumed. Dimensionless variables can be introduced in a macroscopic or "inertial-range" description as and the resulting non-dimensionalized equations are The velocity field in these variables becomes Reynoldsdependent, or explicitly for the linear model: The problem to be solved in these variables is to consider the limit Re, P e → ∞ with random initial data selected from a distribution and with all other quantities held fixed. Another non-dimensionalization corresponds to a microscopic or "dissipation-range" description (III.14) where υ ν := A h is the "dissipation-range velocity". We take these variables as the default, without superscript, expressed in terms of which the dynamics becomes In this picture the velocity field is Re-independent except for the extent of its power-law "inertial range", as indicated explicitly for the linear regularization: The two sets of dimensionless variables are related by implying that (III.18) The limit Re → ∞ with inertial-range variablesx,t fixed thus becomes in dissipation-range variables a long-time, large-distance limit x → ∞, t → ∞ for some specified probability distribution P 0 (x 0 ) of initial data. It should be noted here that the extremal solutions (I.5) are invariant under the change of variables (III.18), so that x ± (t) = x ± (t).

B. The RG Flow and Fixed Points
To investigate the long-time, large-distance behavior of the solutions x(t) of the dissipation-range stochastic dynamics (III.15) we shall appeal to the analogy with critical phenomena discussed in the Introduction. Recall that in near-critical equilibrium lattice systems with spacing a there is scaling behavior at length-scales r in the range a r ξ = the correlation length, which is universal within a broad class of microscopic Hamiltonians H. This universal behavior is described by effective Hamiltonians H b of Kadanoff block-spin variables over b lattice sites, obtained by integrating out the microscopic degrees of freedom between scales a and ba. When rescaled from lattice constant ba back to a, these effective Hamiltonians approach an RG fixed point H * for b 1 [23][24][25]. In classical dynamics with spontaneous stochasticity there is expected to be a similar universality of the long-time statistics, which should be independent of the specific initial data and stochastic perturbations of the dynamics, within broad universality classes. This issue can be studied with path-integrals such as (I.8). Because the early-time statistics and dynamics are supposed to be irrelevant, one can integrate out the prehistories before some arbitrary time-window t w , altering the "bare" Onsager-Machlup action S[ξ] to some new window-dependent effective action Γ tw [ξ]. Varying the window size as bt w and then rescaling back to t w yields a b-dependent effective-action Γ tw,b [ξ]. This effective action can be expected to flow for b 1 to a scale-invariant fixed-point Γ tw, * [ξ] which will describe the spontaneous statistics within some broad universality class.
Carrying out this program for general classes of stochastic perturbations will require methods of functional renormalization group [62,63]. However, for the specific white-in-time perturbations in the Langevin equation (III.15) the stochastic dynamics is Markovian and the single-time probability density function P (x, t) satisfies a closed partial differential equation, the Fokker-Planck equation [64] with Markov operatorL: This permits a more elementary approach, because "integrating out" the early times t < t w corresponds simply to solving this equation with initial distribution P 0 (x) over the time interval 0 < t < t w . By varying the timewindow as bt w , the goal in this simplified RG approach is to find new effective initial data and effective stochastic dynamics for b 1. Analogous to "block-spins" on a rescaled lattice with fixed spacing a, we may define the random "effective initial data" at fixed time t = t w as (III.20) The multiplicative field-renormalization' b −1/δ is here motivated by the scaling symmetry (I.10) of the singular deterministic dynamics and we shall see that it guarantees that fixed point distributions of the "block-spin" variables x where P (x, t) is the solution at time t of the Fokker-Planck equation (III. 19). This formulation of the RG method is very similar to that used previously to study self-similar solutions, fronts and blow-ups of nonlinear PDE's [25,[65][66][67]. Using (III.21) and the formula (III.18), one can express the probability density of solutionx Re (t) of the inertial-range stochastic dynamics (III.11) as This key "bridging relation" directly connects the large-Re limiting behavior of the inertial-range dynamics to the large-b asymptotics of the RG flow.
This equation must be solved with the initial data Q 1 (x w , t w ) = P (x w , t w ). Although this equation is nonautonomous in b, it can be straightforwardly rewritten in autonomous form as This flow equation makes clear that the effective dynamics v b and effective noise ε b change under renormalization, as well as the effective initial data Q b . Note that the (v b , ε b ) sector of the RG flow is decoupled from the Q b sector. With our assumptions on the regularized velocity v(x) The limiting behavior of Q b is also easy to infer from its RG flow equation. Since the effective noise vanishes as b → ∞, the probability Q b in the limit is conserved along the characteristics of the limiting hyperbolic PDE.
In fact, the RG flow equation for Q b is equivalent to a Langevin equation for the random realization It is thus clear that as τ → ∞ the probability must be concentrated on the stable equilibria x eq w of the modified flow, satisfying lim τ →∞ w(x eq w , τ ) = 0, and in our model there are exactly three such points 29) The origin is always strongly linearly unstable for sufficiently large RG time τ since and under our assumptions v (0) > 0. Hence, any stable RG fixed point Q * will satisfy Q * (0) = 0. One might worry that an initial Q 1 (x w ) with a part ∝ δ(x w ) will have that part conserved in time under the RG flow. However, recall that Q 1 (x w ) = P (x w , t w ) and, even if the initial distribution P 0 (x 0 ) contained a piece ∝ δ(x 0 ), it is diffused away by integrating over any small window-time t w > 0. We give a careful analytic proof that Q * (0) = 0 in the following section III C, where we consider the rate of approach to the RG fixed points. The other two equilibrium points x ± w are the stable attractors of the RG drift and, in fact, independent of τ. All of the probability Q b thus flows as b = e τ → ∞ into the two equilibria x ± w . Their large-τ stability is due, obviously, to the additional linear term −x w /δ that appears in the RG drift (III.24).
The conclusion of this analysis is that the stable RG fixed points must have the form for some p in the range 0 ≤ p ≤ 1 or, equivalently, in terms of random realizations In section III D we shall completely characterize the domains of attractions of these RG fixed points, which are non-empty for each p and govern the long-time statistics of suitable initial distributions P 0 (x 0 ), velocity regularizations v(x) and Schmidt numbers Sc. This result is trivial for Sc = ∞, since in that case the integrals are invariants of the RG flow. Such continuous lines of RG fixed points occur in some equilibrium statistical mechanics models, such as the eight-vertex model, and are associated there with a breakdown of universality, in the sense that critical scaling exponents become continuous functions of the interaction parameters in the microscopic Hamiltonian [68,69]. There is a comparable breakdown of universality in the present case since the probabilities p, 1 − p asymptotically assigned to x ± w depend upon fine details of P 0 (x 0 ), v(x) and Sc (see section III D). This non-universality is associated with the lack of strong "chaotic" properties for the limiting deterministic dynamics, where every orbit has zero Lyapunov exponent except for the unstable fixed point at x = 0 which has infinite Lyapunov exponent. There is robust universality only for the case of symmetric initial distributions satisfying P 0 (−x 0 ) = P 0 (x 0 ) and symmetric regularized velocities satisfying v(−x) = −v(x). These symmetry properties are preserved by the stochastic dynamics (III.15) and by the RG flow (III.26) so that, within this class, there remains only the symmetric fixed point The convergence of the RG flow to a stable fixed point Q (p) * yields the limiting statistics of time-histories. Inverting the definition (III.20) of the "effective initial data" for b = t/t w gives and thus (III.33) implies that as t → ∞ (III.37) We can then obtain also the Re → ∞ limit ofx Re (t) in the inertial-range formulation of the problem by means of the scaling relation (III.18), which yields as Re → ∞ x Re (t) ∼x ± (t) = ±(δt) 1/δ with probabilities p,1 − p.
(III. 38) In particular, all of these conclusions hold with p = 1/2 if P 0 (x 0 ) = δ(x 0 ), i.e. if x(0) = 0 with probability 1. It follows that these models exhibit spontaneous stochasticity of a very simple type, with exactly two random "ground states" selected from the continuum of solutions (I.3) of the singular, deterministic ODE. Just as for the idealized model (I.4) previously considered in the probability theory literature [29][30][31][32], which in our language has vanishing Schmidt number Sc = 0, the solutions selected in the zero-noise limit are the two extremal solutions. It is worth mentioning, however, that this result is regularization-dependent and relies upon our assumption that v (x) ≥ 0. If instead the regularization of the velocity makes the origin very strongly stable, then it is possible that the equilibrium x o w = 0 of the RG drift persists with positive probability in the limit τ → ∞.
(We thank A. Mailybaev for this observation.) Because it is instructive, we show in Appendix A explicitly how the origin may be stabilized by such an "unnatural" regularization. In that case, the spontaneously stochastic limit can include with positive probability the identically zero solution x ∞ (t) ≡ 0 (infinite delay) in addition to the extremal solutions x ± (t) with zero delay.

C. Approach to the Fixed Point and Singular Large Deviations
The preceding discussion has identified the limiting statistical distribution of solutions of the singular ODE (I.2) obtained in the "critical limit" Re → ∞, associated to an RG fixed point. A general question of interest is corrections to critical scaling which may be obtained from the RG flow near the fixed point, e.g. the low-temperature behavior of the 1D Ising model near its analogous T = 0 critical point [36]. The limiting statistics in the 1D models we consider have great simplicity, because all histories other than the two extremal solutions having vanishing probabilities in the limit and become rare events. Thus, the finite Re-corrections that describe these rare fluctuations have the form of probabilistic large-deviations results. This issue was first studied for the idealized Sc = 0 version of the model in [30], which derived "singular large deviations" estimates that differ in scaling from those in the Freidlin-Wentzell theory. We here derive similar results for our more general models using the RG approach.

Deductions from the RG Flow
The relation (III.21) connects large-b asymptotics of Q b (x w , t w ) with large-x, large-t asymptotics of the Fokker-Planck solution P (x, t), where exponential decay is generally expected. We shall rigorously justify this exponential decay analytically in the following section, but here we simply make the ansatz of exponential decay at leading-order and then derive consequences from the RG flow (III.23). We expect that φ(x w , t w ) > 0 for |x w | = (δt w ) 1/δ and φ(x w , t w ) = 0 for |x w | = (δt w ) 1/δ , to be consistent with the limiting fixed-point behavior (III.32),(III.33). It is also plausible that This equation is obviously consistent with the conjec- with C = ln φ(0, t w ) an unknown constant, or, equivalently, This functional form is quite general, independent of Sc, the velocity regularization v(x) (so long as Q * (0, t w ) = 0), and the initial distribution P 0 (x 0 ). We shall see in section III C 2 that the value of prefactor 0 does depend upon some of these details. These large-b asymptotics are directly transformed to large-Re asymptotics in the inertial-range description by means of the "bridging relation" (III.22). Substituting the ansatz (III.39), one obtains a local large-deviations result where the second line uses (III.43). It should be noted that the equation (III.41) determining φ(x w , t w ) is a "quasi-steady" Hamilton-Jacobi equation: with a zero-noise Freidlin-Wentzell Hamiltonian [33] for the RG drift (III.24): (III.47) By combining (III.46) and the first line of (III.45) it is straightforward to see that the actionΦ likewise satisfies a time-dependent Hamilton-Jacobi equation with zero-noise Freidlin-Wentzell Hamiltonian for the singular velocity v * (x) in the limit equation (I.2). The Re-scaling in (III.44) is "anomalous" with respect to standard Freidlin-Wentzell asymptotics [33] for the weak-noise limits of Langevin equations, like the inertialrange model (III.11). The standard predictions would follow from the corresponding path-integral if naively evaluated by the Laplace method for Re 1. This would lead one to expect instead that for Re 1, with a different Re-scaling than in (III.44) and with the single-time action function (III.53) Equivalently,Ψ would naively be expected to be the solution of the time-dependent Hamilton-Jacobi equation is Re-dependent and tends to the singular velocity v * (x) as Re → ∞, so that the limiting action has non-unique "ground-states" ξ * satisfying S * [ξ * ] = 0. Because 1−h 1+h < 1 the probabilities of rare events decay more slowly with Re than would be naively expected [70].
The non-standard large-deviations prediction (III.44) can, however, be valid only for |x| < (δt) 1//δ , since the explicit expression in the second line of (III.45) shows thatΦ(x,t) < 0 for |x| > (δt) 1//δ . Likewise, the basic exponential-decay ansatz (III.39) cannot be valid for |x w | > (δt w ) 1/δ because it leads to the conclusion that φ(x w , t w ) < 0. The assumption of a stretched-exponential decay with b → b σ leads to similar inconsistency for σ < 1+h 1−h , but the ansatz with σ = 1+h 1−h , or explicitly (III.57) The new dominant balance in (III.57) as b → ∞ is between the terms proportional to b 1+h 1−h so that ψ(x w , t w ) must satisfy the equation Note that this equation is consistent with the requirement that ψ(x w , t w ) is positive and increasing in |x w | for |x w | > (δt w ) 1/δ . Furthermore, note that the equation (III.58) determining the action function ψ is a "quasi-steady" Hamilton-Jacobi equation In contrast to (III.47), the noise contributes now a quadratic term ∝ p 2 w . The intuition is that effects of singularity are so strong for |x w | < (δt w ) 1/δ that external noise does not influence the probabilities of rare events, but the dynamics is smooth for |x w | > (δt w ) 1/δ and external noise reasserts its role in producing rare fluctuations. In fact, the modified ansatz (III.56) in conjunction with the "bridging relation" (III.22) implies a standard largedeviations result of the form (III.52) with the expected Re-scaling and with the specific action function It is furthermore straightforward using (III.59) to show that the actionΨ in (III.61) satisfies the time-dependent Hamilton-Jacobi equation (III.54) with the inertialrange Freidlin-Wentzell Hamiltonian (III.55). We therefore conclude that standard Freidlin-Wentzell largedeviations theory must hold for the space-time region |x| > (δt) 1//δ exterior to the extremal solutions. This conclusion is entirely plausible, because optimal histories that produce the rare events in that region must depart the origin even faster than do the extremal solutions, through the agency of external noise.
Our results from analysis of the RG flow equations are completely consistent with the "singular largedeviations" theory of [30] for the idealized special case with Sc = 0. That paper rigorously obtained large-deviations results corresponding to our anomalous/singular estimate (III.44) in the interior region and our standard/regular estimate (III.52) in the exterior region. The analysis of [30] showed further that the prefactor 0 in the large-deviations rate functionΦ in (III.45) is the ground-state energy eigenvalue of a certain Schrödinger operator. In the following section we show for our more physical models with Sc > 0 that the constant prefactors 0 are likewise ground-state energies of suitable Schrödinger operators, using a more elementary argument than that employed in [30] and recovering their results in the limit Sc → 0.

Derivation of the Exponential-Decay Ansatz
We here derive carefully the exponential-decay ansatz (III.39) for the distribution Q b (x w , t w ) of effective initial data in the interior region |x w | < (δt w ) 1/δ , which is linked by (III.21) to the decay of P b 1/δ x w , bt w for b → ∞. We make a direct analysis of the long-time, large-distance limit of the Fokker-Planck solution by using a standard method of orthogonal eigenfunction expansions. See [64], section 5.4. As discussed there, the transformation Note that U (x) ≥ 0 when v (x) ≥ 0 and that outside the regularization region (III.65) Setting Sc = 1, the potential (III.64) thus coincides with the singular potential U * (x) of [30] for |x| > 1 (up to a trivial change of notations [71]) but differs from their potential in the regularized region |x| < 1.
Because U (x) ≥ 0, the eigenvalues n , n = 0, 1, 2, ... of the 1-particle HamiltonianĤ are all positive, with corresponding eigenfunctionŝ HΨ n = n Ψ n (III. 66) related to the eigenfunctions ofL andL * (III.68) The transition probability density can be expanded in these eigenfunctions as See [64], section 5.4, for all of these results. From the expansion (III.69) one can see that the solution P (x, t) of the Fokker-Planck equation (III.19) for any initial data P 0 (x 0 ) is given by the series which converges uniformly on compact sets of x for all t > 0. This result suffices to give large-time asymptotics of P (x, t) for fixed x, but to study Q b (x w , t w ) for b → ∞ we need as well a result on the large-distance asymptotics of the eigenfunctions Φ n (x), as follows: We briefly sketch here the proof. Since the derivation is independent of n, we suppress that index in the following. Starting with the eigenvalue equation we substitute Ψ (x) = exp(S(x)) (III. 73) to obtain the equation (III.74) For |x| → ∞ the dominant balance is between the pair of terms involving S (x) and U (x) ∼ Sc 4 |x| 2h , yielding the asymptotic result (III.75) Only the negative sign is consistent with squareintegrability dx |Ψ (x)| 2 < ∞, so we conclude that However, in order to determine the asymptotics of Φ(x) from that of Ψ (x) via (III.68) the leading-order asymptotic formula in (III.76) is inadequate, since and the leading-order term is exactly cancelled. Thus, we must determine a correction: Substituting into (III.74) gives with U 1 (x) = h/2|x| δ from (III.65). The dominant balance here is between the term linear in S 1 (x) and the constant , so that (III.80) yielding as claimed that We thus obtain from the eigenfunction expansion (III.70) an exact series representation for the effective distribution and also an asymptotic formula for large-b: The series converges for all b uniformly on compact subsets of the interior region |x w | < (δt w ) 1/δ and the first term in the series (III.82) dominates as b → ∞ whenever Φ 0 , P 0 = 0. Recall that the ground state wavefunction Ψ 0 of the 1-particle Hamiltonian has no nodes and can always be chosen positive, and it has two more derivatives than the potential U (x) defined in (III.64) from the regularized velocity [72]. Furthermore, from (III.76), Ψ 0 decays rapidly for |x| → ∞. The corresponding eigen-functionΦ 0 ofL * inherits all of these properties of Ψ 0 via the formula (III.68) and, therefore, 0 < Φ 0 , P 0 < ∞ for all possible choices of P 0 . We conclude that as b → ∞ (III.83) in the region |x w | < (δt w ) 1/δ , thereby establishing the exponential ansatz (III.39) and also identifying 0 as the ground-state energy of the 1-particle Hamiltonian defined in (III.63)-(III.64) [73] 3. Discussion of the Results The conclusion of our argument is that the approach to fixed points Q (p) * (x w , t w ) is the same for all p, independent of P 0 and with identical functional form for all Sc and for all regularized velocities v(x) satisfying v (x) ≥ 0. However, the prefactor 0 depends upon Sc and v(x) through the potential U (x) in (III.64). It is quite reasonable that the rate of approach to the limiting statistics as Re → ∞ should be regularization-dependent in the interior region |x| < (δt) 1/δ , since the histories that contribute to such rare events must have lingered in the regularization region |x| <¯ , with a macroscopic timedelay near the originτ :=t − |x| δ δ (III. 84) before departing along the corresponding non-extremal solution trajectory defined in (I.3). This argument provides a physical interpretation of the action function (III.45) for singular large-deviations asΦ = 0τ , with rate function directly proportional to the time-delayτ . The random noise plays no direct role in such large-deviations, because these rare events are mediated by the non-unique deterministic solutions. The prior analysis in [30] for the idealized Sc = 0 model made also implicit regularizations at x = 0 [74]. It is therefore useful to discuss in more detail the relation of our work with that in [30] and to discuss possible regularization-dependence of their results. Since they took = 0 from the outset, the "dissipation range" nondimensionalization that we adopted in section III A is unavailable. The corresponding non-dimensionalization for the Sc = 0 model is in "diffusive range" variables based on the diffusive length-scale D = (D/A) 1 1+h and time-scale t D = 2 D /D, and the only dimensionless group is the "Péclet number" P e. The diffusive-range nondimensionalization of the Sc = 0 model yields a Langevin equation of the same form as our dissipation-range equation (III.15) with Sc → 1, while the inertial-range nondimensionalization remains the same as in our section III A and yields the Langevin equation (III.11) with v Re (x) →v * (x). The macroscopic weak-noise limit is now P e → ∞ and, in our language, [30] obtain largedeviations results in the parameter P e. To formally derive their results as a special case of ours, we must transform our formulas in "dissipation-range" x instead to "diffusion-range"x usingx = Sc 1 1+h x, since / D = Sc 1 1+h , and then take Sc → 0. With this transformation, one easily obtains from our 1-particle Hamiltonian in the idealized limit Sc → 0. Note that the Hamiltonian in the square bracket coincides with that in [30] (see footnote [46]). We thus obtain the relationship of our ground-state eigenvalue 0 = 0 (Sc) with the corresponding eigenvalue of [30],ˆ 0 , as for Sc → 0. As long as Sc 1 while P e = Re Sc 1 then our singular large-deviation estimates in Re yield the large-deviations results of [30] in P e to good accuracy and show how to realize those predictions in a physical manner. We shall not attempt here to make this Sc → 0 limit rigorous, but our formal argument does suggest that the results of [30] should hold for all natural regularizations v(x) satisfying v (x) ≥ 0 in the limit → 0.
The asymptotics of the opposite case Sc 1 are a bit more subtle. This is a singular limit in which the kinetic term of the 1-particle Hamiltonian (III.63) formally vanishes and one expects a sort of "boundary-layer" of thickness 1/ √ Sc. Making the change to the stretched variable y = √ Sc x in the dissipation-range Langevin model (III.15) yields the modified equatioṅ for the variable y fixed and γ := v (0) > 0. In the limit Sc → ∞, one thus obtains the Langevin model with unstable linear drift corresponding to an inverted parabolic potential − γ 2 y 2 (see [64], §5.5.2). This has the form of our starting model (I.4) for A → γ, h → 1, D → 1. Although we have shown that spontaneous stochasticity appears for any finite value of Sc, however large, this Sc = ∞ limiting equation has smooth exponent h = 1 and thus no spontaneous stochasticity. In fact, the transition probability is given explicitly for our model (I.4) with A → γ, h → 1 by [64], Eq.(5.70) as: a Gaussian unimodal distribution centered at the unique deterministic solution x = e γt x 0 for all times t, with position width ∼ (D/γ) 1/2 e γt also growing exponentially.
Applied to the linear model (III.87), the result (III.88) implies that for any initial distribution P 0 (y 0 ) with position spread ∼ 1 There is no contradiction with our claim of spontaneous stochasticity for Sc 1, because the approximation (III.87) breaks down as soon as positions |y| √ Sc are attained with appreciable probability and this occurs within a time of order t ∼ (1/2γ) ln(Sc).
Finally, we make a methodological remark about our RG analysis. In traditional RG calculations for spin systems, the asymptotic approach to a fixed point follows from the linearized RG flow near the fixed point [23][24][25]. Similar linear analysis does not suffice to derive the corresponding results for our problem, such as the singular large-deviations (III.44). In fact, the linearization of the autonomous RG flow (III.26) for small departures (Q b ,ṽ b ,ε b ) from any fixed point (Q * , v * , 0) has the form The crucial observation is that the first term vanishes in the linearized flow equation for the perturbationQ b , or is supported entirely in that region (at least when x o w = 0 is made unstable by the regularization). Thus, the linearized RG equations for (Q b ,ε b ) andṽ b are decoupled, and, in particular, the rate at whichQ b → 0 as b → ∞ must be independent ofṽ b := v b − v * according to the linearized equations. We have seen, however, that the asymptotic approach of Q b to any fixed point Q * (and even which fixed point is attained) depends upon the precise form of v b in the regularization region wherẽ v b = 0. The linearized RG flow near each fixed point is thus insufficient to describe the approach to that fixed point. This is a fundamental difference from RG analyses of equilibrium critical systems, where the RG flow is generally well-described by an autonomous ODE in a finite-dimensional space of "relevant variables". For our problem, however, the RG flow is intrinsically in a space of functions (Q b (x w ), v b (x w ), ε b ) and governed by the set of PDE's (III.26). In that case, even the question whether (Q b , v b , ε b ) is "close" to the fixed point (Q * , v * , 0) becomes subtle, because there are many inequivalent norms, and a global analysis of the RG flow equations is required to study the asymptotic approach.

D. Domains of Attraction
Each of the RG fixed points Q (p) * for p ∈ [0, 1] has its own domain of initial distributions P 0 that it attracts, which depends sensitively upon both velocity regularization v(x) and Schmidt number Sc. Fortunately, for this simple 1D model, we can completely characterize the domains of attraction a priori, as we show now.
Because the fixed points assign non-vanishing probabilities to only the two states x − w = −x + w , it is clearly sufficient to consider the probabilities to be right or left of the origin [75], which are defined as and to determine their long-time asymptotics for t t 0 . Using the the fact that the transition probability p(x , t |x, t) satisfies the backward Kolmogorov equation in its initial space-time arguments (see [64], section 4.7) and using p(x , t|x, t 0 ) = p(x , 0|x, t 0 − t) by timehomogeneity, it then follows that (III.95) Since the adjoint Fokker-Planck operatorL * has all of its spectrum in the left half of the complex plane, the long-time limits exist and are independent of t 0 satisfying the stationary backward equation where V (x) as before is any anti-derivative of v(x). For the velocity regularizations satisfying v (x) ≥ 0 that we consider, p 0 = lim t→∞ P (x(t) = 0) = 0, so that the normalization N is fixed by Note that p ± (x) are non-trivial zero modes of the backward Kolmogorov operatorL * , in addition to the trivial constant zero mode [76]. For any solution P (x, t) of the Fokker-Planck equation, therefore, the integrals are exact constants of the motion and give the asymptotic probabilities to be right or left of the origin. If Sc → ∞ then p ± (x) converge to the Heaviside step functions θ(±x), consistent with the fact that the initial val- The argument made above for the Fokker-Planck equation fails to carry over directly to the RG Fokker-Planck flow equation (III.23), because the latter is nonautonomous and also the integrals such as (III.100) diverge for the RG-drift (III.24). Similar zero modes and associated conservation laws do not apparently exist for the RG flow. It follows directly from definition (III.21) of the distribution of "block-spins" that but these integrals all change with b = t/t w . On the other hand, after taking the joint limit t, b → ∞ the above relation implies that (III.103) Therefore, the fixed-point Q (p) * that attracts any initial distribution P 0 is fully characterized by the constants of motion (III.101), with p = p + , 1 − p = p − .
The existence of these exact integrals of motion help to justify our RG strategy for this degenerate model, with its continuous line of fixed points. The RG method is most successful when statistics are universal, e.g. when the details of the microscopic Hamiltonian are irrelevant to the long-distance critical scaling and simplified coarsegrained models belong to the same universality class. Our strategy of considering "effective initial data" after some arbitrary time window t w > t 0 is premised on the idea that the early-time statistics and dynamics are similarly irrelevant and may be ignored. Although universality breaks down in this specfic 1D model, fortunately the true initial distribution P 0 (x 0 , t 0 ) and the "effective initial distribution" Q 1 (x w , t w ) = P (x w , t w ) have the same integrals of motion (III.101). One can thus infer that even though the latter integral is not an invariant quantity under the RG flow. We are therefore justified in ignoring early times t < t w , since the universality class of the initial distribution P 0 is fully encoded in Q 1 .
It is only for the completely symmetric case, with P 0 (−x 0 ) = P 0 (x 0 ) and v(−x) = −v(x), that there is robust universality in this simple 1D problem, of the same sort encountered in most short-range, equilibrium critical spin systems. For any symmetric velocity v(x) from (III.98),(III.99), so that any symmetric initial distribution P 0 (x 0 ) then has p = 1 − p = 1/2. The symmetric RG fixed-point Q * = Q (1/2) * in (III.35) governs the long-time dynamics for all initial distributions and velocity regularizations in this symmetric class. It is worth remarking, on the other hand, that the domain of attraction of the symmetric fixed point is larger than the class of symmetric data and for any Sc < ∞ there are nonsymmetric P 0 and v that flow into (Q * , v * ). We construct in Appendix B explicit examples of non-symmetric initial distributions P 0 with invariants p + = p − = 1/2, which thus belong to the basin of attraction of Q * .

IV. EMPIRICAL STUDY OF THE 1D MODEL
The rather exhaustive analytical investigation of our model problem (III.1) in the preceding section yields detailed predictions that can be verified by direct numerical simulations. We present such numerical results here, because they render concrete our rich set of mathematical results and also motivate some additional discussion of the theory. Furthermore, it is worthwhile to explore physical systems that could realize our 1D model in a controlled laboratory setting, even if such experiments do not correspond directly to any naturally occurring circumstances. The scarcity of direct experimental evidence for spontaneous stochasticity is vexing, especially considering its presumed ubiquity in Nature.

A. Possible Experimental Realizations
The simplest set-up to realize our 1D model (III.1) is a mechanical experiment using a bead sliding with strong friction on a bent thin wire, as pictured in Figure 2. The equation for arclength s(t) as a function of time t can be obtained by ignoring inertia ms and balancing friction force −γṡ and parallel component of gravitational force −mg dy/ds. Note that the maximum velocity obtained when the wire is bent vertically is v * = mg/γ and any velocity field satisfying |v(s)| ≤ v * can be achieved by choosing a wire in the shape y = y(x) given implicitly by the equations Figure 2 shows the shape required so that the equatioṅ is satisfied for a wire of length 2L and h = 1/3. In any physical experiment this shape will be regularized somehow for |s| < , where the length scale will be no smaller than the radial diameter of the wire. A Langevin whitenoise force η(t) can be applied to the bead, for example by using commercially available white-noise voltage generators [77,78] and a transducer device. In fact, one could add noise in many other ways (for example, randomly jiggling the wire) and, as long as the noise is small, it should not matter for the spontaneous stochasticity phenomenon. A practical alternative to experiments with a mechanical apparatus would be an electric-circuit realization of the system (IV.2), based on the standard impedance analogy with mechanics. Large friction in a mechanical experiment (or large resistance in an electrical circuit) yields the first-order system (I.2) only as a singular limit and, as well-known, there will generally be transient, early-time deviations ( [79], section 3.5). It may therefore be desirable to realize instead the second-order Hamiltonian equation (II.1) with A = mg/L h by the opposite limit of negligible friction. In this case, the equation for arclength s(t) as a function of time t is obtained by balancing inertia ms and and parallel component of gravitational force −mg dy/ds. Similar to before, the maximum acceleration obtained when the wire is bent vertically is g and any acceleration field satisfying |a(s)| ≤ g can be achieved by choosing a wire in the shape y = y(x) given by the equations (IV.1) withv(s) →â(s), whereâ(s) := a(s)/g. Because the equations are identical in form, the wire shape pictured in Figure 2 is also that required to realize the second-order model (II.1) with h = 1/3, in the limit of low friction. Similar mechanical models have, in fact, been invoked previously in the natural philosophy literature in discussions of causality and determinism for classical dynamics [80,81]. In particular, the equations (IV.1) can be explicitly integrated for exponent h = 1/2, as with |s| < L and the resulting curve has been called "Norton's dome" [80,81]. However, these works by philosophers and also those by mathematical probabilists [29][30][31][32] have not recognized the physical necessity in any real experiment of a regularization of the singularity in the dynamics (II.1) at some non-zero length-scale . In particular, uniqueness of solutions holds for any finite regularization and some external source of noise, such as a Langevin force √ 2D η(t) added to (II.1), is required to supply stochasticity, which persists in the joint limit D → 0, → 0. Although we shall not consider here further this stochastic perturbation of the Hamiltonian model (II.1), spontaneous stochasticity does occur for this model in the joint limit of vanishing noise and regularization just as for the quantum version studied in [58].

B. Numerical Simulations
In lieu of laboratory experiments, we present here the results of numerical simulations of our model problem, for Hölder exponent h = 1/3. We solve the Langevin model using the simplest Euler-Maruyama discretization [82], employing either the dissipation-range formulation (III.15) or the inertial-range formulation (III.11) for numerical integration, depending on the situation. Empirical averages were then constructed by averaging over N independent samples. Full details on numerical methods and convergence tests are described in Appendix C. We consider only the deterministic initial distribution P 0 (x 0 ) = δ(x 0 ), which lies in the domain of attraction of the symmetric RG fixed-point (III. 35).
FIG. 3 Probability density functions P (x, t) for particle positions at four different times, all quantities in dissipation-range scaling units. The asterisks indicate the positions of the two extremal solutions at the same four times.
In Figure 3 we plot the probability density functions P (x, t) for particle positions in dissipation-range scaling, obtained for Sc = 1 at times t = 0.5, 3, 5.5, 8, starting at x = 0 at time t = 0. The PDF's are close to Gaussian at early times, but then flatten at time t . = 0.95 and thereafter become bimodal. The PDF's very near this splitting time are shown in Figure 7 in Appendix C. The two peaks of the bimodal PDF's P (x, t) approach the extremal solutions x ± (t) of the singular ODE (I.2) as time advances. A plot of the local extrema (maxima and minima) of the PDF's in the x-t−plane, Figure 4, shows the clear signature of a pitchfork bifurcation, as should be expected for the symmetry x ↔ −x [79]. In our simple model with just two extremal solutions selected in the zero-noise limit, this pitchfork bifurcation can be regarded as the "onset" of spontaneous stochasticity.

FIG. 4
Numerical bifurcation diagram, with local extrema of the empirical PDF P (x, t) plotted versus time t, local maximima in blue and local minima in red. The light orange markers indicate the region within which the values of the PDF vary by less than 0.1% of its extremal value. The solid cyan line is a least mean square fit of a parabola to the outer prongs of the diagram, the good agreement being consistent with the normal form of a pitchfork bifurcation [79].
Further insight can be obtained by plotting the vari- ance x 2 (t) = dx x 2 P (x, t) versus time t, in log-log coordinates. For Sc = 1 we obtain the data plotted with blue dots in Figure 5, which shows as straight lines with distinct slopes the two power-laws for a crossover time t c (Sc) ∼ (1/Sc) 1−h 1+h , with Sc → 1. The same behavior is also observed for Sc 1 as illustrated by the results for Sc = 10 −4 plotted as green dots in Figure 5. The only effect of the increased noise is that the onset of spontaneous stochasticity is delayed. Consistent with our asymptotic analysis for Sc 1 in section III C 3, this onset time is alwayst c ∼ 1 when expressed in the diffusive scalingt = Sc that the variance grows at early times instead as ln Sc (IV.5) and then crosses over to the universal power-law ∼ (δt) 2/δ at longer times. This behavior for Sc 1 is verified with γ = 1 by the data for Sc = 10 4 plotted with red dots in Figure 5. Note that (IV.5) yields also a linearly-growing variance ∼ (2/Sc)t at very early times t 1/2γ but an exponential growth ∼ (1/γSc)e 2γt at the intermediate range of times 1/2γ t (1/2γ) ln Sc. Both of these growth laws are a short transient behavior before the universal power-law scaling ∼ (δt) 2/δ at long times, which is a key signature of spontaneous stochasticity.
With these results as backdrop and with the asymptotic analyses of section III C 3 for Sc 1 and Sc 1 at hand, we now can explain the phase-boundaries plotted in Figure 1 in the Introduction. We consider first the crossover between the "deterministic phase" D and the "spontaneously stochastic phase" SS. As is obvious from the "bridging relation" (III.22), inertial-range position variablesx(t) will remain random in the limit Re → ∞, P e → ∞ only if the bifurcation in the dissipation-range PDF P (x, t) occurs at a time less than or equal to t ∼ Re for an order unity constant c or, with P e = Sc Re, ln(P e) ln(Re) + c 1(2) exp(a ln(Re)), a = 1 − h 1 + h (IV.7) These are the boundary lines plotted in Figure 1 with the constant c 1(2) chosen to provide a least-squares fit to the Υ−isoline for the value 0.85Υ * (0.15Υ * ), with Υ * the fixed-point value. For the other crossover between the "noise-driven phase" N and the "spontaneously stochastic phase" SS, we note from the asymptotic analysis for Sc 1 in section III C 3 that spontaneous stochasticity will be observed as long as P e 1, no matter how small Sc may be. Thus, we obtain now as an approximate phase-boundary ln(P e) ln(P e c ) := c 3(4) (IV. 8) where P e c can be regarded as a "transitional Péclet number" for spontaneous stochasticity in the Re → ∞ limit. This corresponds to the second set of boundary lines plotted in Figure 1 with constants c 3(4) selected for a best fit to the Υ-isoline for the value 1.15Υ * (1.85Υ * ). Finally, by simulating the model (III.11) in inertialrange units for Sc = 1, we have calculated the corresponding position PDF'sP (x,t) at the fixed timet = 1 These PDF's are plotted in Figure 6 for our default Hölder exponent h = 1/3 at a sequence of increasing values of Re = P e. It is seen that the PDF's converge to a mixture of two delta-function spikes located at the extremal solutionsx ± (t), each with probability 1/2. This figure thus illustrates the most basic effect of spontaneous stochasticity, that the historiesx(t) remain random for Re, P e → ∞. As a "control experiment" we have performed the same analysis for our model with h = 1, where insteadx(t) converges to the deterministic limitx ∞ (t) ≡ 0. See Appendix D, Figure 8. Only for singular dynamics with h < 1 can solutions be spontaneously stochastic. Also predicted by our theory is the rate of vanishing of probabilities to observe particles at positionsx of non-extremal solutions in the limit Re, P e → ∞. To make a quantitative comparison, we have plotted in Figure 6 over the interior interval x − (t) <x <x + (t) the singular large deviations estimate (III.44). The constant prefactor of the rate function Φ(x,t) in Eq.(III.45) was obtained numerically using the MATLISE software package [83] to solve for the groundstate energy of the 1-particle Hamiltonian (III.63), which gave 0 . = 0.5545 (see Appendix C). The theoretical prediction clearly matches the empirical PDF's better as Re, P e increase.

C. What is Special, What is General?
The 1D equation (I.2) provides perhaps the simplest possible model of spontaneous stochasticity but it also suffers from a number of closely related special features that distinguish it from more generic cases expected to occur in Nature. These special properties must be taken into account when using the results of the previous section to form expectations about observational signatures of spontaneous stochasticity in general.
The first special feature of the model (I.2) is that only a discrete subset of "ground-states" is selected from the continuum in the limit of increasing singularity and vanishing noise; in fact, precisely two. Most histories therefore have vanishing probability and the approach to the limit has the form of a large-deviations result, but with the novel feature that there is more than one minimum of the action and the usual "law of large numbers" breaks down. Very similar behavior occurs for Lagrangian spontaneous stochasticity in shock solutions of Burgers equation, where only the two extremals are selected from all generalized characteristics at the shock and the probabilities for all other histories vanish exponentially [61]. More generally, however, one expects that a continuous infinity of "ground-state" solutions will be selected, with a non-trivial limiting distribution over them. Numerical studies such as [8,18,22] exhibit no tendency for empirical distributions to be supported on only a finite number of solutions. Nevertheless the limiting probability distributions are expected to be supported on (generalized) solutions of the limiting singular equations and there should quite generally be a large-deviations-type estimate on vanishing probability of non-solutions.
A second special feature of the model (I.2) is that the dynamical vector field v * (x) has a singularity at just one point, x = 0, and the corresponding Cauchy problem has a non-unique solution only if the initial datum is chosen at that singular point. This feature is obviously related to the first, since a space-filling set of singularities would allow the solutions to "branch" at every point and the zero-noise limit could select an infinite ensemble of solutions. C.f. the example in [28], section II.5. There may be some problems of physical relevance where singular points are isolated in the dynamical state space, e.g. multi-body collisions in the classical N -body problem [84]. For typical situations in fluid turbulence, however, the limiting dynamics is singular at essentially every point and non-uniqueness should be generic. For example, the fluid velocity vector field that governs the motion of Lagrangian particles for incompressible turbulence has as Re → ∞ a space-filling set of Hölder singularities with exponent lim p→0 ζ p /p := h * 1/3 (with ζ p the velocity structure-function exponent) [85,86]. Likewise, Eulerian spontaneous stochasticity is expected to be generic in high-Re incompressible fluid turbulence. As nicely reviewed in [18], smooth (strong) solutions v(t) of the incompressible Navier-Stokes equation are unique for fixed initial data v 0 and even a Leray (weak) solution v (t) with nearby initial condition v 0 satisfies the inequality (IV.9) with E(t; v) := ν ∇v(t) 2 the volume-integrated energy dissipation per mass for solution v. This inequality implies that v (t) → v(t) as v 0 → v 0 for fixed viscosity ν, but because of the turbulent "dissipative anomaly" the bound (IV.9) is expected to diverge as ν → 0 and the smooth solution v becomes more nearly singular. Although the inviscid ν → 0 limit is still very poorly understood mathematically, the arguments of [14,53] on turbulent unpredictability suggest that limiting dissipative Euler solutions should be typically non-unique and exhibit the "Nash non-rigidity" phenomenon [15][16][17].
Finally, a third special feature of the model (I.2) is that it is "almost non-chaotic" since the non-fixed orbits must have zero Lyapunov exponent, as for any continuous-time dynamical system, and there are only two such orbits together tracing the entire right half-space {x > 0} and left half-space {x < 0}. However, the spontaneous stochasticity phenomenon requires not only positive Lyapunov exponents as in standard deterministic chaos but in fact infinite Lyapunov exponents, since initial-data arbitrarily close to each other must separate to the same distance in finite time. In the model (I.2) the unstable fixed point x = 0 has the Lyapunov exponentλ = ∂v ∂x (0) ∼ Re generically for any natural regularization and thus all of the "chaotic" dynamics arises from x = 0 whereλ → ∞ as Re → ∞. On the other hand, fluid flows where spontaneous stochasticity is expected to arise as Re → ∞ are known to possess much more robust forms of both Lagrangian and Eulerian chaos, even at moderate Reynolds numbers [41,87]. Here it is important to stress that spontaneous stochasticity requires not only positive Lyapunov exponents but also increasingly singular dynamics in some limit, such as Re → ∞ for fluids. The implications of spontaneous stochasticity are also much stronger. Whereas standard chaos theory for smooth dynamical systems implies the universality of infinite timeaverages described by a natural invariant measure on a strange attractor in the weak-noise limit [41,88], spontaneous stochasticity produces universal statistics in a finite time for vanishingly small noise. In order to exhibit such finite-time universality, spontaneously stochastic solutions do not have to lie on a chaotic attractor. In fact, a steady-state attractor does not even have to exist, as exemplified by such common flows as a decaying turbulent wake or a growing turbulent mixing layer.

V. FUTURE DIRECTIONS
We have developed in this paper a novel renormalization group (RG) approach to spontaneous stochasticity and applied it to the Langevin model (I.2), which is probably the simplest system conceivable which exhibits the basic phenomena. However, the general RG theory can be applied as well to much more complex and realistic systems with a limiting singular dynamics exhibiting spontaneous stochasticity. We here briefly indicate some directions for future RG work, first for Lagrangian particle histories and related time-histories governed by finite-dimensional ODE's, then for Eulerian space-time histories governed by infinite-dimensional ODE's/PDE's.

A. Lagrangian Spontaneous Stochasticity
The simplest generalizations of the model (I.2) are multi-dimensional extensions with point-singularities as considered in [89,90], where the singular radial dynamics is similar to that in the 1D model, while the angular dynamics in hyperspherical coordinates is smooth but otherwise arbitrary. As discussed in detail in [89,90], these mathematical model problems already exhibit a rich variety of behaviors including robust spontaneous stochasticity. These models have also the same fundamental scaling symmetry (I.10) as the 1D model (I.2) for a specified Hölder exponent h. Thus, our RG theory carries over straightforwardly to this class of systems, although many new features now appear, such as continuous lines of RG limit cycles corresponding to spontaneous statistics with broken O(2) symmetry. The class of models developed in [89,90] provide important testbeds to understand better the role of chaotic dynamics in producing robust spontaneous stochasticity. As discussed already in section II, chaotic dynamics and positive Lyapunov exponents are necessary to produce spontaneous stochasticity, but are not by themselves sufficient. The singularities in the dynamics play a crucial role in amplifying infinitesimal noise to macroscopic scales in finite times.
A step closer to realistic turbulence, our RG method should provide a useful framework to study Lagrangian particle histories also in the Kraichnan model of turbulent advection [1,2,5,6]. Lagrangian spontaneous stochasticity was first discovered in this model and there are even rigorous demonstrations of its existence [59,60]. However, little is known still in the Kraichnan model about (almost sure) properties of spontaneously stochastic particle-history distributions for fixed velocity realizations [91]. The Kraichnan velocity is a fractional Brownian random field in space with specified Hölder exponent h, so that it enjoys a statistical scaling symmetry analogous to (I.10) and a version of our RG method should be applicable. In addition to the Kraichnan model which is white-noise (delta-correlated) in time, our RG method should apply also to particle advection by other self-similar Gaussian velocity fields with finite time-correlations [3]. Because all of these Gaussian advection models have velocity fields that are everywhere singular in space and time, similar to real fluid turbulence, one expects that continuous branching will lead to spontaneously stochastic ensembles supported on an uncountably infinite set of particle histories.
Real fluid turbulent flows and Lagrangian particles therein can be treated also by our RG theory. Although there are presently formidable difficulties to carrying out such a study analytically, the RG analysis can be implemented numerically following schemes similar to those used previously for front and self-similar solutions of deterministic PDE's [92]. A very interesting feature of incompressible fluid turbulence is that the Euler fluid equations expected to describe the Re → ∞ limit have infinitely-many scaling symmetries of the form for all real exponents h, each of which gives an action of the one-dimensional group of dilatations, and these symmetries are expected to be statistically realized by multifractality of the velocity field [37,93]. In that picture, each spacetime event (x, t) of the turbulent flow has its own "local Hölder exponent" h = h(x, t) which is determined dynamically. In the previous numerical study of spontaneously stochastic particle distributions in isotropic turbulence [8], it was in fact observed that an intermediate-time nearly self-similar regime was attained for individual events, without any need to average over (x, t). This regime can be further explored with a numerical RG study, where now the field-renormalization exponent analogous to 1/δ in (III.20) must be determined self-consistently for each (x, t) [92]. The distribution of "effective initial data" that leads to self-similar growth backward in time is presumably universal for each h-exponent, since those singularities are expected to be created by dominant self-similar instanton solutions [94].

B. Eulerian Spontaneous Stochasticity
Our RG approach can be applied also to study Eulerian spontaneous stochasticity. Very attractive cases for numerical RG analysis are the self-similar "equilibrium" states in turbulent shear-flows and gravitationallyunstable layers, where evidence has already emerged for robust spontaneous stochasticity [21,22,95]. In these flows the Hölder exponent h that governs the self-similar scaling is set by the singular initial data, so that no self-consistent determination of h is required. Numerical RG in this application will allow the long-time selfsimilar regime to be more accurately probed, alleviating the restriction on direct numerical simulations imposed by finite domain size. These canonical flows are simple enough that analytical RG analysis may also be possible. An interesting feature here is that both the singular initial data and also the limiting ideal Euler dynamics are P T -invariant, e.g. under the transformation u(x, y, t) → −u(x, −y, −t), v(x, y, t) → v(x, −y, −t) (V.2) for the 2D singular vortex sheet in [22]. However, the linear fluid instability breaks this symmetry [96] [97] and the spontaneous statistics associated to a growing turbulent layer are obviously not P T -invariant. Here the viscosity ν of the fluid plays the role of a symmetry-breaking field in a phase-transition, whose effect does not disappear in the inviscid, zero-noise limit.
An ultimate aim of the RG theory of spontaneous stochasticity must be to treat singular Euler solutions which describe turbulent inertial-ranges and which, for generic initial data, should experience the "inverse errorcascade" that was predicted by Lorenz [14] and verified by following work [18,20,[51][52][53]. This inverse cascade is characterized by an error field that grows self-similarly in scale [14,53] or, in logarithmic wavenumber and time variables, by a "stochastic front" that leaves universal spontaneous statistics in its wake [20]. RG analysis of this type of multiscale spontaneous stochasticity presumably requires not only windowing out early times but also eliminating high wave-number modes. The effective dynamics for long times and low wavenumbers at the RG fixed point should contain not only "eddy-viscosity" effects that account for dissipative anomaly but also "eddy noise" associated to the stochastic anomaly [98,99]. One suspects that some type of turbulent fluctuationdissipation relation will connect these two effects, which would be great practical interest in numerical modelling [19]. A very remarkable numerical observation in [20] is that the spontaneous statistics achieved in a finite time for individual turbulent velocity realizations is identical to that observed for infinite-time attractors in forced steady-states. This result underlines the fundamental role of spontaneous stochasticity in achieving universal statistics for general turbulent flows. It presents also an imposing challenge, because the anomalous scaling and multifractality that characterize the turbulence state have so far defied theoretical analysis. Functional renormalization group methods that are non-perturbative and formally exact [62,63] have already yielded some novel predictions for high-Re steady-states of Navier-Stokes turbulence [100,101] and these methods should be capable also to describe spontaneously stochastic solutions.
Work on many of these directions is currently in progress and will be reported in following publications. that the exit time T (γ) a,x for the general equation (A.6) can be related to that studied in [102] as and thus We can then estimate the probability that the particle exits the interval (−a, a) before time t by the Chebyshebv-Bienaymé inequality for any λ > 0 To exploit this inequality for the strong-drift limit, we can use a standard asymptotic estimate of the Kummer functions ( [103], §13.1.4) to infer that If we take in the inequality (A.9) the fixed value λ = (a 2 − x 2 )/2, then we infer that Because we therefore get for all x 0 ∈ (−a, a) that lim γ→∞,γ−t→∞ a −a dx P (γ) (x, t|x 0 , 0) = 1, (A.15) that is, the stochastic particle in this strong-drift limit γ → ∞ never escapes from the interval (−a, a) up to time t, even if t → ∞, as long as γ − t → ∞.
We can apply this result to the 1D model of SS in its dissipation-range description by taking a < 1, so that the interval (−a, a) lies in the regularization region |x| < 1. If we regularize the problem by taking so that the probability is positive to remain in a vanishingly small neighborhood of the origin as Re → ∞.
The velocity regularization required to stabilize the particle at x = 0 is admittedly "unnatural", corresponding to a potential well which is increasingly deep as Re → ∞. However, from the macroscopic inertial-range point of view the regularization region |x| ≤ Re −1/(1+h) disappears entirely from view as Re → ∞ and only the unstable singular drift field v * (x) is observed. In such "unnatural" models the exact solutionx ∞ (t) ≡ 0 with infinite waiting-time may occur with positive probability in addition to the two extremal solutionsx ± (t) with zero waiting-time. Note that these three solutions of the singular initial-value problem are exactly those that are invariant under the basic scaling symmetry (I.10).
We show here that the symmetric RG fixed point Q * in (III.35) attracts some non-symmetric initial probability distributions P 0 . For simplicity we assume that the velocity v(x) is symmetric, v(−x) = −v(x), although it should not be too difficult to remove that restriction.
(C.1) whereD = 1/Sc for (III.15) andD = 1/P e for (III.11). Here N i ∼ N (0, 1) is a pseudo-random sequence of standard normal random variables, which we obtained by generating uniform random variables with the Mersenne twister algorithm [104] and applying a 256-step Ziggurat method to obtain normal variables [105]. Statistics presented are empirical averages over N independent samples and empirical probability density functions P (x, t) of particle positions are further smoothed using Gaussian kernel-density estimators [106].
For each average calculated, weak convergence with respect to step size ∆t was verified by comparing the average for two step sizes that differ by a factor of 2. The result was deemed to be converged with respect to step size if the error was dominated by the number of samples used. For Figures 1, 3, 4, 5, 6, 8 we chose ∆t = 0.001, 0.0156, 0.005, 0.00333, 0.00333, 0.00333, respectively. The number of samples N was chosen to be large enough for the goal at hand. For Figures 3, 6, 8 of position PDF's N = 10 7 was found to be sufficient so that the error is smaller than the line width of the plot.
This number had to increased to N = 10 9 in order to evaluate the bifurcation diagram in Figure 4, since the PDF becomes nearly flat near the bifurcation and this greatly increases the error in estimation of the location of local extrema. See Figure 7. To identify those extrema, we split the x − t plane into parts such that each must have only one local extremum and then found the global maximum/minimum on each part. The error was determined by finding the intervals where the values of the PDF differ by less than 0.1% from the extremal value.
To arrive at the phase diagram in Figure 1 we integrated the Langevin model in the inertial-range form (III.11) up to time 1 for the set of Re and P e values shown and the variance for each was evaluated by Monte Carlo average with N = 2680 samples. Convergence was checked by increasing N and assessing by eye the difference in the phase diagram. Figure 5 corresponds conceptually to slices of the phase diagram in Figure 3 along lines of constant Sc, each with slope 1, and rescaled from inertial to dissipative units. However, the range of Re and P e was modified to better illustrate all of the scaling regimes. Furthermore, the number of samples was increased to N = 5 × 10 5 in order to ensure that the error in the estimated variance is smaller than the width of the plotted curves.
To generate Figures 3, 4, 6, 8, 7 we approximated the PDF's using kernel-density estimators [106]. In this scheme the PDF is estimated by the N -sample average: ) is a filter kernel with bandwidth h. We chose a Gaussian filter kernel and for each PDF the bandwidth was chosen as in [8] via the least sensitivity principle [107]. According to this principle one selects the bandwidth that has the least influence on the estimated PDF. For each PDF the bandwidth was selected individually by evaluating the L 1 norm of the increments of the PDF for varying h. Then the bandwidth that minimizes the increments was chosen.
Finally, the constant 0 in the theoretical curves in Figure 6 was calculated with the MATSLISE 2.0 software, which uses a high-order piecewise constant perturbation method to calculate both eigenvalues and eigenfunctions [83]. The value obtained was 0 . = 0.55447399023 with estimated error 1.3529×10 −11 . We also verified that the numerically obtained eigenfunction Ψ 0 (x) (not shown here) satisfies the asymptotics (III.76) for |x| 1. To confirm that the persistent stochasticity observed in Figure 6 is a consequence of the Hölder regularity exponent h < 1, we have carried out an identical analysis for the same model with only the exponent changed to h = 1 and with all other parameters the same. Note that the regularized velocity field (I.6) for h = 1 reduces to a linear drift v(x) = Ax, corresponding to diffusion in an inverted parabolic potential. In Figure 8 we plot for this h = 1 case the inertial-range PDF's of particle po-sitionsP (x,t) for the same increasing sequence of values Re = P e as in Figure 6. Unlike there, however, we see now that the PDF's converge to a single delta-function at x = 0, corresponding to a deterministic limit. This null result verifies that the stochasticity in the infinite-Re, P e limit is not a trivial consequence of the instability of the fixed point x = 0 but requires the quasi-singular velocity with exponent h < 1. The results obtained in our model with h = 1 are typical of smooth dynamical systems, which may exhibit deterministic chaos with exponential error growth but for which determinism is restored at any finite timet in the zero-noise limit.