text:
# Resonant nucleation of spatio-temporal order via parametric modal amplification

## I Introduction

The emergence of spatio-temporal ordered structures in nonlinear systems is an ideal laboratory for investigating the trend toward complexification observed in nature at the physical, chemical, and biological level walgraef . It is known that localized, or spatially-bound, order emerges in systems which are out of thermodynamic equilibrium, at the expense of growing overall disorder bar-yam . Examples can be found in hydrodynamics, in networks of chemical reactions cross , in phase transitions and critical phenomena phase ; magnets ; degennes , and in living organisms schnerb . For these localized ordered structures to survive, they must interact with an external environment, which maintains nonequilibrium conditions: in general, an isolated nonlinear system will eventually reach equilibrium, maximizing its entropy to the detriment of localized ordered structures.

In other words, if localized ordered structures emerge during the approach to equilibrium of a closed system, they will eventually disappear as energy becomes increasingly equipartitioned amongst the system’s many degrees of freedom. The emergence of localized ordered structures is thus a nonequilibrium phenomenon, usually corresponding to an unequal partitioning of energy: if they are long-lived, they will serve as bottlenecks to equipartition. \[An exception occurs for 1d field theories with SDW, where configurations such as kink-antikink pairs correspond to the equilibrium.\] Following the pioneering work of Fermi, Pasta, and Ulam in the early 1950s ford , several recent works have investigated the possibility that discrete kink-like structures may act as bottlenecks for equipartition in chains of nonlinear oscillators discrete ; bottlenecks .

In field theory and cosmology, most of the interest on ordered configurations has focused on topological (kinks, vortices, monopoles) rajaraman ; vilenkin  or nontopological (nontopological solitons, Q-balls) NTSs  static solutions of the equations of motion. Those configurations can be boosted towards each other with a certain velocity and made to interact, to model nonlinear composite particle interactions. A well-known example are 1-dimensional breathers, time-dependent bound states that have been shown to result from interactions of kink-antikink pairs campbell . However, this is not surprising: any two or more static solitonic configurations with short-range attractive interactions can be prepared to generate spatially-bound time-dependent structures digal . This is to be contrasted with the spatially-bound, time-dependent structures of interest here, that emerge spontaneously during the nonlinear evolution of a system. We have recently found such structures in 2+1-dimensional simulations of nonlinear scalar field models gleiser-howell . Such structures were identified with oscillons, long-lived time-dependent localized field configurations that have been previously found deterministically in a variety of physical systems, ranging from field theory models with amplitude-dependent nonlinearities gleiser ; bettinson , vertically vibrated grains umbanhowar , and stellar interiors umurhan .

In this paper, we will examine the emergence of spatially-bound, time-dependent field configurations during the approach to equilibrium of (2+1)-dimensional scalar field models with typical double-well interactions. These interactions are controlled by a tunable external parameter, which can be interpreted as an extra interaction that may be switched on or off at will. In the parlance of Ising magnetic systems, the external interaction may be a spatially-homogeneous magnetic field, $`(t)`$. We will investigate both the symmetric double-well (SDW) – extending previous results of ref. gleiser-howell  – as well as provide new results for the asymmetric case (ADW). The interest in examining both cases stems from the fact that the ensuing dynamics will be quite different: while for the SDW the degeneracy in free energy density dictates that no critical nucleation is possible, this will not be the case for the ADW. In fact, for the SDW we observe the synchronous emergence of long-lived localized configurations (oscillons), which eventually disappear as the system reaches equipartition gleiser-howell . We will argue that these configurations act as bottlenecks for equipartition, concentrating the field’s energy in a narrow band of long-wavelength modes. For the ADW, and with the system initially quenched to the metastable state, the synchronous emergence of oscillons may trigger the nucleation of a critical droplet of the lower free-energy phase, which grows to complete the phase transformation. We will show that, for a range of temperatures, the time-scales associated with this process are orders of magnitude faster than the typical statistical homogeneous nucleation, where a system adiabatically cooled into its metastable phase decays to its lowest free energy phase with typically exponentially-large time scales langer ; domb . In fact, the present treatment, with an “instantaneous” quench (faster than other typical time-scales in the system), and the arbitrarily slow cooling implicit of homogeneous nucleation calculations, represent the two extreme cases of a whole spectrum of possible cooling rates. Since we will argue that the emergence of coherent structures is due to the resonant amplification of parametric excitations of certain field modes, we will refer to this nucleation mechanism as resonant nucleation. Anticipating one of our results, for field models described by amplitude-dependent scalar nonlinearities, resonant nucleation of spatio-temporal order will be relevant whenever the cooling time scale $`\tau _{\mathrm{cool}}`$ is faster than the equilibration time-scale of the longest wavelength in the system, $`\tau _{\mathrm{eq}}^0`$. \[The superscript “0” is a reminder that often the longest wavelength is associated with the zero mode of the system, which has the slowest relaxation time-scale.\] We note that we will always work far from the critical point where, of course, $`\tau _{\mathrm{eq}}^0\mathrm{}`$. For our recent work on the approach to criticality please see ref. howell .

The paper is organized as follows. In the next section we discuss the general properties of the model (a general Ginzburg-Landau model), the initial conditions, and the numerical implementation of the dynamics. In section 3 we present a mean field (homogeneous Hartree) approximation adequate to describe the system in the presence of small fluctuations about equilibrium. The fact that the homogeneous Hartree approximation breaks down as nonperturbative excitations become progressively more important will provide us with useful information about the emergence of ordered structures. In section 4 we argue that the emergence of spatio-temporal order is a consequence of parametric resonance, induced by the nonlinear coupling between the field’s zero mode – the forcing agent – and smaller-$`k`$ modes. We also provide entropic arguments in support of the statement that the emergence of coherent spatio-temporal structures delays equipartition of the energy and, thus, the approach to equilibrium. In section 5, we extend our investigation to ADWs. For a range of temperatures \[or initial energies, see below\], we measure the time scales for the decay of the metastable state as a function of barrier height, contrasting them with those computed from homogeneous nucleation theory. We identify two possible mechanisms for the decay: for small asymmetries, a critical nucleus is produced by the coalescence of two or more oscillons, while for larger asymmetries a single oscillon evolves into a critical nucleus. We show that, for a range of temperatures associated with the presence of oscillons in the system, the time scale associated with the decay obeys a simple power law, $`\tau _{\mathrm{RN}}[E_b/k_BT]^B`$, where $`E_b/k_BT`$ is the free-energy barrier of a critical fluctuation, and the exponent $`B`$ satisfying $`2.44<B(T)<3.36`$. We conclude in section 6 with a summary of our results and outlook for future work.

## II The nonlinear model and its implementation

We start by introducing the model we will work with and its statistical and numerical implementation. We consider a (2+1)-dimensional real scalar field (or scalar order parameter) $`\varphi (𝐱,t)`$ evolving under the influence of an on-site potential $`V(\varphi )`$. If the system does not interact with an external bath, the continuum Hamiltonian is conserved and gives the total energy of the system, for a given field configuration,

$$H[\varphi ]=d^2x\left[\frac{1}{2}(_t\varphi )^2+\frac{1}{2}(\varphi )^2+V(\varphi )\right],$$
(1)
where

$$V(\varphi )=\frac{m^2}{2}\varphi ^2\frac{\alpha }{3}\varphi ^3+\frac{\lambda }{8}\varphi ^4$$
(2)
is the potential energy density. The parameters $`m`$,$`\alpha `$, and $`\lambda `$ are positive definite and temperature independent. Before we go any further, a note on our choice of potential. One could as easily have chosen a typical Ginzburg-Landau (GL) potential, $`V(\varphi )=\frac{a}{2}\varphi ^2+\frac{\lambda }{8}\varphi ^4(t)\varphi `$, with $`(t)`$ an external homogeneous magnetic field. The GL potential can be transformed into the potential of eq. 2 by a field shift, $`\varphi ^{}=\varphi +C(;a,\lambda )`$, where $`C`$ is the solution of the cubic equation $`C^32\frac{a}{\lambda }C+\frac{2H}{\lambda }=0`$, $`\alpha =\frac{3}{2}\lambda C`$, and $`m^2=a+\frac{3}{2}C^2\lambda `$. We chose the potential of eq. 2 for convenience. Our control parameter is $`\alpha `$. We note that when $`\alpha =0`$ the potential is an anharmonic single well. In this case, there are no spinodal instabilities. (By spinodal instabilities we mean $`V^{\prime \prime }(\varphi )<0`$ for certain values of $`\varphi `$, which may result in the unstable growth of modes.) When $`\alpha 3/2`$ (in scaled units, see below) the potential can assume a typical double-well shape \[SDW for equality and ADW for $`\alpha >3/2`$\], and spinodal instabilities are possible. We will thus refer to the interactions modeled by $`\alpha `$ as destabilizing interactions (DIs). Changing $`\alpha `$ from an initial value of zero to a finite positive value can also be thought of as an “instantaneous” quench: the system, initially prepared in a single well, is “tossed” into a double-well potential within a time-scale much faster than any equilibration time-scale for longer wavelengths in the system, those with $`k^22\pi V^{\prime \prime }(\varphi =0)`$, where $`k`$ is the mode’s wave number. We will say more about this later on.

It is helpful to introduce the dimensionless variables $`\varphi ^{}=\varphi \sqrt{\lambda }/m`$, $`x^{}=xm`$, $`t^{}=tm`$, and $`\alpha ^{}=\alpha /(m\sqrt{\lambda })`$ (henceforth dropping the primes). When $`\alpha \alpha _c=3/2`$ the two minima are degenerate (SDW). If $`\alpha 0`$, the $`𝒵_2`$ symmetry is explicitly broken. We will only be concerned here with the cases modeled by $`\alpha \alpha _c`$. For $`\alpha >\alpha _c`$ the minimum at $`\varphi =0`$ becomes the local (or metastable) minimum (ADW) and $`\varphi _+=\alpha +\sqrt{\alpha ^22}`$ is the global minimum.

### II.1 Statistical Implementation and Initial Conditions

A statistical description of the state of the system is given by a set of time-dependent observables $`\{O^{(1)}(t),O^{(2)}(t),\mathrm{}\}`$, each one calculated as an ensemble average over $`N`$ realizations, or microstates. For example,

$$O^{(1)}(t)=N^1\underset{i=1}{\overset{N}{}}O_i^{(1)}(t).$$
(3)
In most cases, an observable within a particular microstate is itself an average over the area of the system, such as the average value of the field $`\varphi _i(t)=A^1_Ad^2x\varphi _i(𝐱,t)`$ for the particular microstate $`i`$. Most of the observables considered will have this “double” averaging. Such is the case for $`\varphi (t)`$, $`\pi (t)`$ (the field’s conjugate momentum), and the corresponding hierarchy of correlation functions. It will be noted when this is not the case, particularly when quantities are no longer assumed to be translationally invariant, such as those describing localized configurations.

We prepare the field such that it is in thermal equilibrium at temperature $`T`$ in the single-well ($`\alpha =0`$) potential at some time prior to the quench at $`t=0`$, when $`\alpha `$ is set to a nonzero value. When one is faced with a nonlinear equation of motion such as eq. 4 below, the thermal probability distribution for $`\varphi (𝐱,0)`$ cannot be obtained exactly. To generate the initial conditions for each of the microstates in the ensemble, we instead couple the field in the single-well potential to a heat bath via a generalized Langevin equation with additive noise,

$$\ddot{\varphi }+\gamma \dot{\varphi }^2\varphi =\varphi \frac{1}{2}\varphi ^3+\xi ,$$
(4)
where the viscosity coefficient $`\gamma `$ is related to the stochastic force of zero mean $`\xi (𝐱,t)`$ by the fluctuation-dissipation relation,

$$\xi (𝐱,t)\xi (𝐱^{},t^{})=2\gamma T\delta ^2(𝐱𝐱^{})\delta (tt^{}).$$
(5)
Eq. 5 defines the two-point correlation function of the stochastic force, assumed here to be Markovian. Since this force will only be used during the preparation of the initial thermal state – after the quench the coupling with the bath is turned off ($`\gamma 0`$) – this assumption is sufficient. \[However, it would be quite interesting to examine the effects of multiplicative noise in the nonequilibrium dynamics of this system when the contact with the thermal bath is kept throughout its evolution. \[The role of multiplicative noise on the onset of synchronization in 1+1-dimensional fields has been recently investigated by Munõz and Pastor-Satorras munoz .\] We note that even though we are preparing each microstate in an initial thermal state, there is a one to one correspondence between the temperature and the average energy of the field. In fact, we can refer to the ensemble prepared as above as a microcanonical ensemble, since each microstate has energy $`E\pm \delta E`$, where $`\delta E/E<10^3`$. Such numerical accuracy in energy matches that obtained with other methods used preparing the initial state, including the popular Monte Carlo sampling. It simply expresses the fact that, for large enough lattices that approach the thermodynamic limit, there is an equivalence between canonical and microcanonical ensembles chandler .

### II.2 Numerical Implementation

The full evolution of the scalar field for each microstate is simulated using a staggered leapfrog integration of the equation of motion. It is implemented on a square lattice with periodic boundary conditions. Integrations are performed with $`\mathrm{\Delta }x=0.2`$ and $`\mathrm{\Delta }t=0.02`$ with 1024 lattice sites per side. The field is then defined over an area with sides of length $`L=204.8`$. The leapfrog scheme replaces second-order derivatives with finite-difference expressions that are accurate to $`𝒪(\mathrm{\Delta }x^2,\mathrm{\Delta }t^2)`$. The lattice spacing and system size are chosen so that computations can be carried out in a timely manner, while the numerical values of any observables and their corresponding continuum values are in good agreement. For example, after the quench, the energy in this closed system is conserved by more than 1 part in $`10^3`$. Results are independent of lattice size, provided large enough lattices are used. Finally, all of the observables are then averaged over an ensemble of $`N=100`$ microstates.

## III Dynamics within the perturbative realm: Homogeneous Hartree Approximation

### III.1 Setting up the initial configuration

We start our investigation using a familiar perturbative method, the leading order or Homogeneous Hartree Approximation (HHA). An excellent account can be found in bonini , where successive truncations within the hierarchy of correlation functions were used to simplify the equation of motion, all the while capturing the essential aspects of the dynamics. We will borrow many of the same methods and notation found in this reference. The HHA assumes that the statistical distribution of the fluctuations in the field and its momentum remain Gaussian throughout the evolution of the system. We will see that this is indeed an excellent approximation in the following cases: i) just prior to the quench (that is, the initial state of the field in a single well) for all temperatures considered; ii) at all times after the quench when the temperature is low; iii) at very early times after the quench when the temperature is high. We start by considering the system before the quench is implemented ($`\alpha =0`$), that is, the dynamics with the potential $`V(\varphi )=\frac{1}{2}\varphi ^2+\frac{1}{8}\varphi ^4`$ boyanovsky:03 .

A simple way to implement the HHA in this situation is to factor the cubic term in eq. 4 as $`\varphi ^3(𝐱,t)3\varphi ^2\varphi (𝐱,t)`$, where $`\varphi ^2`$ is the translationally invariant mean-square variance of the field. Since the field is in thermal equilibrium, $`\varphi ^2`$ is time-independent. In this approximation, the equation of motion becomes linear,

$$\ddot{\varphi }^2\varphi =\left(1+\frac{3}{2}\varphi ^2\right)\varphi .$$
(6)
The Fourier transform of the field is a complex quantity given by

$$\overline{\varphi }(𝐤,t)=\frac{1}{(2\pi )^2}d^2x\varphi (𝐱,t)e^{i𝐤𝐱}.$$
(7)
The corresponding equation of motion for $`\overline{\varphi }(𝐤,t)`$ in $`k`$-space is \[with $`k=(k_x^2+k_y^2)^{1/2}`$\]

$$\ddot{\overline{\varphi }}=\omega ^2(k)\overline{\varphi },$$
(8)
with the effective frequency squared

$$\omega ^2(k)=k^2+1+\frac{3}{2}\varphi ^2=k^2+m_H^2.$$
(9)
We have introduced the Hartree mass $`m_H^2=1+\frac{3}{2}\varphi ^2`$. Notice that when there are fluctuations in this field ($`\varphi ^2>0`$) the effective Hartree mass is greater than the “bare” mass ($`m^2=1`$) given by the potential in eq. 2. In this approximation, the modes evolve independently of one another, just as they would in a linear truncation of the equation of motion, only now each interacts with the homogeneous mean-field background. A higher-order approximation is needed, however, to implement direct scattering between different modes.

If the system is in thermal equilibrium, the momentum modes in $`k`$-space will satisfy

$$|\overline{\pi }(𝐤)|^2=T.$$
(10)
It then follows from eq. 8 that the corresponding power spectrum for the field is

$$|\overline{\varphi }(𝐤)|^2=\frac{T}{k^2+m_H^2}.$$
(11)
To obtain the equilibrium quantity $`\overline{\varphi }^2`$, one integrates eq. 11 over the entire $`k`$-space and divides by the area. (Equivalently, one can integrate the corresponding two-point correlation function over $`x`$-space to obtain $`\varphi ^2`$, since Parseval’s theorem guarantees that the two quantities are the same.) One must be careful here, since the lattice results depend logarithmically on the ultraviolet cutoff (lattice-spacing). We refer the reader to ref. borrill , where a method is provided to handle such dependence perturbatively.

In figures 1 and 2 we show the probability distributions for the initial state of the field and its momentum for various temperatures. These distributions are obtained numerically by collecting the values of $`\pi (𝐱,0)`$ and $`\varphi (𝐱,0)`$ of every microstate into bins of width $`\mathrm{\Delta }\pi =\mathrm{\Delta }\varphi =0.01`$, and then normalizing so that the sum over each distribution is unity. Not all data points are shown, making these figures more readable. Also shown are the Gaussian probability distributions (solid lines) that are predicted from the Boltzmann-Gibbs formalism,

$`P(\pi )=`$ $`{\displaystyle \frac{1}{\sqrt{2\pi T}}}e^{\frac{\pi ^2}{2T}}`$ (12)
$`P(\varphi )=`$ $`{\displaystyle \frac{1}{\sqrt{2\pi aT}}}e^{\frac{\varphi ^2}{2aT}}.`$ (13)
The fitting parameter $`a0.51`$ is obtained by comparing the area (and ensemble)-averaged two-point corrrelation function of $`\varphi `$. The variance in figure 2, $`\varphi ^2`$, can be fitted linearly, $`\varphi ^2=aT`$. The agreement between the numerical results and the theoretical predictions confirms the accuracy of the lattice dynamics and the applicability of the HHA for low enough temperatures. The probability distribution of $`\pi (𝐱,0)`$ is exact, while that of $`\varphi (𝐱,0)`$ is indeed close to Gaussian for the temperatures considered. We also show, in figure 3, the power spectrum of the field modes, $`|\overline{\varphi }(𝐤,0)|^2`$, and compare it with the prediction from eq. 11, valid for the HHA.

### III.2 Dynamical evolution and breakdown of Hartree approximation

In this subsection we examine the dynamics of the field after the quench. We will limit our investigation to the case where $`\alpha \alpha _c=3/2`$ at $`t=0`$, that is, the potential is switched from a single well to a SDW. The dynamics for ADW potentials will be examined in section 5.

At $`t=0`$ the coupling to the bath is set to zero, and the dynamics is thus conservative thereafter. The system is closed, and can be described by a microcanonical ensemble of energy $`E\pm \delta E`$, where $`\delta E/E1`$. As we remarked before, for large enough systems in equilibrium, as is the case here for $`t<0`$, a description with fixed $`E`$ is equivalent to one with fixed $`T`$. We will use this initial temperature to characterize the system. It is interesting to note that in the thermodynamic limit, an instantaneous quench to any positive value of $`\alpha `$ will not change the energy of the system. The field is initially symmetric in its distribution and the quench introduces an odd term into the Hamiltonian. Instead, the energy is redistributed between its kinetic, gradient, and potential parts, and as a result the final equilibrium temperature will necessarily differ from the initial temperature. In this investigation, however, this change is always less than $`0.15`$ percent. For open systems, where the system remains in contact with a heat bath, we refer the reader to ref. gleiser-howell .

Due to the initial thermal distribution, defined entirely by the temperature $`T`$, the sudden introduction of the cubic term into the potential actually has a minor effect on most of the field. In general, there is little net change in the forces present throughout the system, since the curvature of the single- and double-well potentials are the same at $`\varphi =0`$ and the majority of the field is indeed close to this value (see figure 2). That is, the linear equation of motion remains unchanged for the small-amplitude fluctuations that represent the peak of the field’s probability distribution. The few locations where the field is more responsive to the perturbation correspond to regions with relatively large absolute values of $`\varphi `$. These values define the tails of the probability distribution. The new nonlinear term present in the potential breaks the symmetry in the equation of motion, and this second-order effect causes larger-amplitude field fluctuations with negative (positive) values to experience a larger (smaller) restoring force than before. Overall, a net positive force is imposed on the symmetrically prepared field and, as a consequence, its average starts at $`\varphi =0`$ and then begins to move toward $`\varphi _{\mathrm{max}}`$, the maximum of the potential barrier between the two states.

We will see that $`\varphi (t)`$ is initially driven more out of equilibrium than all of the fluctuations about it. This is not surprising, as longer wavelength modes have longer relaxation time-scales langer . It is natural, then, to distinguish the dynamics of the two from each other. First, define $`\varphi _{\mathrm{ave}}(t)\varphi (t)`$ and let $`\delta \varphi (𝐱,t)`$ represent the fluctuations about this average. Then,

$$\varphi (𝐱,t)=\varphi _{\mathrm{ave}}(t)+\delta \varphi (𝐱,t),$$
(14)
where $`\delta \varphi (𝐱,t)=0`$ by definition. In the HHA, one assumes that these fluctuations and their momenta remain Gaussian in their distribution. One then obtains an equation of motion for $`\varphi _{\mathrm{ave}}(t)`$,

$$\ddot{\varphi }_{\mathrm{ave}}=\frac{V_{\mathrm{eff}}}{\varphi _{\mathrm{ave}}},$$
(15)
where we have introduced the effective potential for $`\varphi _{\mathrm{ave}}`$,

$$V_{\mathrm{eff}}(\varphi _{\mathrm{ave}},\delta \varphi ^2)=\left[1m_\mathrm{H}^2(t)\right]\varphi _{\mathrm{ave}}+\frac{1}{2}m_\mathrm{H}^2(t)\varphi _{\mathrm{ave}}^2\frac{\alpha _\mathrm{c}}{3}\varphi _{\mathrm{ave}}^3+\frac{1}{8}\varphi _{\mathrm{ave}}^4,$$
(16)
that depends explicitly on the mean-square variance of the fluctuations through the time-dependent Hartree mass $`m_\mathrm{H}^2(t)=1+\frac{3}{2}\delta \varphi ^2(t)`$.

$`V_{\mathrm{eff}}`$ at $`t=0`$ is shown in figure 4 for various temperatures. For sufficiently low initial temperatures (how low will become clear soon), $`\varphi _{\mathrm{ave}}(t)`$ will oscillate about the LHS minimum of the effective potential, $`\varphi _{}=1\sqrt{13\delta \varphi ^2}`$. It will then begin to relinquish its energy to higher modes through nonlinear scattering, and its amplitude will decay. It is possible, then, to approximately bound this region of oscillation from the set of turning points for $`\varphi _{\mathrm{ave}}(t)`$ \[with $`\delta \varphi ^2`$ retaining its $`t=0`$ value\],

$$\varphi _{\mathrm{ave}}^{\mathrm{tp}}=\{0,\mathrm{\hspace{0.17em}1}\pm \sqrt{16\delta \varphi ^2},\mathrm{\hspace{0.17em}2}\}.$$
(17)
These values give the location of $`\varphi _{\mathrm{ave}}(t)`$ where $`V_{\mathrm{eff}}=0`$, and $`\pi _{\mathrm{ave}}(t)=_t\varphi _{\mathrm{ave}}(t)`$ vanishes. When $`\delta \varphi ^2<1/6`$ there are four such turning points because of the symmetry about $`\varphi _{\mathrm{max}}`$. The primary region of interest is when $`\varphi _{\mathrm{ave}}^{\mathrm{tp}}\varphi _{\mathrm{inf}()}`$, when fluctuations begin to probe the spinodal region of $`V_{\mathrm{eff}}`$. By using eq. 17, we can predict that $`\varphi _{\mathrm{ave}}(t)`$ will probe this region when $`\delta \varphi ^2=1/9`$, or $`T0.22`$. Clearly, this is also where the leading-order Hartree approximation fails.

It is interesting to consider also the case when $`\delta \varphi ^2=1/3`$, or $`T0.65`$. There are only two turning points for $`\varphi _{\mathrm{ave}}`$, since $`V_{\mathrm{eff}}`$ is now a single-well potential symmetric about $`\varphi =\varphi _{\mathrm{max}}`$. The fluctuations are large enough that the local maximum in the bare potential has little influence in the dynamics, and $`\varphi _{\mathrm{ave}}=\varphi _{\mathrm{max}}`$. The fluctuations actually probe beyond both minima of the potential. This case is similar to the phenomena of symmetry restoration in continuous phase transitions, and so in this context $`T0.65`$ would be considered the critical temperature, $`T_c`$. For all higher temperatures, $`\varphi _{\mathrm{ave}}=\varphi _{\mathrm{max}}`$ phase ; gleiser-mix . In this work, we are only concerned with temperatures well below $`T_c`$, in fact below $`T=0.32`$, where $`\varphi _{\mathrm{ave}}^{\mathrm{tp}}=\varphi _{\mathrm{max}}`$ (see figure 4).

In order to further quantify the dynamics in the HHA, we derive the set of differential equations for the time evolution of the coupled two-point correlation functions. We define the translationally invariant quantities

$`G_{\varphi \varphi }(𝐱𝐲,t)`$ $`=`$ $`\delta \varphi (𝐱,t)\delta \varphi (𝐲,t)`$
$`G_{\pi \pi }(𝐱𝐲,t)`$ $`=`$ $`\delta \pi (𝐱,t)\delta \pi (𝐲,t)`$ (18)
$`G_{\pi \varphi }(𝐱𝐲,t)`$ $`=`$ $`{\displaystyle \frac{1}{2}}\left(\delta \pi (𝐱,t)\delta \varphi (𝐲,t)+\delta \varphi (𝐱,t)\delta \pi (𝐲,t)\right).`$
Upon taking the time-derivative of each expression, the equations of motion can then equivalently be written in $`k`$-space. They are

$`\dot{G}_{\varphi \varphi }(k,t)`$ $`=`$ $`2G_{\pi \varphi }(k,t)`$
$`\dot{G}_{\pi \pi }(k,t)`$ $`=`$ $`2\overline{\omega }^2G_{\pi \varphi }(k,t)`$ (19)
$`\dot{G}_{\pi \varphi }(k,t)`$ $`=`$ $`G_{\pi \pi }(k,t)\overline{\omega }^2G_{\varphi \varphi }(k,t),`$
where $`\overline{\omega }^2(k)=k^2+V_{\mathrm{eff}}^{\prime \prime }[\varphi _{\mathrm{ave}}(t),\delta \overline{\varphi }^2(t)]`$.

These equations, coupled with

$`\dot{\pi }_{\mathrm{ave}}`$ $`=`$ $`{\displaystyle \frac{V_{\mathrm{eff}}}{\varphi _{\mathrm{ave}}}}`$ (20)
$`\dot{\varphi }_{\mathrm{ave}}`$ $`=`$ $`\pi _{\mathrm{ave}},`$
give the full dynamics of $`\pi _{\mathrm{ave}}(t)`$ and $`\varphi _{\mathrm{ave}}(t)`$ and their corresponding two-point correlation functions $`\delta \pi ^2(𝐤,t)`$ and $`\delta \varphi ^2(𝐤,t)`$ in the HHA. These two groups of expressions are likewise coupled through the quantities $`\overline{\delta \varphi }^2(t)=\delta \varphi ^2(t)`$ and $`\varphi _{\mathrm{ave}}(t)`$. The initial conditions are given by eqs. 10, 11, and also $`\pi _{\mathrm{ave}}(0)=\varphi _{\mathrm{ave}}(0)=0`$.

These five first-order equations were solved numerically using a fourth-order Runge-Kutta algorithm. The same lattice spacing from the leap-frog algorithm was used. The calculation of $`\overline{\delta \varphi }^2(t)`$ is obtained by numerically averaging the radial function $`G_{\varphi \varphi }(k,t)`$ over the disk in $`k`$-space with radius $`k^2=k_x^2+k_y^2`$, bounded by $`k_{\mathrm{max}}=\pi /\mathrm{\Delta }x`$. There is a slight discrepancy between the area of rectangular phase space used in the two-dimensional numerical integrations and the circular phase space used in this method. In light of the previous discussion regarding the dependence of the correlation functions on the lattice spacing, this difference could possibly introduce error in the calculations of $`\overline{\delta \varphi }^2(t)`$. To compensate for the smaller area used in the HHA, the initial conditions for $`\overline{\delta \varphi }^2(k,0)`$, given by eq. 11, are multiplied by a normalization factor of 1.041 so that the radial integration gives the expected average $`\overline{\delta \varphi }^2=aT`$ at $`t=0`$. The results will confirm the validity of this normalization method.

Shown in figure 5 is the numerical time evolution of $`\varphi _{\mathrm{ave}}(t)`$ compared to that obtained through the HHA, for various initial temperatures. The agreement is clearly excellent at lower temperatures, where $`\varphi _{\mathrm{ave}}(t)`$ undergoes almost-linear oscillations about the minimum of $`V_{\mathrm{eff}}`$: there is very little scattering between modes. At $`T0.18`$ scattering begins to take place, allowing for a slow transfer of energy from low-$`k`$ modes to higher-$`k`$ modes not accounted for in the HHA. We will soon show that this is also the temperature where the synchronous emergence of spatio-temporal order begins to take place. Notice also the dependence of the oscillation frequencies on the initial thermalization temperature.

Figures 6 and 7 compare the numerical time-evolution of $`\delta \pi ^2(t)`$ and $`\delta \varphi ^2(t)`$ with the corresponding HHA predictions. The quantities are normalized to unity at $`t=0`$ through $`\pi ^2=T`$ and $`\varphi ^2=aT`$, respectively. The agreement is again excellent at lower temperatures. At $`T=0.18`$, however, it is evident that energy is being fed into the fluctuations of the field. The full numerical evolution depicts the amplitude of oscillation of $`\delta \pi ^2(t)`$ decreasing at a rate larger than that predicted by the HHA. Similarly, the amplitude of oscillation of $`\delta \varphi ^2(t)`$ increases at a larger rate. Again, this discrepancy is due primarily from the approximation’s inability to account for direct scattering between modes.

In figure 8 we show the transfer of energy between different $`k`$-modes. Displayed is the time evolution of $`|\overline{\delta \varphi }(k,t)|^2`$ for indicative values of $`k`$. The initial thermalization temperature is $`T=0.22`$, corresponding to the value in which the majority of the field probes the spinodal region (cf. section IV). This is well above the region of temperatures in which the HHA is accurate. From the figure it is clear that scattering between low $`k`$ modes modes is taking place. The approximation is only accurate for shorter wavelength modes, due to the fact that they remain in equilibrium even after the quench. The breakdown of the HHA at longer wavelengths implies that important nonlinear effects are strongly influencing the behavior of the field. In the next section we investigate the emerging dynamics of these modes.

## IV Emergence of spatio-temporal order in a SDW model

### IV.1 Beyond the Homogeneous Hartree Approximation

The results of the last section should have made it clear that at high temperatures nonlinear scattering leads to an efficient energy transfer between modes. This energy transfer will eventually lead to the thermalization of the field in the SDW potential. As predicted by the HHA, at $`T0.18`$ there is consistent energy transfer from the zero-mode to these higher-$`k`$ modes. At these initial temperatures, $`\varphi _{\mathrm{ave}}(t)`$ undergoes damped oscillations about its effective minimum. For $`T=0.22`$, the amplitudes of these oscillations are large enough that $`\varphi _{\mathrm{ave}}(t)`$ just crosses the inflection point of the bare potential: roughly half of the field probes the spinodal region.

Figure 9 illustrates the transfer of energy from $`\varphi _{\mathrm{ave}}(t)`$ to the fluctuations $`\delta \varphi ^2(t)`$ for $`T=0.22`$. The quantity

$$C(t)\frac{\delta \varphi ^2(t)aT}{aT}$$
(21)
is introduced so that a visual comparison can be made between these two components of the field. Notice the correlation between the decaying amplitude of oscillation for $`\varphi _{\mathrm{ave}}(t)`$ and the increasing amplitude of the fluctuations. At $`t60`$ most of this energy transfer is complete, after which time the two quantities display smaller oscillations that are almost exactly out of phase with each other. For earlier times, $`0t60`$, a highly nonlinear scattering process occurs, as the zero-mode of the field pumps its energy into its neighboring small-$`k`$ modes.

Figure 10 shows the numerical time-evolution of $`|\overline{\delta \varphi }(k,t)|^2`$ for $`0k4`$ and $`0t60`$, at $`T=0.22`$. First, notice that growth occurs mainly for modes with $`0k0.5`$. Their maximum amplitude reaches a value that is up to 15 times greater than at $`t=0`$. The oscillations then begin to decay and, for $`t>60`$ (discussed later), the amplitudes of the modes take on a constant value, albeit still many times larger than their initial value. It is also apparent that for the time interval shown modes with $`k>1`$ are nearly time-independent, maintaining their initial thermal spectrum given by eq. 11. This is true of all large-$`k`$ modes present in the numerical simulation, up to $`|k|_{\mathrm{max}}15.7`$, which also indicates that the lattice spacing used has no effect on the emerging dynamics of the system.

We would now like to search for possible spatio-temporal ordering during the field’s approach to equilibrium. However, due to the noisy initial conditions, it is quite difficult to extract information about any emergent pattern, be it spatial and/or temporal. An added advantage to having an analytical expression for the initial thermal spectrum for all temperatures of interest (cf. eq. 11) should be clear: a Wiener filter is extremely successful at removing undesirable noise whenever its spectrum is known, all the while maximally preserving the nonperturbative fluctuations that emerge numericalrecipes .

### IV.2 Synchronous emergence of spatio-temporal order

With the filtered field, one can catalog and track the location of all local extrema in the system at each instant in time and for all initial temperatures. To distinguish the large-amplitude fluctuations from the long-wavelength thermal noise that is also passed by the filter, only the values of these extrema that are greater than the characteristic amplitude of the noise, $`\delta \varphi ^2(𝐱,t)>aT`$, are considered. With ample sorting, a library is compiled containing the description of every large-amplitude fluctuation present during the entire evolution of the system (the maximum time for the simulations is $`t_{\mathrm{max}}=2000`$). For each fluctuation the following attributes are recorded: its core location $`𝐱_\mathrm{f}`$, nucleation time $`t_{\mathrm{nuc}}`$, core amplitude $`\varphi _\mathrm{a}`$ (measured at $`𝐱_\mathrm{f}`$), radius $`R`$ (measured at half maximum), its total number of oscillations $`N_{\mathrm{osc}}`$, the period of each oscillation $`𝒯_N`$, and its lifetime $`\tau `$.

Each of these attributes is time-dependent, with the exception of $`N_{\mathrm{osc}}`$ and $`\tau `$. Almost every large-amplitude fluctuation oscillates anharmonically about the background value $`\varphi _{\mathrm{ave}}(t)`$, which is recognized by the fact that most of them reappear at approximately the same location with regular frequency. Therefore, in general it is possible to establish the relationship $`N_{\mathrm{osc}}=\tau /𝒯`$, where $`𝒯𝒯_N`$ is the average period of oscillation. The set of observables is further reduced by averaging $`𝐱_f`$, $`\varphi _\mathrm{a}`$, and $`R`$ over their corresponding values at the maximum of each oscillation for the lifetime of the fluctuation. Other than the nucleation time $`t_{\mathrm{nuc}}`$, each attribute is now time-independent. Hereafter, we will omit the $`\mathrm{}`$ notation and assume that these values are indeed averages.

With this information, a probability distribution function for the large-amplitude fluctuations is obtained for each initial temperature. This function, denoted as $`F_T(𝐱_f,t_{\mathrm{nuc}},\varphi _\mathrm{a},R,𝒯,\tau )`$, can be used to calculate averages over any of its variables. Similarly, one can integrate this function over a subset of variables to obtain a reduced distribution. In all cases, $`F_T`$ is normalized so that a numerical integration over all variables evaluates to unity. The bin widths of each variable (as any one variable of $`F_T`$ will be visualized as a histogram) will be specified below.

For the moment, we focus our attention on the characteristics of fluctuations early in the evolution of the field, $`0t100`$, and at an initial temperature $`T=0.22`$. The quantity

$$f_T(\varphi _\mathrm{a},\tau )=_Ad^2x_\mathrm{f}_0^{100}𝑑t_{\mathrm{nuc}}_0^{L/2}𝑑R_0^{t_{\mathrm{max}}}𝑑𝒯F_T$$
(22)
gives the reduced distribution function for the core amplitude and lifetime, and is shown in figure 11. The bin widths are $`\mathrm{\Delta }\varphi _\mathrm{a}=0.05`$ and $`\mathrm{\Delta }\tau =1`$. The darkest regions in the figure correspond to the most probable fluctuations. Particularly noticeable is that fluctuations with core amplitudes $`0.55\varphi _\mathrm{a}1.0`$ and lifetimes $`35\tau 65`$ dominate the distribution. One can also see that fluctuations with small amplitudes $`0.2\varphi _\mathrm{a}0.5`$ have mostly short lifetimes. There is a group of fluctuations, however, that have lifetimes $`80\tau 100`$ and thus remain present in the system for substantially long periods of time. In addition, their core amplitudes are confined to a narrower region $`0.7\varphi _\mathrm{a}0.9`$. Also marked on the horizontal axes is the location $`\varphi _\mathrm{a}=0.81`$, which corresponds to the difference between the local maximum and the left minimum of the effective potential for $`T=0.22`$, given by eq. 16. It follows that these long-lived oscillatory fluctuations have large enough amplitudes that their cores probe a narrow region centered about $`\varphi _{\mathrm{max}}`$. In fact, the large majority of long-lived fluctuations probe beyond the inflection point of $`V_{\mathrm{eff}}`$ at $`\varphi _{\mathrm{inf}()}0.4`$. They are clearly nonperturbative in nature and would not be seen within the HHA.

In figure 12 we show a sequence of snapshots of the filtered field within the time interval $`21.0t29.4`$. There are several features to be noted. \[Clearly, they are more striking with a full animation. We invite the reader to view simulations in ref. website .\] First, large-amplitude fluctuations seem to be of similar size (spatially-bound configurations) and spread throughout the lattice. Second, they emerge in near synchrony. As we will see below, this global ordering persists for the first few oscillation cycles. Third, as can be seen from figure 11, these configurations can be quite long-lived.

We identify these spatially-bound, long-lived configurations with oscillons, which have been originally found as deterministic solutions of models with amplitude-dependent nonlinearities gleiser ; bettinson . In these studies, the field was prepared with a bubble-like profile (typically a Gaussian or a tanh), and for large enough initial radius and amplitude, observed to fall into an oscillon configuration, characterized by a rapid oscillation of its core, as shown in figure 12. Using harmonic perturbation analysis about an oscillon solution, it has been shown that such configurations require a minimum radius (in 2 spatial dimensions for the SDW potential used here) sornborger ,

$$R_{\mathrm{crit}}=\frac{\sqrt{6}}{m}.$$
(23)
We will now extend this result to the general ADW potential. First, we parameterize the radially-symmetric oscillon configuration by a Gaussian gleiser ; sornborger

$$\varphi _{\mathrm{osc}}(r,t)=\varphi _a(t)\mathrm{exp}[r^2/R^2],$$
(24)
where $`\varphi _a(t)`$ gives the time behavior at the core of the configuration and $`R`$ is its radial dispertion, assumed here to be constant. (This assumption is equivalent to reducing the dynamics to one degree of freedom; it has been shown to provide results in excellent agreement with numerical simulations gleiser ; sornborger .) With this ansatz, a Lagrangian can be obtained and, from it, the equation of motion for $`\varphi _a(t)`$. Linear perturbations about the solution, $`\delta \varphi =\varphi _a(t)\varphi _0(t)`$, oscillate with frequency

$$\omega ^2(\varphi _0,\alpha ,R)=\frac{2}{R^2}+1\frac{4}{3}\alpha \varphi _0+\frac{3}{4}\varphi _0^2.$$
(25)
Note that $`\omega ^2`$ has a minimum at $`\varphi _0^c=\frac{8}{9}\alpha `$. There is an unstable bifurcation at

$$R_{\mathrm{crit}}^2=\frac{2}{\frac{16}{27}\alpha ^21}.$$
(26)
For the SDW, $`\alpha _c=3/2`$ and the instability is given by eq. 23. Oscillons are only possible for values of $`R`$ larger than $`R_{\mathrm{crit}}`$. Larger asymmetries imply smaller values of $`R_{\mathrm{crit}}`$.

In order to strengthen our claim that oscillons are nucleated during the nonequilibrium evolution of $`\varphi (𝐱,t)`$, the inset in figure 14 shows the distribution functions $`f_T(R)`$ and $`f_T(𝒯)`$ for the radii and period of oscillation, respectively. To make a comparison with deterministic oscillons, only fluctuations with $`\varphi _\mathrm{a}0.81`$ (i.e. those that probe beyond the local maximum of the potential) and $`\tau 80`$ are considered. That is, these are the lower bounds placed on the corresponding integrations over $`F_T`$. The bin widths are $`\mathrm{\Delta }R=\mathrm{\Delta }𝒯=0.1`$. The fitted curves are Gaussian functions with centers $`R_{\mathrm{ave}}=4.3`$ and $`𝒯_{\mathrm{ave}}=8.5`$ and widths $`\sigma _R=0.4`$ and $`\sigma _𝒯=1.8`$, respectively. To complete the comparison, consider a deterministic oscillon evolving in a time-independent effective potential that is representative of this case. From figure 9 a good approximation would be to take $`C(t)0.4`$, so that $`\delta \varphi 0.15`$. The effective mass can thus be obtained from the curvature of the left-hand well, given by eq. 16, and evaluates to $`m_{\mathrm{eff}}=0.73`$. Finally, from eq. 23, we find that the effective minimum radius for which the deterministic oscillon remains stable is $`R_{\mathrm{crit}}^{\mathrm{eff}}=\sqrt{6}/m_{\mathrm{eff}}3.37`$. Comparing this result with figure 14 we see that, indeed, the emerging oscillons satisfy this requirement remarkably well. \[Note that a proper comparison between oscillons that evolve in different potentials (i.e. for fields at different initial temperatures) always requires a suitable scaling of $`\varphi `$, $`𝐱`$, and $`t`$.\]

Also shown in figure 14 is the distribution of nucleation times for the emerging oscillons. This quantity gives the number of oscillons nucleated between $`t`$ and $`t+\mathrm{\Delta }t`$, with $`\mathrm{\Delta }t=1`$. It is obtained by calculating the reduced distribution function $`f_T(t_{\mathrm{nuc}})`$ integrated over all the remaining variables. The same lower limits are used for the integrals over $`\varphi _\mathrm{a}`$ and $`\tau `$, so that only long-lived large-amplitude fluctuations are considered. If there are $`N_{\mathrm{tot}}`$ such fluctuations, then the final expression for the number of oscillons nucleated as a function of time is

$$N(t)=N_{\mathrm{tot}}f_T(t_{\mathrm{nuc}}).$$
(27)
The sharp peaks in $`N`$ at early times, $`t<60`$, correspond to the synchronous emergence of oscillons. For $`t>60`$, this global ordering gives way to only local ordering, in which oscillons emerge at arbitrary times with similar probability. Even this local emergence disappears after approximately $`t500`$ (more on this later). Notice the correlation between the synchronous nucleation of oscillons and the energy loss from the zero-mode (c.f. figure 9) at early times, $`t<60`$. This phenomenon was observed for temperatures within $`0.16T0.26`$ (see below). Below this temperature range the HHA is valid and no coherent structure emerges. Above this temperature range the field separates into large, slowly evolving thin-walled domains, signalling the approach to criticality.

### IV.3 Synchronous emergence of order is due to parametric resonance

In order to understand the mechanism behind the observed synchronous emergence, we decompose the field as in eq. 14. The linearized equation satisfied by the Fourier modes of the fluctuations is

$$\ddot{\overline{\delta \varphi }}=\overline{\omega }^2(k)\overline{\delta \varphi },$$
(28)
where

$$\overline{\omega }(k)^2=k^2+V_{\mathrm{eff}}^{\prime \prime }\left[\varphi _{\mathrm{ave}}(t)\right].$$
(29)
In this approximation, the fluctuations obey a linear equation of motion with a time-dependent frequency determined by the value of $`\varphi _{\mathrm{ave}}(t)`$. Equations of this type, generalized Mathieu equations, are known to exhibit parametric resonance, which can lead to exponential amplification ($`\mathrm{exp}\eta t`$) in the oscillations of $`\delta \varphi (k,t)`$ at certain wavelengths. They have been the focus of much recent attention in studies of reheating of the universe after an inflationary expansion phase kofman . To verify that this is the mechanism behind the synchronous amplification of oscillon modes early in the field’s evolution, we use in eq. 29 the anharmonic solution for $`\varphi _{\mathrm{ave}}(t)`$ that would be obtained if the effective potential remained time-independent. That is, we require that $`\varphi _{\mathrm{ave}}(t)`$ oscillates in the effective potential as determined immediately after the quench (i.e., without energy loss), just as the parametric driving force in generalized Mathieu equations oscillates with constant amplitude. In this sense, we are considering the initial response of the fluctuations to the oscillations of $`\varphi _{\mathrm{ave}}(t)`$, without accounting for the subsequent transfer of energy between the two. Only with this approach can one achieve a constant exponential amplification of $`\overline{\delta \varphi }`$ and thus a good measure of resonance. Even though this approximation only strictly holds for early times, the results obtained will confirm its validity in describing the essential physics behind the observed synchronous emergence.

Figure 15 shows lines of constant exponential amplification rate $`\eta `$ of the fluctuations $`\delta \varphi (k,t)`$ for various $`k`$ and $`T`$. At low temperatures, $`T<0.14`$, no modes are ever amplified. As the temperature is increased, so is the amplitude of oscillation in $`\varphi _{\mathrm{ave}}(t)`$, gradually causing the band $`0<k0.48`$ to resonate. It is helpful at this point to consider what modes comprise a typical oscillon. To do so, we use a Gaussian ansatz for the oscillon with a fixed radius $`R`$, $`\varphi _{\mathrm{osc}}(r,t)\mathrm{exp}(r^2/R^2)`$ gleiser . The Fourier transform is itself a Gaussian with a characteristic radius $`K=2/R`$. Thus, a Gaussian in real space is generally comprised of the band $`0k2/R`$ in $`k`$-space. Earlier \[cf. figure 14\] it was found that the average radius of the emergent oscillons is $`R_{\mathrm{ave}}=4.3`$ at an initial temperature $`T=0.22`$, which in this approximation corresponds to the band $`0k0.47`$. This agrees very well with what is seen in figure 15. Although in this linear approximation the entire band does not resonate until higher temperatures, it is expected that in the full equation of motion the nonlinear coupling between modes will broaden the resonance band at lower temperatures.

We have observed that the synchronous emergence of oscillons occurs for initial temperatures $`0.16T0.26`$. The lower bound is directly related to the nonlinear scattering between modes, which is absent at low temperatures. For temperatures above $`T0.26`$, the appearance of large domains precludes the identification of oscillons. This is typical of the approach to criticality, where large spatial correlations in the field become more probable, destroying any local ordering. Thus, the observed quench-induced synchronous emergence of order is limited to temperatures sufficiently below $`T_c`$.

### IV.4 Oscillons as bottlenecks for equipartition

In the introduction, we mentioned that long-lived locally-ordered structures may serve as bottlenecks to equipartition, as they slow down the distribution of energy between the system’s various degrees of freedom. Below, we present evidence that oscillons do precisely that. A more formal proof is left for future work.

A good observable to follow is the time evolution of $`\overline{\pi }^2(k,t)`$ for different modes, which should become flat when thermal equilibrium is achieved and the equipartition theorem is satisfied. In highly dissipative systems (such as those governed by a Langevin equation of motion with a large viscosity coefficient), energy exchange between each mode and the thermal bath tends to dominate the dynamics. In the nonequilibrium closed system considered here, however, energy is only exchanged between the modes themselves, and the thermalization process turns out to be quite different.

In figure 16 we show the time evolution of $`\overline{\pi }^2(k,t)`$ for the values $`k`$=0.25, 0.5, 0.75, 1.0, 1.25, and 1.5 (top to bottom lines, respectively). The curves are time-averaged over a window $`\mathrm{\Delta }t=10`$ so that the large oscillatory behavior at early times ($`t<60`$) is suppressed. Also shown in the inset is the time evolution of $`\pi _{\mathrm{ave}}(t)`$, which exhibits the thermalization trend of the zero mode. $`\pi _{\mathrm{ave}}(t)`$ displays damped oscillatory motion before reaching its equilibrium value at $`t160`$. This is roughly the time in which the mode $`k=0.25`$ reaches its maximum average kinetic energy. The mode $`k=0.5`$ also increases from $`t=0`$ to a maximum value at about the same time before decaying. For modes $`k>0.5`$, however, the general trend is to slowly absorb energy as they asymptotically approach their equilibrium values from below on time-scales that can be much longer than the lower-$`k`$ modes. There are two clear trends: small modes (with $`k0.5`$) grow very quickly and then relax to their final equilibrium values, while larger modes (with $`k>0.5`$) slowly approach theirs. Thus, it is reasonable to conclude that most of the zero-mode’s initial energy goes into its neighboring modes, those responsible for oscillon-like configurations \[cf. figure 15\], and only later on it is shared amongst larger $`k`$ modes. This seems to suggest that there is a two-step energy cascade, a fast one at early times to small-$`k`$ modes, followed by a slower one to larger $`k`$ modes. The early emergence of oscillons is consistent with this scenario.

We can further quantify this argument by introducing a measure of the entropy of the system based on the partitioning of the kinetic energy, $`K(𝐤,t)=1/2|\overline{\pi }(𝐤,t)|^2`$:

$$S(t)=_Ad^2kp(𝐤,t)\mathrm{ln}p(𝐤,t),$$
(30)
where

$$p(𝐤,t)=\frac{K(𝐤,t)}{d^2kK(𝐤,t)}$$
(31)
is the probability distribution of the kinetic energy in $`k`$-space. $`S(t)`$ attains its maximum value when equipartition is satisfied, and $`p(𝐤,t)=A^1`$ where $`A`$ is the area of this space. This occurs both at the initial thermalization $`(t=0)`$ and final equilibrium states, since in these cases all modes carry the same fractional kinetic energy. At these times we recover the thermodynamic entropy of the system, $`S_{\mathrm{eq}}=\mathrm{ln}A`$, which is a measure of the number of accessible states in the system. \[The numerical implementation gives $`S_{\mathrm{eq}}=\mathrm{ln}N`$, where $`N`$ is the number of lattice points, or degrees of freedom.\] In figure 17 we show the change of $`S(t)`$ from the initial state, $`S(t=0)S(t)`$, for the system at $`T=0.22`$. At late times ($`t>100`$), we have found that the system equilibrates exponentially in a time-scale $`\tau 500`$. At early times, the localization of energy at lower $`k`$ modes, corresponding to the synchronous emergence of oscillons, prolongs this approach to equipartition. The inset of figure 17 shows the large oscillations in the early evolution of $`S(t)`$ (dotted line). Also shown (solid line) is the average between successive peaks of $`S(t)`$, with a plateau at $`30t100`$ that coincides with the maximum oscillon presence in the system. This is to be contrasted with an exponential approach to equilibrium. Thus, we see evidence that oscillon configurations indeed serve as early bottlenecks to equipartition, temporarily suppressing the cascading of energy from low-$`k`$ modes to higher modes.

## V Resonant Nucleation

We will now investigate the dynamics with the ADW potential, that is, for values of $`\alpha >3/2`$. Homogeneous nucleation predicts that, had the system been prepared in the metastable state at $`\varphi =0`$ at some temperature $`T`$, its decay rate per unit area would be controlled by the Arrhenius factor, $`\mathrm{\Gamma }_{\mathrm{HN}}\mathrm{exp}[E_b/T]`$, where $`E_b`$ is the free-energy barrier for the nucleation of a critical nucleus or bounce langer ; domb . (We set $`k_B=1`$.) It is often overlooked that this result is quite sensitive to how the initial state is prepared. It assumes that the bounce is nucleated within a homogeneous background, by adopting a perturbative Gaussian approximation for the evaluation of the partition function. In other words, it assumes that no large-amplitude perturbations are present in the system gleiser-heckler ; tetradis . It is also well-known that there are discrepancies between the theoretical prediction from homogeneous nucleation and a large array of experimental results concerning nucleation rates, from He<sup>3</sup> leggett  and nematic liquid crystals degennes  to martensitic materials and polycrystals offerman  to fluids bartell  and numerical simulations gleiser\_nuc . More often than not, theoretical estimates greatly underestimate nucleation rates. Surely, many of these discrepancies are related to limited knowledge of the nonequilibrium processes leading to nucleation at the atomic or mesoscopic scale for various materials. However, we propose here that in certain cases the problem is rooted on the detailed set up of the initial metastable state. In particular, the quench we propose in this work, which may be mimicked in experimental situations, generates an instability capable of greatly accelerating the decay of the metastable state. In fact, we have observe that, for a range of temperatures, the nucleation rate per unit area is not that of homogeneous nucleation, $`\mathrm{\Gamma }_{\mathrm{HN}}=C\mathrm{exp}[E_b/T]`$, but a power law, $`\mathrm{\Gamma }_{\mathrm{RN}}=C(E_b/T)^B`$, where $`E_b(\alpha )`$ is the critical nucleus free energy and the exponent $`B`$ is weakly temperature-dependent. RN stands for Resonant Nucleation. The power law applies precisely within the temperature range where we observe the synchronous emergence of oscillons at early times, $`0.18T0.24`$, as we show next.

The change in $`\alpha `$ will not modify the basic parametric amplification mechanism discussed above in the context of the SDW: for a range of temperatures, oscillons are still nucleated in synchrony at early times. The difference is that, for $`\alpha >\alpha _c`$, the system must reach the global free-energy minimum at $`\varphi _+(\alpha )`$. \[Of course, for large enough temperatures, $`T>0.26`$, and/or small energy barriers the system will go straight into the global minimum, a cross-over decay. We are not interested in these cases.\]

In the usual treatment, the critical nucleation is described as the saddle-point in field-configuration space, the unstable direction dictated by the negative eigenvalue associated with the growth of the critical nucleus. This is given by the solution of the equation $`\varphi ^{\prime \prime }+(1/r)\varphi ^{}=V^{}`$, with appropriate boundary conditions, which determines the 2d bounce langer . In the absence of forcing, which is the case for homogeneous nucleation, the system will randomly probe phase space until it eventually “finds” the saddle point. The presence of quench-induced large-amplitude fluctuations in the field will drastically accelerate its decay. In figure 18 we show the evolution of the order parameter $`\varphi _{\mathrm{ave}}(t)`$ as a function of time for several values of asymmetry, $`1.518\alpha 1.56`$ for $`T=0.22`$. Not surprisingly, as $`\alpha \alpha _c=1.5`$, the field remains longer in the metastable state, since the nucleation energy barrier $`E_b\mathrm{}`$ at $`\alpha _c`$. However, a quick glance at the time axis shows the fast decay time-scale, of order $`10^{12}`$. For comparison, for $`1.52\alpha 1.56`$, homogeneous nucleation would predict nucleation time-scales of order $`10^{56}\tau _{\mathrm{HN}}10^{21}`$ (in dimensionless units), respectively. \[The nucleation barriers are, $`E_b(\alpha =1.158)=33.67`$ and $`E_b(\alpha =1.56)=13.10`$.\] Clearly, while for smaller asymmetries $`\varphi _{\mathrm{ave}}(t)`$ displays similar oscillatory behavior to the SDW case before transitioning to the global minimum, as $`\alpha `$ is increased (asymmetry in potential is increased) the number of oscillations decreases sharply.

In figure 19 we show the log-log plot of the ensemble-averaged nucleation time-scales as a function of the nucleation barrier, or critical droplet free energy, $`E_b/T`$ for different temperatures, $`T=0.18,0.20,`$ and $`0.22`$. The error bars are computed from the dispersion of the measured time scales within the ensemble. The best fit is a power law, $`\tau _{RN}(E_b/T)^B`$, with $`B=2.440\pm 0.008`$ for $`T=0.20`$ and $`T=0.22`$, and $`B=3.36\pm 0.04`$ for $`T=0.18`$. This simple power law holds for the same range of temperatures where we have observed the synchronous emergence of oscillons. \[The range is similar for both SDW and ADW.\] It is not surprising that the exponent $`B`$ is larger for $`T=0.18`$, where the synchronous emergence of oscillons is marginal: at lower temperatures, the oscillations of the field’s zero mode have lower amplitude and oscillons become rarer. In fact, as we remarked before, we have not observed any for $`T<0.16`$ for the SDW. For lower temperatures, we should thus expect a smooth transition into the exponential time-scales of homogeneous nucleation. We intend to analyse this transition and to obtain an analytical explanation for the power law in a forthcoming publication. As a first step in this direction, we present below what we believe is the mechanism by which the transition completes for different nucleation barriers.

First, note that the error bars in figure 19 increase with nucleation barrier. This spread in time scales is related to the dynamics of the system. For $`\alpha \alpha _c`$ the radius of the nucleation bubble diverges, $`R_b\mathrm{}`$. Clearly, due to the disparity between the power law and the exponential time scales, the mechanism by which the system approaches the global minimum is not through a random search in configuration space as is the case for homogeneous nucleation. Instead, we argue that oscillons will act as seeds for the nucleation of a critical fluctuation. The way in which this happens depends on the nucleation barrier. If $`\alpha \alpha _c`$, the critical nucleus will have a radius much larger than a typical oscillon \[cf. figure 14\]; it will appear as two or more oscillons coalesce. As $`\alpha `$ is increased further, the radius of the critical nucleus decreases, approaching that of an oscillon. In this case, a single oscillon may become the critical nucleus and promote the decay of the metastable state. This explains the small number of oscillations on $`\varphi _{\mathrm{ave}}(t)`$ as $`\alpha `$ is increased \[cf. figure 18\]. To corroborate our argument, in figure 20 we contrast the critical nucleation radius with that of oscillons for different values of energy barrier at $`T=0.22`$. The average oscillon radius \[cf. figure 14\] is $`R_{\mathrm{ave}}=4.30\pm 0.4`$, and is identified by the dashed line, with the dotted lines denoting its dispersion. The critical nucleus radius, $`R_b`$, is equal to $`R_{\mathrm{ave}}`$ for $`\alpha 1.575\pm 0.0075`$, or $`E_b/T=53.9\pm 0.4`$. Thus, for $`\alpha 1.568`$ one cannot distinguish between a single oscillon and a critical bubble. Since, as we have seen, oscillons appear very early on, the decay happens quite fast. For $`1.5<\alpha <1.568`$, approximately, the critical nucleus is the result of the coalescence of two or more oscillons. They diffuse through the lattice and scatter, forming bound states, somewhat as in kink-antikink breathers in 1d field theory campbell ; hanggi . We emphasize that these values are only approximate; they do not include the increase in coalescence rates due to diffusive motions on the lattice or an enhancement on the unstable growth of oscillons due to perturbations. We expect that the one-oscillon decay will ensue for smaller values of $`\alpha `$. In website  the reader can see a few representative simulations.

Clearly, these arguments can be refined. In particular, the details of oscillon coalescence remain unexplored. However, they encapsulate the basic mechanism by which metastable states decay via resonant nucleation, as seen in numerical simulations.

## VI Summary and Outlook

We have investigated the nonequilibrium dynamics of a double-well model with a 2d scalar order parameter. The system was initially prepared in thermal equilibrium at a fixed temperature in a anharmonic single-well potential, $`V(\varphi )=\frac{1}{2}\varphi ^2+\frac{1}{8}\varphi ^4.`$ We then turned on a cubic interaction with tunable strength so that the resulting double-well potential can be symmetric (SDW) or asymmetric (ADW). The turning-on of the interaction can be interpreted as a “quench,” as it occurs within a time-scale faster than any in the system. The system is thus tossed into a double-well potential and, through its nonlinear interactions, redistributes its energy until it reaches equipartition. For the SDW case, we have observed that this approach to equilibrium occurs as a two-step energy cascade: quickly at first, as the zero-mode relinquishes its energy to nearby (in $`k`$-space) long-wavelength modes, and slower later, as the energy is partitioned with the short-wavelength modes. Although we have adopted the language of nonequilibrium field theory, we remark that one could equally well treat the system as a microcanonical ensemble with fixed energy $`E`$.

We showed that, for a range of initial temperatures, the initial stages of the dynamics are dominated by damped oscillations of the field’s zero mode, which, via parametric resonance, excite the modes associated with coherent structures named oscillons. We also showed that the initial emergence of oscillons occurs at the breakdown of the homogeneous Hartree approximation. Furthermore, this emergence occurs in synchrony, characterizing a global emergence phenomenon. As oscillons are long-lived coherent structures, they act as bottlenecks for equipartition. The observed synchronous emergence dies away with the damping of the zero-mode oscillations. Oscillons then appear randomly throughout the lattice, until such a time when there isn’t enough energy in long-wavelength modes.

We then investigated the dynamics in ADWs, where the final equilibrium state is at the global minimum of the potential. The system was thermalized as for the SDW, in its metastable well. Also as before, the quenching induces damped oscillations of the field’s zero mode and oscillons emerge. We observed that, approximately for the same range of temperatures where oscillons emerged in the SDW case, the decay rate was greatly enhanced, when compared to homogeneous nucleation results. In fact, we obtained an excellent fit to a power law decay time scale, $`\tau _{\mathrm{RN}}[E_b/T]^B`$, where $`E_b`$ is the energy of the critical nucleus or bounce and $`B`$ is a numerical exponent, $`2.44B3.36`$, for the temperatures we investigated. We argued that there were two main mechanisms for what we called resonant nucleation: for small asymmetries, the critical nucleus appears as two or more oscillons coalesce. For larger asymmetries, a single oscillon grows to become the critical nucleus, promoting the decay of the metastable state.

The greatly accelerated decay is due to the resonant oscillations of the zero mode, induced by the quenching. It is thus plausible that such behavior may be observed in several related systems with a scalar order parameter. The emergence of oscillons in granular materials appears due to the sinusoidal vibrations of their container umbanhowar ; this is mimicked here by the oscillation of the field’s zero mode, which corresponds to the periodic oscillation of the whole lattice. As for granular materials, patterns emerge above a certain amplitude. We do not claim to have modelled the behavior of granular materials here, but believe to have captured some of the essential physics. Also, since the model here falls in the Ising universality class, we expect that similar qualitative behavior should occur, for example, in ferromagnets and possibly binary liquids and metal alloys, if long wavelength oscillations can be induced by the quenching process.

###### Acknowledgements.

MG was support in part by a National Science foundation grant PHY-0099543. We thank J. D. Gunton and J. A. Krumhansl for their comments and insights.


format: text

text:
# On the angular momentum of the neutron star

## 1

The Crab pulsar NP0532 with the smallest pulsation period .033 sec, if it is connected with 1054 year SN, is also one of the youngest pulsars, neutron stars (n.s.), remnants of SN explosion. If the pulsation period NP0532 is the rotation period  then n.s. with rotation period $`T=0.033`$ sec, mass $`MM_{{\scriptscriptstyle }}`$ and radius $`R10^6`$ cm has rotational momentum less than one-tenth the solar momentum. Evidently or n.s. lost his momentum at birth or the momentum loss took place during n.s. life. In , it is claimed that the large momentum and mass loss is possible during the catastrophic hydrodynamical instability. We reconsider the problem and show that momentum loss is small and mass loss is negligible .

## 2

We assume that the mass ejection of the rotating and contracting star takes place at the star equator. Then for the momentum loss we have expression

$$dL=\alpha ^2R^2\omega dM,$$

where $`\omega `$ is the rotational angular velocity, $`R`$ is the equatorial radius, and we assume that the mass element leaves the star at distance $`\alpha R>R,`$ due to e.g. magnetic forces confining the matter.
Moreover we assume that all the way during the contraction, the injection condition

$$\omega ^2=GM/R^3$$

holds, where $`G`$ is the gravitational constant. Also we write the angular momentum relative the rotation axis in the form

$$L=kMR^2\omega ,$$

where $`k`$ is structural parameter depending on density distribution inside the star, and get finally:

$$\frac{M}{M_0}=\left(\frac{R}{R_0}\right)^\beta ,$$

$$\frac{L}{L_0}=\left(\frac{R}{R_0}\right)^{3\beta /2+1/2},$$

$$\frac{\omega }{\omega _0}=\left(\frac{R}{R_0}\right)^{\beta /23/2},$$

$$\beta =\frac{k}{2\alpha ^23k}.$$

## 3

The values of $`k`$ were calculated for spherical configurations in , and for configurations at the state of rotational instability in . In all cases $`k<1`$; for instance, for $`n=3`$ polytrope case (which is likely the case for star losing his gravitational stability), at $`\alpha 1`$ and taking into account the small deviations from the spherical symmetry, we have $`k=0.038,\beta =0.02.`$

More particularly for a star with $`M\mathrm{2\hspace{0.17em}10}^{33}`$ g, and $`L=\mathrm{4\hspace{0.17em}10}^{48}`$ erg sec, the condition of rotational instability commences first at $`R_0=\mathrm{2.1\hspace{0.17em}10}^7`$ cm, $`\omega =\mathrm{1.2\hspace{0.17em}10}^2\text{sec}^1`$. After contracting to state of n.s. with $`R=\mathrm{1.2\hspace{0.17em}10}^6`$ cm, we get (assuming $`\alpha =1`$):

$$\omega _{n.s.}=\mathrm{8.8\hspace{0.17em}10}^3\text{sec}^1,T_{n.s.}=\mathrm{7\hspace{0.17em}10}^4\text{sec}.$$

## 4

For n.s. formed from a star with mass $`M_0`$, $`\omega _{n.s.}\omega _0R_0^{3/2}M_0^{1/2}`$, that is the angular velocity of ”new-born” n.s. increases with increasing $`M_0`$ and does not depends on $`L_0`$. Also $`M_{n.s.}M_0`$ because the mass loss is negligible in all cases. Taking into account the magnetic forces ($`\alpha >1`$) leads to even lesser values of the momentum and mass loss.
We conclude that the rotational instability can not be an effective mechanism of the momentum and mass loss at n.s. birth, see also .
Therefore the neutron star (pulsar) is formed with a strong rotation and the momentum loss (spin-down) takes place during the pulsar’s life. Of course we do not take into account the possibility of other mechanisms of momentum and mass loss due to e.g. explosions or any other mechanisms  due to uncertainty about their efficiency.
As it is known , the different mechanisms of pulsar spin-down lead to different dependences of $`dT/dt`$ on time $`t`$. If we assume for evaluation purposes that $`dT/dt=cT`$, then we get $`c=10^{10}\text{sec}^1`$, which is close to value observed at the present.
Shemakha Astrophysical observatory, Jan 1970   Z.F. Seidov

## 5 References

1. T. Gold 1968, 1968Nature..218..731G
2. O.Gusejnov, F. Kasumov 1969, 1969ATsir…539..2G
ATsir = Astronomicheskij Tsirkulyar, Moscow, AN SSSR
3. Z. Seidov 1969, 1969Afz…5..503S
4. R. James 1964, 1964ApJ…140..552J
5. L. Auer, N. Woolf 1965, 1965ApJ..142..182A
6. V. Porfiryev 1969, 1969AZh…46..560P
7. I. Shklovsky 1969, 1969Atsir..494..2S
8. T. Gold 1969, 1969Natur..221…25G
9. I. Gunn, J. Ostriker 1969, 1969Natur..221..454G
10. V. Safronov 1951, 1951AZh ..28..244S
11. Z. Seidov, 1970, 1970ATsir..567..3S
published 01 Jun 1970


text:
# Effects of forcing in three dimensional turbulent flows

## Abstract

We present the results of a numerical investigation of three-dimensional homogeneous and isotropic turbulence, stirred by a random forcing with a power law spectrum, $`E_f(k)k^{3y}`$. Numerical simulations are performed at different resolutions up to $`512^3`$. We show that at varying the spectrum slope $`y`$, small-scale turbulent fluctuations change from a forcing independent to a forcing dominated statistics. We argue that the critical value separating the two behaviours, in three dimensions, is $`y_c=4`$. When the statistics is forcing dominated, for $`y<y_c`$, we find dimensional scaling, i.e. intermittency is vanishingly small. On the other hand, for $`y>y_c`$, we find the same anomalous scaling measured in flows forced only at large scales. We connect these results with the issue of universality in turbulent flows.

The effects of both external forcing mechanisms and boundary conditions on small-scale turbulent fluctuations have been the subject of many theoretical, numerical and experimental studies frisch ; SA97 . The 1941 theory of Kolmogorov frisch  is based on the assumption of local isotropy and homogeneity, that is any turbulent flow, independently on the injection mechanism, recovers universal statistical properties, for scales small enough (and far from the boundaries). Indeed, experiments and numerical simulations give strong indications that Eulerian and Lagrangian isotropic/anisotropic small-scales velocity statistics are pretty independent of the large-scale forcing mechanisms SW02 ; GFN02 ; boden ; BCTT03 ; BBCLTV03 . Still, we lack a firm understanding for these evidences. From the theoretical point of view, precious hints arise from linear problems, like passive scalar or passively advected magnetic fields. For the class of Kraichnan models kr94 , anomalous scaling has been shown to be associated to statistically stationary solutions of the unforced equations for correlation functions FGV01 . Scaling exponents are consequently universal with respect of the injection mechanisms. Concerning non-linear problems, as the Navier-Stokes case, analytical results have been often pursued by means of the Renormalization Group (RG) FNS77 ; FF78 . In the RG framework, turbulence is stirred at all scales by a self-similar Gaussian field, with zero mean and white-noise in time. The two-point correlation function in Fourier space is given by

$`f_i(𝒌,t)f_j(𝒌^{},t^{})`$ (1)
$`D_0k^{4d}(k_0^2+k^2)^{y/2}P_{ij}(𝒌)\delta (𝒌+𝒌^{})\delta (tt^{}).`$
Here $`1/k_0L`$ is the largest length in the system (infrared cut-off), $`D_0`$ is the forcing intensity, $`P_{ij}(𝒌)`$ is the projector assuring incompressibility and $`d`$ is the spatial dimension (always assumed to be $`d=3`$ hereafter). The influence of the stirring mechanism at small scales is governed by the value of the slope $`y`$. We go from a situation when the forcing has a strong input at all scales, $`y0`$ originally investigated in FNS77 , to a quasi large-scale forcing when $`y\mathrm{}`$. Renormalization Group calculations, based on a $`y`$-expansion, predict a power-law energy spectrum $`E(k)k^{12y/3}`$, in the domain $`\eta k^1L`$, where $`\eta `$ is the viscous scale of the system, and for $`y1`$. Notice that the Kolmogorov value, $`E(k)k^{5/3}`$, describing experimental turbulent flows stirred by a large-scale forcing, is obtained for $`y=4`$, i.e. quite far from the perturbative region where the RG calculations are under control. The Kolmogorov spectrum can be obtained, however, by means of a simple dimensional analysis, still within the same framework DDM79 . Extension of the RG formalism to finite $`y`$ values, up to $`y=4`$, have been attempted with different kind of approximation YO86 ; AAKV03  altough in a range where convergence of the RG expansion is not granted anymore EY94 . As for the numerical simulations, in SMP98  the problem has been investigated for various $`y`$ values, and it has been shown that, for $`y=4`$, results are in good agreement with the picture of large-scale forced turbulence, while for $`y<4`$ the situation becomes less clear. However, the low numerical resolution used in SMP98  makes these results far from being conclusive.

Beside the issue connected to the RG approach, there exists a whole set of interesting questions concerning turbulent flows with a power-law forcing. To which extent small-scale fluctuations are sensitive to the injection mechanism? Does it exist a critical value $`y_c`$ separating different regimes? Can we observe anomalous scaling in the forcing dominated case $`y0`$? Let us suppose, for example, that it exists a finite $`y_c`$, beyond which small-scale statistics is forcing independent: this would rule out any attempt to control intermittency analytically, by means of a perturbative approach which starts from forcing dominated turbulent solutions at $`yy_c`$.

Hints on the problem can also come from the study of the one-dimensional Burgers equation in presence of a power law forcing. In HJ97 , a numerical study was presented showing that there is a critical value of the forcing slope, such that the velocity field passes from the usual bifractal statistics (observed in large-scale forced Burgers flows), to a statistics influenced from the forcing. Surprisingly, also in the forcing dominated regime, a non-trivial (multifractal?) scaling was observed in HJ97 . A rigorous understanding of the mechanism leading to this result is, however, still missing: we will come back to this point later on, after having presented our results.

In the sequel, we address the problem of the small-scale statistical properties of three dimensional turbulent flows at varying the parameter $`y`$. We will show that at crossing $`y_c=4`$ there exist two regimes for the velocity field statistics with well defined and different scaling properties.

We solved the Navier-Stokes equations with a second-order hyperviscous dissipative term $`\nu \mathrm{\Delta }^2`$. Temporal integration has been carried over for about $`20÷30`$ large-eddy turnover times. We performed variuos experiments, at resolutions $`128^3`$, $`256^3`$ and $`512^3`$, corresponding to a maximum Taylor’s Reynolds number equal to $`Re_\lambda =220`$ for the $`512^3`$ run. As for the stirring force, we specialized in the two following cases, one for each regime: the first with $`y=3.5<y_c`$, the second with $`y=6>y_c`$. We also show some results obtained with an analityical large-scale forcing, i.e. the equivalent of $`y\mathrm{}`$. The range of the forcing, in Fourier space, extends down to the maximum resolved wavenumber. As we are always confined in a finite box, we neglect here possible subtle behaviours due to infrared divergences in the injection mechanism.

We start considering what happens to the system when the slope of the forcing is changed. It is instructive to consider the equation for the energy flux through the wavenumber $`k`$: $`\mathrm{\Pi }(k)\pi k^2\mathrm{}[(𝒌\widehat{𝒗}(𝒒))(\widehat{𝒗}(𝒌)\widehat{𝒗}(𝒑))+(𝒌\widehat{𝒗}(𝒑))(\widehat{𝒗}(𝒌)\widehat{𝒗}(𝒒))]\text{d}𝒑\text{d}𝒒`$, where the three wavevectors satisfy $`𝒌+𝒑+𝒒=0`$, and the symbol $`\mathrm{}`$ stands for the imaginary part. Such equation is equivalent to the Kármán-Howarth equation in physical space; it states that in a stationary, isotropic and homogeneous flow, the contribution to the energy flux $`\mathrm{\Pi }(k)`$ due to the non-linear terms balances the total energy input from the injection mechanism (see my ; frisch ):

$$\mathrm{\Pi }(k)_{k_0<|𝒌|<k}\mathrm{}(𝒇(𝒌)𝒗(𝒌))\text{d}𝒌,$$
(2)
where $`k_01/L`$, the symbol $`\mathrm{}`$ stands for the real part and where we have neglected dissipative effects. For the special class of forcings (1), the rhs of (2) can be further simplified to: $`\mathrm{\Pi }(k)_{k_0<|𝒌|<k}|𝒇(𝒌)|^2\text{d}𝒌`$. From (1) the forcing spectrum is $`E_f(k)=|𝒇(𝒌,t)|^2k^{3y}`$. It follows that for $`y4`$, the energy flux is constant in Fourier space for $`kL1`$ (up to logarithmic corrections for $`y=4`$). In other words, the energy injection is dominated by the small wave-number region in the integral (2). In this case we expect to be very close to the typical experimental situation of turbulence with a large-scale, analytical forcing: energy is transferred down-scale via  an intermittent cascade. Coherently, the third order longitudinal structure function, $`S^{(3)}(r)\left[(𝒗(𝒙+𝐫)𝒗(𝒙))\widehat{𝒓}\right]^3`$, follows a linear behaviour in $`r`$ as predicted by the $`4/5`$ law frisch . For $`y<4`$, the energy flux no longer saturates to a constant value as a function of $`k`$. The integral in (2) now becomes ultraviolet dominated. The direct input of energy from the forcing mechanism affects inertial range statistics in a self-similar way, down to the smallest scales where dissipative terms start to be important. In this situation, we get for the energy flux $`\mathrm{\Pi }(k)k^{4y}`$, with a constant prefactor which depends on the ultraviolet cut-off. The corresponding scaling behaviour for the third order structure function is now given by $`S^{(3)}(r)r^{y3}`$.

What about higher order statistics? One is tempted to guess that for $`y4`$ the fluctuations induced by the injection mechanism are always sub-leading, anomalous scaling being the result of the cascade mechanism driven by the non-linear terms of the equations of motion. If so, for $`y>4`$ we should fall in the same class of “universality” of turbulence generated with large-scale forcing, i.e. small-scale velocity fluctuations should be universal. Therefore, as far as $`y>4`$, longitudinal structure functions should scale as:

$$S^n(r)\left[(𝒗(𝒙+𝐫)𝒗(𝒙))\widehat{𝒓}\right]^nr^{\zeta _{\mathrm{}}^{(n)}}.$$
(3)
In (3), we have denoted with $`\zeta _{\mathrm{}}^n`$ the scaling exponents measured with a smooth large-scale forcing.

On the other hand, for $`y<4`$, energy is directly injected in the inertial range. Here, a dimensional matching with the forcing gives a scaling behaviour which is always leading with respect to what predicted in the $`y>4`$ range (3). We expect now that anomalous scaling disappears, everything being dominated by the Gaussian energy input at all scales. By the simple dimensional argument connecting the scaling of structure functions to that of the external forcing, for the range $`y<4`$ we have:

$$S^{(n)}(r)r^{\zeta _y^{(n)}}\text{with}\zeta _y^{(n)}=\frac{n}{3}(y3).$$
(4)
In Fig. 1 we show the sixth order structure function, $`S^{(6)}(r)`$ for the two cases $`y=3.5`$ and $`y=6`$, compensated with the dimensional prediction given by the matching with the forcing (4).

As it is clear, only for $`y=3.5`$ the statistics follows the forcing injection obtaining a nice compensation. On the other hand, for $`y=6`$ the statistics is much closer to what usually measured with an analytical large-scale forcing. This is quantitatively confirmed by the inset picture, where we plot the logarithmic derivatives of $`S^{(6)}(r)`$ vs. $`S^{(3)}(r)`$, for the two cases $`y=3.5`$ and $`y=6`$, together with the results of the simulation with a large-scale, analytical forcing corresponding to $`y\mathrm{}`$. Here the local slopes for $`y=6`$ and $`y\mathrm{}`$ fluctuate around the same value, compatible with those reported in literature SA97 ; GFN02 , while the local slope for $`y=3.5`$ is different and tends to the expected dimensional value. Values of all scaling exponents obtained in the simulations are summarized in Table 1. Let us notice that the measured exponents for the $`y=3.5`$ case are very close to the non-intermittent, dimensional prediction. Only for high order moments (i.e. $`n=6`$) there is a small deviation from the expected value $`\zeta ^{(6)}/\zeta ^{(3)}=2`$.

To quantify the level of intermittency at changing the scale, we also plot the Probability Density Function (PDF) of velocity increments at different scales, normalized to have unit variance. Figure 2 shows the PDFs, in the case $`y=6`$, for three different separations $`r_1=34\eta `$ and $`r_2=74\eta `$ in the inertial range, and $`r_3=114\eta `$ in the energy containing range. The three curves have larger tails than a Gaussian distribution and have an intermittent behaviour, i.e. they cannot be superposed. In Fig. 3 we show the PDFs, for the case $`y=3.5`$, at the same separations $`(r_1,r_2,r_3)`$. Now, the three curves are almost indistinguishable, and show a very good rescaling, a signature of the absence of intermittent effects. Only for negative increments, a very tiny discrepancy is measured. It is hard to say whether this is a robust effect or a spurious Reynolds dependent phenomenon. We will come back to this issue later in the conclusions.

A dramatic difference at crossing the $`y=y_c`$ value is also observed in the energy dissipation statistics. For both cases $`y=3.5`$ and $`y=6`$, we measured the PDF’s of the coarse-grained energy dissipation $`\epsilon _r(𝒙)=V_r^1_{V_r(𝒙)}\epsilon (𝒙+𝐫)𝑑𝐫`$, where $`\epsilon `$ is the rate of dissipation for unit volume, and the volume $`V_r(𝒙)`$ is centered at $`𝒙`$ and has characteristic length scale $`rL`$. In Fig. 4 we compare the PDF’s $`𝒫(\epsilon _r)`$ at the scale $`r=8\eta `$. Here the results are even more impressive as the shape change is particularly strong.

In this letter, we have presented clear evidences that turbulent small-scale fluctuations, in presence of a direct injection of energy at all scales, undergo a sharp transition for $`y_c=4`$. The first regime, for $`y>y_c`$, is mainly forcing independent: small-scale fluctuations develop anomalous scaling in agreement with what observed in experiments and/or numerics obtained with a large-scale forcing. This is a stringent test of turbulence universality: even if directly affected by the injection of energy, small-scale fluctuations show a robust behaviour. Things change abruptly at the critical value of $`y_c`$, where the direct injection of energy becomes the dominant effect in the inertial range. In this second regime, corresponding to $`y<y_c`$, small-scale fluctuations get closer and closer to a Gaussian statistics and intermittency disappear note .

Before concluding, we discuss two possible mechanisms which could partially disprove the last statement. First, even when $`y<y_c`$, we may imagine that the intermittent energy cascade dominating the statistics for $`y>y_c`$ might show up. For example, we may have that for high order moments the forcing dominated solutions become subleading with respect to those associated to the cascading mechanism. In such case, the loss of rescaling in the PDF’s tail of Fig. 2, may be due to the survival of these rare anomalous fluctuations. Second, even more complex is the scenario proposed in HJ97 , where the possibility to have a forcing dependent multiscaling statistics, when $`y<y_c`$, is conceived. This is not the case for the linear Kraichnan models FGV01 , where forcing dependent solutions are always dimensionally scaling. The main difference is that, in the Navier-Stokes problem, the hierarchy of equations for correlation function is unclosed: one cannot solve it for a single order independently of all the others. In NS, being the low order moments always forcing dominated for $`y<y_c`$, one may observe some forcing dependency also on high order fluctuations, via their coupling with low order correlation functions. This may be a possible explanation of the multifractal behaviour observed in Burgers equation HJ97 .

We are grateful to N. V. Antonov, J. Bec, A. Celani, U. Frisch, M. Sbragaglia, and M. Vergassola for useful discussions. This research was supported in part by the EU under the Grants No. HPRN-CT 2000-00162 “Non Ideal Turbulence” and by the INFM (Iniziativa di Calcolo Parallelo).