text:
# A Statistical Mechanics Model of Isotropic Turbulence Well-Defined within the Context of the ϵ Expansion

## 1 Introduction

Turbulence remains an outstanding problem in classical mechanics. Early on, several self-consistent or closure-based statistical models were introduced, with perhaps Kraichnan’s Direct Interaction Approximation the most well known . Within the statistical mechanics literature, attention has focused upon the random force model that Dominicis and Martin  generalized from Forster, Nelson, and Stephen’s model of a randomly stirred equilibrium fluid . This model has received quite a bit of attention, and its transport properties have been further examined by Yakhot and Orszag and Avellaneda and Majda and co-workers . Various treatments, up to and including field-theoretic $`ϵ`$ expansions, have been applied to the random force model of turbulence. As pointed out by Eyink, however, even the $`ϵ`$ expansion does not lead to a controlled calculation in the random force model because the term representing the effects of high Reynolds numbers must be calculated to arbitrarily high order in perturbation theory, even for an $`O(ϵ)`$ calculation . More recently, a variety of novel renormalization group techniques has been applied to the problem of Navier-Stokes turbulence .

The experimentally observed scaling behavior of turbulent energy dissipation, often called the Kolmogorov energy cascade, suggests that there may be a strong analogy between critical phenomena and turbulence. Indeed, the search for such an analogy motivated much of the field-theoretic work on fluid turbulence . Just as there are wild fluctuations in particle density near a critical point, so to are there large fluctuations in the fluid velocity at high Reynolds number turbulence. This similarity suggests that the effects of turbulence may be modeled by a random, velocity-dependent force, just as the effects of critical fluctuations can be modeled by a random, density-dependent force in the standard $`\varphi ^4`$ model. While there are likely random forces that are independent of the fluid velocity in the context of turbulence, there are also very likely random forces that are velocity dependent. Thus, the turbulent force should really depend on the gradient of the velocity, as it is not simply large velocities that lead to turbulence, but rather regions of high gradient, such as walls or boundaries, that lead to turbulent behavior. Such boundary roughness is one mechanism for breaking the symmetry from laminar to turbulence fluid flow.

Indeed, random boundary roughness along the walls generates velocity-gradient-dependent forces in the bulk that are quenched in time. Moreover, in practical, engineering-type calculations, the turbulent forces are often related to velocity gradients by a constitutive relation that contains a “turbulent viscosity” parameter, the simplest of these relations being linear . Following this line of reasoning, we introduce a new statistical mechanics model for turbulence that renormalizes the effect of turbulent stresses into a velocity-gradient-dependent term. This model will turn out to be well-defined within the $`ϵ`$ expansion. That is, the renormalization group theory takes into account all physical effects of this model, consistently to order $`ϵ`$. Due to its similarity with turbulent viscosity type models, this approach may lead to a closer connection with practical turbulence calculations.

We introduce our velocity-gradient-dependent random force model in Sec. 2. The problem is cast in the framework of time-dependent field theory in Sec. 3. Renormalization group flow equations are also calculated in this section, and three appendices provide details of the renormalization group calculations. The behavior of these flow equations in two and three dimensions, and the predictions for the Kolmogorov energy cascade, are described in Sec. 4. Nontrivial intermittency corrections to the single-time structure functions are calculated by an operator product expansion in Sec. 5. We conclude in Sec. 6.

## 2 Velocity-Gradient-Dependant Random Force Model

Our goal is to write a form of the Navier-Stokes equation that contains a random piece, the random piece representing the statistical effects of the turbulence. The Navier-Stokes equation with a random force is

$$_tv_i+\underset{k}{}_{ik}_jv_j_jv_k=\nu ^2v_i+f_i,$$
(1)
where $`\nu =\mu /\rho `$ is the kinematic viscosity, and $`f_i=_k\mathrm{\Pi }_{ik}(F_k_kP)/\rho `$ is the total body force on the fluid. The presence of the projection operator $`\widehat{\mathrm{\Pi }}_{ik}(𝐤)=\delta _{ik}k_ik_k/k^2`$ in these formulas ensures that the incompressibility condition $`𝐯=0`$ is maintained . The Fourier transform is defined by $`\widehat{f}(𝐤)=𝑑𝐱f(𝐱)\mathrm{exp}(i𝐤𝐱)`$.

We choose the random force to depend on the gradient of the velocity:

$$f_i(𝐱,t)=\gamma _{ijk}(𝐱)_jv_k(𝐱,t).$$
(2)
We assume that $`\gamma _{ijk}`$ is symmetric under exchange of $`j`$ and $`k`$, so that the force looks like a turbulent stress. We average over the statistics of this force using a field-theoretic representation of the Navier-Stokes equation. We choose the correlation function to be

$$\widehat{\gamma }_{ijk}(𝐤_1)\widehat{\gamma }_{lmn}(𝐤_2)=(2\pi )^d\delta (𝐤_1+𝐤_2)|𝐤_1|^yD_{ijk}^{lmn},$$
(3)
where

$`D_{ijk}^{lmn}=D_\alpha (\delta _{il}\delta _{jm}\delta _{kn}+\delta _{in}\delta _{jm}\delta _{kl}`$
$`+\delta _{il}\delta _{jn}\delta _{km}+\delta _{im}\delta _{jl}\delta _{kn}+\delta _{im}\delta _{jn}\delta _{kl}+\delta _{in}\delta _{jl}\delta _{km})`$
$`+D_\beta (\delta _{ij}\delta _{lm}\delta _{kn}+\delta _{ik}\delta _{jm}\delta _{ln}+\delta _{ik}\delta _{jn}\delta _{lm}+\delta _{ij}\delta _{km}\delta _{ln}).`$
(4)
Initially, we treat this as a mathematical problem, taking $`y`$ to be arbitrary. Later, we determine $`y`$ by requiring that the transport properties of turbulence are reproduced. This scaling form of the correlation function, Eq. (3), applies only in the inertial, Kolmogorov regime, for wavevectors below an upper cutoff related to the inverse of the dissipation length scale and above a lower cutoff related to the inverse of the so-called integral length scale. It is this Kolmogorov scaling regime that is of interest in the present work. Given the form of Eq. 2, the parameters $`\sqrt{D}_\alpha `$ and $`\sqrt{D}_\beta `$ can be viewed as modeling gradients of the turbulent viscosity.

## 3 Renormalization Group Calculations

We write the Navier-Stokes equation in field-theoretic form so that the renormalization group can be applied systematically within the $`ϵ`$ expansion . Within the field-theoretic formalism, any observable can be calculated. The average velocity, for example, is given by $`v_i(𝐱,t)=b_i(𝐱,t)`$, where the average over the $`b`$ field is taken with respect to the weight $`\mathrm{exp}(S)`$. Using Eqs. 14, we arrive at the following action:

$`S`$ $`=`$ $`{\displaystyle _𝐤}{\displaystyle 𝑑t\widehat{\overline{b}}_i(𝐤,t)[_t+\nu k^2+\delta (t)]\widehat{b}_i(𝐤,t)}`$
$`+i\lambda {\displaystyle _{𝐤_1𝐤_2𝐤_3}}{\displaystyle 𝑑t(2\pi )^d\delta (𝐤_1+𝐤_2+𝐤_3)}`$
$`\times k_{1_j}\widehat{\overline{b}}_k^{}(𝐤_1,t)\widehat{b}_k^{}(𝐤_2,t)\widehat{b}_j^{}(𝐤_3,t)`$
$`+{\displaystyle \frac{1}{2}}{\displaystyle _{𝐤_1𝐤_2𝐤_3𝐤_4}}{\displaystyle 𝑑t_1𝑑t_2(2\pi )^d\delta (𝐤_1+𝐤_2+𝐤_3+𝐤_4)}`$
$`\times 𝐤_{2}^{}{}_{j}{}^{}𝐤_{4}^{}{}_{m}{}^{}D_{ijk}^{lmn}|𝐤_1+𝐤_2|^y`$
$`\times \widehat{\overline{b}}_i^{}(𝐤_1,t_1)\widehat{b}_k^{}(𝐤_2,t_1)\widehat{\overline{b}}_l^{}(𝐤_3,t_2)\widehat{b}_n^{}(𝐤_4,t_2).`$
The notation $`_𝐤`$ stands for $`𝑑𝐤/(2\pi )^d`$, the integrals over time are from $`t=0`$ to some large time $`t=t_f`$, and the summation convention is implied. This action is written in terms of the divergence-free part of the field, $`\widehat{b}_i^{}(𝐤)=_k\widehat{\mathrm{\Pi }}_{ik}(𝐤)\widehat{b}_k(𝐤)`$. We have used the Feynman gauge, adding in a curl-free component in the quadratic terms to make later calculations easier. Initially $`\lambda =1`$. We have used the replica trick  to incorporate the statistical disorder, but have suppressed these details since they do not enter in a one-loop calculation.

We now apply renormalization group theory to this action. It is important to note that the fields must all have zero average value before the renormalization group is applied , otherwise it would not be correct to truncate perturbation theory at any finite order. If there were an average velocity, the action in (LABEL:5) would be different, containing a term of the form $`i𝐤𝐯`$ in the propagator. The vertices in the theory are shown in Fig. 1.

From power counting, the upper critical dimension for this theory is $`d_\mathrm{c}=2+y`$. Note that this upper critical dimension is exactly defined once the model is specified . The deviation of the physical dimension from the upper critical dimension is parameterized by $`ϵ=d_\mathrm{c}d`$. We use the momentum shell procedure, where fields on a shell of differential width $`\mathrm{ln}a=dl`$ are integrated out, $`\mathrm{\Lambda }/a<k<\mathrm{\Lambda }`$. Note that the combination $`dl`$ invariably means a differential on $`l`$; in all other cases, the factor $`d`$ denotes the physical dimension. As usual, we rescale time by the dynamical exponent $`t^{}=a^zt`$ and distance by $`𝐤_{}^{}=a𝐤_{}`$. The $`b`$ field is scaled as $`\widehat{b}^{}(𝐤^{},t^{})=a^{z1d+\alpha }\widehat{b}(𝐤,t)`$. To maintain dimensional consistency, so that the $`b`$ field scales as a velocity, one must set $`\alpha =0`$ . To keep the time derivative in $`S`$ constant , the $`\overline{b}`$ field is scaled as $`\widehat{\overline{b}}^{}(𝐤^{},t^{})=a^{1z\alpha }\widehat{\overline{b}}(𝐤,t)`$. In the loop calculation, we make use of the relation for the reference system averages

$`\widehat{\overline{b}}_i^{}(𝐤_1,t_1)\widehat{b}_j^{}(𝐤_2,t_2)_0`$ $`=`$ $`(2\pi )^d\delta (𝐤_1+𝐤_2)`$
$`\times \widehat{\mathrm{\Pi }}_{ij}(𝐤_2)\widehat{G}_0(𝐤_2,t_2t_1)`$
where $`\widehat{G}_0(𝐤,t)=\mathrm{exp}[\nu k^2t]\mathrm{\Theta }(t)`$, and $`\mathrm{\Theta }(t)=1`$ if $`t>0`$ and 0 otherwise. Note that elimination of modes at one end of the spectrum by perturbation theory is the standard procedure in renormalization group theory . Use of Eq. (LABEL:5aaa) does not imply that the system is somehow Gaussian, as the parameters within the renormalized theory are flowing. The critical properties of the Ising model at the non-Gaussian Wilson-Fisher fixed point, for example, are analyzed in exactly this way . We make use of the rotational averages: $`k_lk_u_\mathrm{\Omega }=\delta _{lu}k^2/d`$, $`k_lk_mk_sk_u_\mathrm{\Omega }=(\delta _{lm}\delta _{us}+\delta _{lu}\delta _{ms}+\delta _{ls}\delta _{mu})k^4/[d(d+2)]`$, and $`k_lk_mk_nk_sk_tk_u_\mathrm{\Omega }=M_{lmn}^{stu}k^6/[d(d+2)(d+4)]`$, where the function $`M_{lmn}^{stu}`$ is equal to all possible couplings of pairs of the arguments:

$`M_{lmn}^{stu}`$ $`=`$ $`\delta _{lt}\delta _{ms}\delta _{nu}+\delta _{lt}\delta _{mn}\delta _{su}+\delta _{lt}\delta _{mu}\delta _{ns}`$ (7)
$`+`$ $`\delta _{lm}\delta _{ts}\delta _{nu}+\delta _{lm}\delta _{nt}\delta _{su}+\delta _{lm}\delta _{tu}\delta _{ns}`$
$`+`$ $`\delta _{ls}\delta _{mt}\delta _{nu}+\delta _{ls}\delta _{mn}\delta _{tu}+\delta _{ls}\delta _{mu}\delta _{nt}`$
$`+`$ $`\delta _{ln}\delta _{mt}\delta _{su}+\delta _{ln}\delta _{ms}\delta _{tu}+\delta _{ln}\delta _{mu}\delta _{st}`$
$`+`$ $`\delta _{lu}\delta _{mt}\delta _{ns}+\delta _{lu}\delta _{mn}\delta _{st}+\delta _{lu}\delta _{ms}\delta _{nt}.`$
The one loop contributions are shown in Fig. 2.

There are 13 diagrams of the form Fig. 2a, detailed calculation of which is described in Appendix A. There are 13 similar diagrams of the form Fig. 2b, detailed calculation of which is described in Appendix B. Finally, there are 100 diagrams of the form Fig. 2c, and a detailed discussion of them is given in Appendix C. These last diagrams require some care in their calculation, as they contribute to the complex tensor structure of $`D_{ijk}^{lmn}`$ in Eq. (LABEL:5). These last diagrams make contributions in the form

$`𝐤_{2}^{}{}_{j}{}^{}𝐤_{4}^{}{}_{v}{}^{}\widehat{\overline{b}}_i^{}(𝐤_1,t_1)\widehat{b}_k^{}(𝐤_2,t_1)\widehat{\overline{b}}_r^{}(𝐤_3,t_2)\widehat{b}_w^{}(𝐤_4,t_2)`$
$`\times M_{lmn}^{stu}D_{ijk}^{lmn}D_{rst}^{uvw}.`$ (8)
There are 15 terms of the form Fig. 2c when the sums over $`l,m,n`$ and $`s,t,u`$ are taken. Each of these terms corresponds to one of the 15 terms in Eq. (7), and each gives a contribution that is exactly of the form in Eq. (LABEL:5). The theory is, therefore, self-consistent in that no terms are generated at this order that are not of the original form, and the symmetry of the original theory is maintained.

The flow equations that result from the one-loop calculation are

$`{\displaystyle \frac{d\mathrm{ln}\nu }{dl}}`$ $`=`$ $`z2+[(d^2+2d2)D_\alpha +(d2)D_\beta ]`$
$`\times {\displaystyle \frac{yK_d}{\nu ^2d(d+2)}}\mathrm{\Lambda }^ϵ`$
$`{\displaystyle \frac{d\mathrm{ln}\lambda }{dl}}`$ $`=`$ $`\alpha [(d^2+2d2)D_\alpha +(d2)D_\beta ]`$
$`\times {\displaystyle \frac{K_d}{\nu ^2d(d+2)}}\mathrm{\Lambda }^ϵ`$
$`{\displaystyle \frac{d\mathrm{ln}D_\alpha }{dl}}`$ $`=`$ $`2z4+ϵ[2D_\alpha +4D_\beta ]`$
$`\times {\displaystyle \frac{2K_d(d^2+2d2)}{\nu ^2d(d+2)(d+4)}}\mathrm{\Lambda }^ϵ`$
$`{\displaystyle \frac{d\mathrm{ln}D_\beta }{dl}}`$ $`=`$ $`2z4+ϵ[2(d^2+2d2)D_\alpha ^2/D_\beta `$ (9)
$`+(d^3+6d^28)D_\alpha +(d^2+2d8)D_\beta ]`$
$`\times {\displaystyle \frac{2K_d}{\nu ^2d(d+2)(d+4)}}\mathrm{\Lambda }^ϵ.`$
The constant $`K_d=S_d/(2\pi )^d`$, and $`S_d=2\pi ^{d/2}/\mathrm{\Gamma }(d/2)`$. In solving these equations, we set $`\alpha =0`$ for dimensional consistency. In the standard model of turbulence , terms higher order in $`\lambda (l)`$ must be kept in the flow equations. In the present model, $`\lambda (l)`$ flows to zero rapidly, and higher order do not contribute at this level.

## 4 Results and Discussion

In two dimensions, both $`D_\alpha `$ and $`D_\beta `$ are relevant. We find a fixed point of $`D_\alpha ^{}=D_\beta ^{}=2ϵ\nu ^2\mathrm{\Lambda }^ϵ/[3K_2(1+y)]`$. The dynamical exponent is given by $`z=2yϵ/[2(1+y)]`$. The Reynolds number term scales as $`\lambda (l)=\lambda _0e^{ϵl/[2(1+y)]}`$.

In greater than two dimensions, $`D_\beta `$ reaches a fixed point, but $`D_\alpha `$ flows to zero. We find

$`D_\beta ^{}`$ $`=`$ $`{\displaystyle \frac{d(d+2)\nu ^2\mathrm{\Lambda }^ϵ}{2(1+y)(d2)K_d}}ϵ`$
$`z`$ $`=`$ $`2{\displaystyle \frac{yϵ}{2(1+y)}}`$
$`\lambda (l)`$ $`=`$ $`\lambda _0e^{ϵl/[2(1+y)]}.`$ (10)
Interestingly, the dynamical exponent is the same in two and greater dimensions, as is the decay of $`\lambda (l)`$.

A particularly beautiful feature of this theory is that $`\lambda (l)`$ decays exponentially to zero. This property is what makes the theory well defined within the epsilon expansion. If $`\lambda (l)`$ had stayed at unity, the vertex in Fig. 1b could be inserted arbitrarily many times in the loop expansion, and terminating the expansion at any finite order would not be justified by any small parameter. The present calculation, on the other hand, is a controlled expansion in $`ϵ`$ and $`\lambda (l)`$. Note that the quenched random forces, which mimic the effects of, say, wall roughness, break statistical Galilean invariance, and this allows $`\lambda (l)`$ to flow, in contrast to the conventional model with random forces delta-correlated in time. Indeed, our explicit calculation shows that $`\lambda (l)`$ decays exponentially to zero.

We now turn to a calculation of the energy spectrum, defined by

$$E(𝐤)=\frac{(d1)}{2}K_dk^{d1}\widehat{C}_{11}(𝐤),$$
(11)
where the velocity-velocity correlation function is given in the field-theoretic language as

$$(2\pi )^d\delta (𝐤_1+𝐤_2)\widehat{\mathrm{\Pi }}_{ij}(𝐤_1)\widehat{C}_{ij}(𝐤_1)=\widehat{b}_i^{}(𝐤_1,t)\widehat{b}_j^{}(𝐤_2,t).$$
(12)
Under the scaling of time and space that occurs within the renormalization group calculation, this correlation function scales as

$`\widehat{C}_{ij}(𝐤)`$ $`=`$ $`{\displaystyle \frac{a^{2(d+1z\alpha )}}{a^d}}{\displaystyle \frac{\widehat{b}_i^{}{}_{}{}^{}(𝐤_1^{},t^{})\widehat{b}_j^{}{}_{}{}^{}(𝐤_2^{},t^{})}{(2\pi )^d\delta (𝐤_1^{}+𝐤_2^{})\widehat{\mathrm{\Pi }}_{ij}(𝐤_1^{})}}`$ (13)
$`=`$ $`a^{d+22z2\alpha }\widehat{C}_{ij}^{}(𝐤^{}).`$
Making the assumption that $`\widehat{C}_{ij}(𝐤)(\mathrm{const})k^\delta `$, we find from Eq. 13 that $`\delta =d+22z2\alpha `$. The energy spectrum, therefore, scales as

$$E(𝐤)(\mathrm{const})k^{2z3+2\alpha }.$$
(14)
For any isotropic statistical theory of turbulence, then, the dimensional consistency condition of $`\alpha =0`$ enforces a relation between the dynamical exponent and the exponent of the energy cascade. In particular, the Richardson separation law implies $`z=2/3`$, and this result is equivalent to enforcing the Kolmogorov energy cascade: $`E(𝐤)(\mathrm{const})k^{5/3}`$. The relation $`z=2/3`$ implies $`y=3.44152`$ in two dimensions and $`y=4.28849`$ in three dimensions.

To calculate the Kolmogorov constant $`C_\mathrm{K}`$ we introduce a source of randomness into the model:

$$\delta S=\frac{D}{2}\underset{i}{}𝑑t_𝐤\widehat{c}(k)\widehat{\overline{b}}_i^{}(𝐤,t)\widehat{\overline{b}}_i^{}(𝐤,t).$$
(15)
In the range for which scaling occurs, we set $`\widehat{c}(k)k^{3z2d}`$. For later convenience, we also require that $`lim_{r0}c(r)=1`$. The randomness expressed in this source term drives the model away from the trivial solution $`v_i(𝐱,t)0`$. Note that the randomness parameter scales as $`D(l)D`$. This randomness parameter does not contribute to $`\nu `$ since $`\lambda (l)`$ flows to zero, and nothing contributes to $`D(l)`$ at one loop. Since $`D(l)`$ contributes to physical properties at higher loops only through $`\lambda (l)`$, $`D_\alpha (l)`$, and $`D_\beta (l)`$, all of which are small, the effects of $`D`$ are controlled within the $`ϵ`$ expansion. Using the matching Eq. 13, the correlation function is given by $`\widehat{C}_{ij}(𝐤)=[D\mathrm{\Lambda }^{z2}/(2\nu )]k^\delta `$. For fully-developed isotropic turbulence $`z=2/3`$, and using the notation

$$E(𝐤)C_\mathrm{K}\epsilon ^{2/3}k^{5/3},$$
(16)
we find $`C_\mathrm{K}\epsilon ^{2/3}=(d1)K_dD\mathrm{\Lambda }^{4/3}/(4\nu )`$. Similarly, the wavevector-dependent viscosity considered in the fluid mechanics literature  is given by $`\nu (𝐤)=\nu \mathrm{\Lambda }^{4/3}k^{4/3}`$. Using the notation

$$\nu (𝐤)N\epsilon ^{1/3}k^{4/3},$$
(17)
we find $`N\epsilon ^{1/3}=\nu \mathrm{\Lambda }^{4/3}`$. The energy dissipation rate is given by

$$\epsilon =_0^\mathrm{\Lambda }𝑑k2\nu k^2E(k).$$
(18)
Using Eqs. (16) and (18), we find $`\epsilon =(27/8)\nu ^3C_\mathrm{K}^3\mathrm{\Lambda }^4`$. Finally, to complete the matching we assume that the wavevector cutoff is one-half the Kolmogorov dissipation number, $`2\mathrm{\Lambda }=k_d(\epsilon /\nu ^3)^{1/4}`$ . Putting these relations together in three dimensions, we find

$`\epsilon `$ $`=`$ $`0.0380D`$
$`C_\mathrm{K}`$ $`=`$ $`1.68`$
$`N`$ $`=`$ $`0.397.`$ (19)
Note that field theory cannot calculate non-universal parameters such as these with precision, as these results depend on the assumption in Eq. (18) and the relation between $`\mathrm{\Lambda }`$ and $`k_d`$. A more detailed matching calculation of these values using the model of turbulence in Eqs. (3) and (15) to refine the matching calculation would be of interest.

## 5 Intermittency

We here address the issue of intermittency in our model. That is, we seek to determine the scaling of the single-time structure function

$$S_{2n}(r)=\left[|𝐛^{}(𝐱+𝐫)𝐛^{}(𝐱)|^2\right]^n(\mathrm{const})r^{\zeta _{2n}}.$$
(20)
From simple dimensional analysis, we find $`\zeta _{2n}=nz`$. From explicit calculation for our model, we find an exponent that differs from this value. This difference is referred to in the fluid mechanics literature as intermittency. That $`\zeta _{2n}nz`$ is an expression of the non-Gaussian nature of the fixed point identified in our model and of the divergence of the single-time structure functions in the limit of an infinite integral length scale.

We use the same arguments about scaling of space and time as in Eq. 13 to express the original correlation function in terms of the renormalized correlation function:

$$S_{2n}(r)e^{n(22z)l}S_{2n}(r(l);l).$$
(21)
Here $`r=\mathrm{exp}(l)r(l)`$. This equation is applied until $`r(l)`$ is of the order of the dissipation length scale $`r(l^{})=2\pi /\mathrm{\Lambda }h`$. At this length scale, we then match the correlation function to a perturbation theory result:

$`S_{2n}(r(l^{});l^{})`$ $``$ $`\left[{\displaystyle \frac{h^2}{2\nu }}c(h;l^{})\right]^n`$ (22)
$`=`$ $`\left({\displaystyle \frac{h^2}{2\nu }}\right)^ne^{(3z2)nl^{}}c^n(h)`$
$`=`$ $`\left({\displaystyle \frac{h^2}{2\nu }}\right)^ne^{(3z2)nl^{}},`$
where we have used the scaling of the small-$`r`$ behavior of the $`c(r)`$ function in Eq. 15 and have used $`lim_{h0}c(h)=1`$. Combining Eqs. 21 and 22 we find

$`S_{2n}(r)(\mathrm{const})r^{nz}f_n(r/L).`$ (23)
We have here introduced the fact that traditional renormalization group arguments can determine asymptotic behavior only up to a scaling function of the ratio $`r/L`$, where $`L`$ is the macroscopic size of the system. This is because $`r`$ and $`L`$ are scaled by the same factor in the renormalization group analysis, and so no dependence on the ratio $`r/L`$ is detectable. For turbulence, $`h`$ is the dissipation length scale, and $`L`$ is the integral length scale. In many applications of renormalization group theory to condensed matter systems, the scaling function $`f_n(r/L)(\mathrm{const})`$ as $`L\mathrm{}`$, and so it does not play a role. In our case, on the other hand, the scaling function gives the corrections to intermittency.

To determine the function $`f_n`$ we use the operator product expansion . A similar strategy has proven successful in the study of turbulent transport of passive scalars . The operator product expansion states that

$$F(𝐱_1)F(𝐱_2)\underset{\alpha }{}c_\alpha (|𝐱_1𝐱_2|)F_\alpha ^R(𝐱).$$
(24)
Here the $`𝐱=(𝐱_1+𝐱_2)/2`$, and $`F_\alpha ^R`$ is the set of all renormalized operators that are generated by the renormalization group flow of $`F`$. Equation 24 is nothing more than a Taylor series expansion, where both the bare terms in the Taylor series and those terms that are generated by the renormalization flow are included. In our particular case, instead of a pair of operators, we have $`S_{2n}`$, which is a product of $`2n`$ factors on the left hand side of this equation. The important point about this expansion is that the functions $`c_\alpha (r)`$ are finite and exhibit no dependence on the system size $`L`$. Any possible system size dependence of this expansion, therefore, is contained within $`F_\alpha ^R(𝐱)`$. By comparison to Eq. 23, we see that the scaling function $`f_n`$ is thus determined by the behavior of these renormalized operators. We first determine the scaling of the renormalized operators:

$$F_\alpha ^R=e^{\mathrm{\Delta }_\alpha l^{}}F_\alpha ^R(l^{}).$$
(25)
We follow the renormalization flows until $`L(l^{})=\mathrm{exp}(l^{})L=r`$, a criterion that automatically ensures the functional form specified in Eq. 23. We thus conclude that

$$F_\alpha ^R=\left(\frac{r}{L}\right)^{\mathrm{\Delta }_\alpha }F_\alpha ^R(l^{}).$$
(26)
In fact, we determine the scaling of these operator averages by using the generating functional $`\delta S=A_\alpha 𝑑td^d𝐱F_\alpha ^R`$. We will find

$$A_\alpha (l)e^{(z+d2nz+\gamma _\alpha )l}A_\alpha (0).$$
(27)
In terms of $`A_\alpha `$ and the partition function $`Z`$, the operator average is given by

$$𝑑td^d𝐱F_\alpha ^R=\frac{d\mathrm{ln}Z}{dA_\alpha (0)}|_{A_\alpha (0)=0}.$$
(28)
This equation makes clear why Eq. 25 has the form that it does and identifies $`\mathrm{\Delta }_\alpha =2nz+\gamma _\alpha `$. The value of $`\gamma _\alpha `$ will be determined by a nontrivial fixed point of the renormalization group flow equations for $`A_\alpha `$. Combining Eqs. 26 and 27, we find

$$F_\alpha ^R\left(\frac{r}{L}\right)^{2nz\gamma _\alpha }F_\alpha ^R(l^{}).$$
(29)
The function $`F_\alpha ^R(l^{})`$ is determined by matching exactly as in Eq. 22. We, thus, find that

$`F_\alpha ^R`$ $``$ $`\left({\displaystyle \frac{r}{L}}\right)^{2nz\gamma _\alpha }\left({\displaystyle \frac{r}{L}}\right)^{(23z)n}`$ (30)
$`=`$ $`\left({\displaystyle \frac{r}{L}}\right)^{n(2z)\gamma _\alpha }.`$
To make use of this result of the operator product expansion, it remains only to calculate the value of $`\gamma _\alpha `$. Once we have this value, we find the scaling function to be

$$f_{2n}(r/L)\left(\frac{r}{L}\right)^{n(2z)\gamma _\alpha }\mathrm{as}L\mathrm{to}\mathrm{}.$$
(31)
In fact, the operator $`S_{2n}`$ will generate several new operators $`F_\alpha `$ in the expansion of Eq. 24. The appropriate value of $`\gamma _\alpha `$ to use in Eq. 31 is the largest one. These generated operators may mix upon the renormalization, in which case the appropriate value of $`\gamma _\alpha `$ is the largest eigenvalue of the flow equation matrix.

We first determine the scaling function $`f_1`$. We limit consideration to the case $`d>2`$, where $`D_\alpha ^{}=0`$. The correlation function $`S_{2n}`$ is given in a Taylor series as $`S_{2n}(r)[r^2_xb_i_xb_i]^n`$. This will generate the symmetrized operator $`F=[_jb_i_jb_i+_ib_j_jb_i]^n`$, which we consider.

For the case $`n=1`$, we consider the generating functional

$`\delta S_{II}`$ $`=`$ $`{\displaystyle }dt{\displaystyle _{𝐤_1𝐤_2}}[A^{(1)}𝐤_1𝐤_2\widehat{b}^{}(𝐤_1,t)\widehat{b}^{}(𝐤_2,t)`$ (32)
$`+A^{(0)}𝐤_1\widehat{b}^{}(𝐤_2,t)𝐤_2\widehat{b}^{}(𝐤_1,t)].`$
We find

$$\frac{dA^{(i)}}{dl}=(z+d2z)A^{(i)}+\frac{2(d1)g}{d}A^{(1)},$$
(33)
where $`g=D_\beta K_d/(\nu ^2\mathrm{\Lambda }^ϵ)`$. We, thus, identify $`\gamma _1=2(d1)g^{}/d`$. Using Eq. 10 we find in three dimensions $`z=2ϵ/4+O(ϵ^2)`$ and $`g^{}=15ϵ/4+O(ϵ^2)`$ and conclude from Eqs. 23 and 31 that

$`S_2(r)`$ $``$ $`(\mathrm{const})r^z\left({\displaystyle \frac{r}{L}}\right)^{2z\gamma _1}`$ (34)
$`=`$ $`(\mathrm{const})r^{2ϵ/4}\left({\displaystyle \frac{r}{L}}\right)^{ϵ/45ϵ}.`$
We now determine the scaling function $`f_2`$. We start with the symmetrized generating functional

$`\delta S^{(2)}`$ $`=`$ $`{\displaystyle 𝑑t_{𝐤_1𝐤_2𝐤_3𝐤_4}\widehat{b}^{}(𝐤_1,t)\widehat{b}^{}(𝐤_2,t)}`$ (35)
$`\times \widehat{b}^{}(𝐤_3,t)\widehat{b}^{}(𝐤_4,t)[A_1^{(2)}𝐤_1𝐤_2𝐤_3𝐤_4`$
$`+A_2^{(2)}(𝐤_1𝐤_3𝐤_2𝐤_4+𝐤_1𝐤_4𝐤_2𝐤_3)].`$
This term generates two additional generating functionals:

$`\delta S^{(1)}`$ $`=`$ $`{\displaystyle 𝑑t_{𝐤_1𝐤_2𝐤_3𝐤_4}\widehat{b}^{}(𝐤_1,t)\widehat{b}^{}(𝐤_2,t)}`$ (36)
$`\times [A_1^{(1)}𝐤_1𝐤_2𝐤_4\widehat{b}^{}(𝐤_3,t)𝐤_3\widehat{b}^{}(𝐤_4,t)`$
$`+A_2^{(1)}(𝐤_1𝐤_3𝐤_4\widehat{b}^{}(𝐤_3,t)𝐤_2\widehat{b}^{}(𝐤_4,t)`$
$`+𝐤_1𝐤_4𝐤_2\widehat{b}^{}(𝐤_3,t)𝐤_3\widehat{b}^{}(𝐤_4,t)`$
$`+𝐤_2𝐤_3𝐤_4\widehat{b}^{}(𝐤_3,t)𝐤_1\widehat{b}^{}(𝐤_4,t)`$
$`+𝐤_2𝐤_4𝐤_1\widehat{b}^{}(𝐤_3,t)𝐤_3\widehat{b}^{}(𝐤_4,t))`$
$`+A_3^{(1)}(𝐤_3𝐤_4𝐤_1\widehat{b}^{}(𝐤_3,t)𝐤_2\widehat{b}^{}(𝐤_4,t)`$
$`+𝐤_3𝐤_4𝐤_2\widehat{b}^{}(𝐤_3,t)𝐤_1\widehat{b}^{}(𝐤_4,t))]`$
and

$`\delta S^{(0)}`$ $`=`$ $`{\displaystyle 𝑑t_{𝐤_1𝐤_2𝐤_3𝐤_4}\widehat{b}_i^{}(𝐤_1,t)\widehat{b}_j^{}(𝐤_2,t)\widehat{b}_k^{}(𝐤_3,t)\widehat{b}_l^{}(𝐤_4,t)}`$ (37)
$`\times [A_1^{(0)}𝐤_{1}^{}{}_{j}{}^{}𝐤_{2}^{}{}_{i}{}^{}𝐤_{3}^{}{}_{l}{}^{}𝐤_{4}^{}{}_{k}{}^{}`$
$`+A_2^{(0)}(𝐤_{1}^{}{}_{j}{}^{}𝐤_{2}^{}{}_{l}{}^{}𝐤_{3}^{}{}_{i}{}^{}𝐤_{4}^{}{}_{k}{}^{}+𝐤_{1}^{}{}_{j}{}^{}𝐤_{2}^{}{}_{k}{}^{}𝐤_{3}^{}{}_{l}{}^{}𝐤_{4}^{}{}_{i}{}^{}`$
$`+𝐤_{1}^{}{}_{l}{}^{}𝐤_{2}^{}{}_{i}{}^{}𝐤_{3}^{}{}_{j}{}^{}𝐤_{4}^{}{}_{k}{}^{}+𝐤_{1}^{}{}_{k}{}^{}𝐤_{2}^{}{}_{i}{}^{}𝐤_{3}^{}{}_{l}{}^{}𝐤_{4}^{}{}_{j}{}^{})`$
$`+A_3^{(0)}(𝐤_{1}^{}{}_{k}{}^{}𝐤_{2}^{}{}_{l}{}^{}𝐤_{3}^{}{}_{i}{}^{}𝐤_{4}^{}{}_{j}{}^{}+𝐤_{1}^{}{}_{k}{}^{}𝐤_{2}^{}{}_{l}{}^{}𝐤_{3}^{}{}_{j}{}^{}𝐤_{4}^{}{}_{i}{}^{}`$
$`+𝐤_{1}^{}{}_{l}{}^{}𝐤_{2}^{}{}_{k}{}^{}𝐤_{3}^{}{}_{i}{}^{}𝐤_{4}^{}{}_{j}{}^{}+𝐤_{1}^{}{}_{l}{}^{}𝐤_{2}^{}{}_{k}{}^{}𝐤_{3}^{}{}_{j}{}^{}𝐤_{4}^{}{}_{i}{}^{})]`$
We will find $`A_2^{(1)}(l)=A_3^{(1)}(l)`$. Although we have included $`A_3^{(0)}`$ for generality, we will find that this term is not generated, and $`A_3^{(0)}(l)=0`$.

A lengthy calculation shows that

$$\frac{d𝐀}{dl}=(z+d4z)I𝐀+\frac{g}{d(d+2)(d+4)}M𝐀$$
(38)
where the vector $`𝐀=(A_1^{(2)},A_2^{(2)},A_1^{(1)},A_2^{(1)},A_3^{(1)},A_1^{(0)},A_2^{(0)})`$, and the matrix $`M`$ is given in Table 1.

We diagonalize the matrix $`M`$, finding the eigenvalues $`(0,0,0,0,d^3+4d^24d16,2(d^3+5d^2+2d8),4(d^3+7d^2+6d))`$. The largest eigenvalue is the last one for all $`d`$. In three dimensions we identify $`\gamma _2=432g^{}/105`$. Using Eqs. 23 and 31 we find that

$$S_4(r)(\mathrm{const})r^{4ϵ/2}\left(\frac{r}{L}\right)^{ϵ/2108ϵ/7}.$$
(40)
We have, therefore, derived the intermittency corrections to dimensional analysis for the present model. The corrections are calculated in a controlled fashion and are proportional to $`ϵ`$. The coefficients of the correction are not small, and for finite $`ϵ`$, higher order terms in the expansion are required for an accurate estimation of the effects of intermittency.

One might wonder whether there are any corrections to the Kolmogorov energy cascade, Eq. 14, that arise from the operator product expansion. More generally, are there any corrections to $`[b_i^{}(𝐱,t)b_i^{}(𝐱,t)]^n`$? There are no such corrections. This type of operator flows under the renormalization group only to operators with more derivatives, such as $`_lb_j^{}_mb_k^{}(b_i^{}b_i^{})^{n1}`$. These operators are less relevant than the original operator, and so they make no contribution to the scaling at leading order. In the language of Eq. 24, $`F_\alpha ^R1`$, and the scaling function $`f_n(r/L)(\mathrm{const})`$ as $`L\mathrm{}`$.

## 6 Conclusion

An alternative, simpler model would have been to take the turbulent forces to be proportional to the velocity, rather than the velocity gradient. The simplest model, moreover, would take the forces to be white noise in time and uncorrelated in each of the spatial dimensions. For this model to be nontrivial, a mean fluid flow must be introduced . Interestingly, when this is done for forces that are random in time as well as space, the resulting theory has the same flow equations as the traditional random force model of turbulence . Also of interest to note is that a theory with velocity-gradient-dependent random forces that are white noise in time would have no renormalization of any parameter, as the diagrams of Fig. 2 would all vanish due to the causality of the bare propagator, Eq. (LABEL:5aaa).

The scaling of the Kolmogorov energy cascade is determined once the value of the dynamical exponent is fixed, *i.e.* $`z=2/3`$ for isotropic turbulence. In random force models such as the present one, the scaling of the energy cascade simply serves to fix the correlation function of the random forcing. The predictive power of models such as these lie in their ability to provide nontrivial predictions of the intermittency corrections. In the present model, we are able to provide these corrections as a systematic expansion in $`ϵ`$.

In summary, we have introduced a new statistical mechanics model for isotropic turbulence. This model makes use of a random, velocity-gradient-dependent force. This model is both consistent with practical, engineering-type calculations and well-defined within the renormalization group $`ϵ`$ expansion. This model makes stronger the analogy between turbulence and critical phenomena.

Owing to the irrelevance of the convection terms at the fixed point, our results may alternatively be viewed as an analysis of transport in a new class of random media.

###### Acknowledgements.

This research was supported by the National Science Foundation and by an Alfred P. Sloan Foundation Fellowship to M.W.D.

## Appendix A: One-Loop contributions to $`\nu `$

We here show how the diagrams of Fig. 2a contribute to the propagator, Fig. 1a. Each of the terms is associated with one of the $`D_{ijk}^{lmn}`$ terms in Eq. (LABEL:5). In the calculation of the averages on the shell, there are five terms associated with $`D_\alpha `$ and three associated with $`D_\beta `$, as the last two terms associated with $`D_\alpha `$ are identical, and the last two terms associated with $`D_\beta `$ are also identical. The first term, associated with $`D_\alpha `$, is

$`I_1`$ $`=`$ $`2\times {\displaystyle \frac{D_\alpha }{2}}{\displaystyle 𝑑t_1𝑑t_2_{𝐤_1𝐤_2𝐤_3𝐤_4}}`$ (41)
$`\times (2\pi )^d\delta (𝐤_1+𝐤_2+𝐤_3+𝐤_4)|𝐤_1+𝐤_2|^y`$
$`\times \mathrm{\Theta }(t_1t_2)(2\pi )^d\delta (𝐤_2+𝐤_3)e^{\nu k_2^2(t_1t_2)}`$
$`\times \left\{𝐤_2𝐤_4\widehat{\overline{b}}_i^{}(𝐤_1,t_1)\widehat{b}_j^{}(𝐤_4,t_2)\widehat{\mathrm{\Pi }}_{ij}(𝐤_2)\right\}`$
$`=`$ $`{\displaystyle \frac{K_d}{\nu ^2}}\left[D_\alpha {\displaystyle \frac{y(d+1)}{d(d+2)}}\right]{\displaystyle _{\mathrm{\Lambda }/a}^\mathrm{\Lambda }}𝑑qq^{dy3}`$
$`\times {\displaystyle }dt{\displaystyle _𝐤}\nu k^2\widehat{\overline{b}}_i^{}(𝐤,t)\widehat{b}_i^{}(𝐤,t)`$
For the remaining contributions, we list the symmetry factor, terms in braces in the integrand of Eq. (41) that change, and the final contribution in brackets that change. The contributions are shown in Table 2, where the dependence of the fields upon time has been suppressed. The term $`e^{\nu k_i^2(t_1t_2)}`$, which has the same momentum argument as the $`\widehat{\mathrm{\Pi }}(𝐤_i)`$ term, has been suppressed. The delta function is also suppressed. Summing all these contributions to $`\nu `$, we get the first flow equation of Eq. (9).

## Appendix B: One-Loop contributions to $`\lambda `$

We here show how the diagrams of Fig. 2b contribute to the convection term, Fig. 1b. Each of the terms is associated with one of the $`D_{ijk}^{lmn}`$ terms in Eq. (LABEL:5). It is convenient to define the convection operator $`\widehat{M}_{ijk}(𝐤)=k_j\widehat{\mathrm{\Pi }}_{ik}(𝐤)+k_k\widehat{\mathrm{\Pi }}_{ij}(𝐤)`$. The convection term of Eq. (LABEL:5) then becomes

$`+`$ $`{\displaystyle \frac{i\lambda }{2}}{\displaystyle _{𝐤_1𝐤_2𝐤_3}}{\displaystyle 𝑑t(2\pi )^d\delta (𝐤_1+𝐤_2+𝐤_3)}`$ (42)
$`\times \widehat{M}_{ijk}(k_1)\widehat{\overline{b}}_i^{}(𝐤_1,t)\widehat{b}_j^{}(𝐤_2,t)\widehat{b}_k^{}(𝐤_3,t).`$
In the calculation of the averages on the shell, there are again five terms associated with $`D_\alpha `$ and three associated with $`D_\beta `$. The first such term is

$`J_1`$ $`=`$ $`2\times \left({\displaystyle \frac{1}{2!}}\right)2\left({\displaystyle \frac{i\lambda }{2}}\right){\displaystyle \frac{D_\alpha }{2}}{\displaystyle 𝑑t_1𝑑t_2}`$ (43)
$`\times {\displaystyle _{𝐤_1𝐤_2𝐤_3𝐤_4}}(2\pi )^d\delta (𝐤_1+𝐤_2+𝐤_3+𝐤_4)`$
$`\times {\displaystyle }dt_3{\displaystyle _{𝐤_5𝐤_6𝐤_7}}(2\pi )^d\delta (𝐤_5+𝐤_6+𝐤_7)`$
$`\times |𝐤_1+𝐤_2|^y\mathrm{\Theta }(t_1t_3)\mathrm{\Theta }(t_3t_2)`$
$`\times (2\pi )^d\delta (𝐤_2+𝐤_5)(2\pi )^d\delta (𝐤_6+𝐤_3)`$
$`\times e^{\nu k_2^2(t_1t_3)}e^{\nu k_6^2(t_3t_2)}\widehat{M}_{lmn}(𝐤_5)`$
$`\times \{k_{2_k}k_{4_k}\widehat{\overline{b}}_i^{}(𝐤_1,t_1)\widehat{b}_j^{}(𝐤_4,t_2)\widehat{b}_n^{}(𝐤_7,t_3)`$
$`\times \widehat{\mathrm{\Pi }}_{jl}(𝐤_2)\widehat{\mathrm{\Pi }}_{im}(𝐤_6)\}`$
$`=`$ $`{\displaystyle \frac{i\lambda }{2}}{\displaystyle \frac{K_d}{\nu ^2}}\left[D_\alpha {\displaystyle \frac{(d+1)}{d(d+2)}}\right]{\displaystyle _{\mathrm{\Lambda }/a}^\mathrm{\Lambda }}𝑑qq^{dy3}`$
$`\times {\displaystyle }dt{\displaystyle _{𝐤_1𝐤_2𝐤_3}}(2\pi )^d\delta (𝐤_1+𝐤_2+𝐤_3)`$
$`\times \widehat{M}_{ijk}(𝐤_1)\widehat{\overline{b}}_i^{}(𝐤_1,t)\widehat{b}_j^{}(𝐤_2,t)\widehat{b}_k^{}(𝐤_3,t)`$
For the remaining contributions, we list the symmetry factor, terms in braces in the integrand of Eq. (43) that change, and the final contribution in brackets that change. The contributions are shown in Table 3, where again the dependence on time has been suppressed. The terms $`e^{\nu k_i^2(t_1t_3)}`$ and $`e^{\nu k_j^2(t_3t_2)}`$, which have the same momentum arguments as the two $`\widehat{\mathrm{\Pi }}(𝐤_i)`$ terms, have been suppressed. The delta functions have also been suppressed. Summing all these contributions to $`\lambda `$, we get the second flow equation of Eq. (9).

## Appendix C: One-Loop contributions to $`D_\alpha ,D_\beta `$

We here show how the diagram of Fig. 2c contribute to the disorder term, Fig. 1c. Due to the symmetry of the $`D_{ijk}^{lmn}`$ term in Eq. (4), the four possible types of diagrams in Fig. 2c contribute the same value. The result is

$`L`$ $`=`$ $`2^3\left({\displaystyle \frac{1}{2!}}\right)\left({\displaystyle \frac{1}{2}}\right)^2`$ (44)
$`\times {\displaystyle }dt_1dt_2{\displaystyle _{𝐤_1𝐤_2𝐤_3𝐤_4}}(2\pi )^d\delta (𝐤_1+𝐤_2+𝐤_3+𝐤_4)`$
$`\times {\displaystyle }dt_3dt_4{\displaystyle _{𝐤_5𝐤_6𝐤_7𝐤_8}}(2\pi )^d\delta (𝐤_5+𝐤_6+𝐤_7+𝐤_8)`$
$`\times |𝐤_1+𝐤_2|^yk_{2_j}k_{4_m}\widehat{\overline{b}}_i^{}(𝐤_1,t_1)\widehat{b}_k^{}(𝐤_2,t_1)`$
$`\times |𝐤_5+𝐤_6|^yk_{6_s}k_{8_v}\widehat{\overline{b}}_r^{}(𝐤_5,t_3)\widehat{b}_w^{}(𝐤_8,t_4)`$
$`\times e^{\nu k_6^2(t_3t_2)}\mathrm{\Theta }(t_3t_2)\widehat{\mathrm{\Pi }}_{lt}(𝐤_6)(2\pi )^d\delta (𝐤_6+𝐤_3)`$
$`\times e^{\nu k_4^2(t_2t_4)}\mathrm{\Theta }(t_2t_4)\widehat{\mathrm{\Pi }}_{nu}(𝐤_4)(2\pi )^d\delta (𝐤_4+𝐤_7)`$
$`\times D_{ijk}^{lmn}D_{rst}^{uvw}`$
$`=`$ $`{\displaystyle \frac{K_d}{2\nu ^2}}[{\displaystyle \frac{2}{d}}\delta _{ms}\delta _{lt}\delta _{nu}`$
$`{\displaystyle \frac{2}{d(d+2)}}\delta _{nu}\left(\delta _{ms}\delta _{lt}+\delta _{ml}\delta _{st}+\delta _{mt}\delta _{ls}\right)`$
$`{\displaystyle \frac{2}{d(d+2)}}\delta _{lt}\left(\delta _{ms}\delta _{nu}+\delta _{mn}\delta _{su}+\delta _{mu}\delta _{sn}\right)`$
$`+{\displaystyle \frac{2}{d(d+2)(d+4)}}M_{lmn}^{stu}]D_{ijk}^{lmn}D_{rst}^{uvw}`$
$`\times {\displaystyle _{\mathrm{\Lambda }/a}^\mathrm{\Lambda }}dqq^{dy3}`$
$`\times {\displaystyle }dt_1dt_2{\displaystyle _{𝐤_1𝐤_2𝐤_3𝐤_4}}(2\pi )^d\delta (𝐤_1+𝐤_2+𝐤_3+𝐤_4)`$
$`\times |𝐤_1+𝐤_2|^y𝐤_{2}^{}{}_{j}{}^{}𝐤_{4}^{}{}_{v}{}^{}`$
$`\times \widehat{\overline{b}}_i^{}(𝐤_1,t_1)\widehat{b}_k^{}(𝐤_2,t_1)\widehat{\overline{b}}_r^{}(𝐤_3,t_2)\widehat{b}_w^{}(𝐤_4,t_2)`$
It is clear that to evaluate this expression, we need to evaluate a term such as

$`{\displaystyle \underset{lmn,stu}{}}`$ $`k_{2_j}k_{4_v}\widehat{\overline{b}}_i^{}(𝐤_1,t_1)\widehat{b}_k^{}(𝐤_2,t_1)\widehat{\overline{b}}_r^{}(𝐤_3,t_2)\widehat{b}_w^{}(𝐤_4,t_2)`$ (45)
$`\times D_{ijk}^{lmn}D_{rst}^{uvw}\left[\delta _{lt}\delta _{ms}\delta _{nu}\right]`$
Fourteen other terms need to be evaluated in order to calculate the total contribution from Eq. (44). These terms each contribute in a form that can be cast as a contribution to $`D_\alpha `$ and $`D_\beta `$ in Eq. (LABEL:5). Shown in Table 4 are the terms and their contributions.

Summing all these contributions to $`D_\alpha `$ and $`D_\beta `$, we get the third and fourth flow equations of Eq. (9).


format: text

text:
# Aging and species abundance in the Tangled Nature model of biological evolution.

## I Introduction

The dynamics and organization of biological ecosystems is a fascinating example of complex interacting systems with many levels of emerging structure and time scales. Biological evolution creates intricate taxonomic hierarchies presumably as an effect of mutation, natural selection and the ensuing adaptation. Taxonomic structures from the level of individuals through species and genera up to kingdoms are generated and vanish again in a never ending succession. Different strata in the hierarchy are described by very different timescales and with very different types of dynamics. At the level of individuals, fairly well defined characteristic lifetimes exist for each specific type (species) and the population dynamics can be considered smooth. This picture changes as one considers the system at the more coarse grained level of species and genera. The lifetime distribution of, e.g., genera is broad (see e.g. Newman and Sibani (1999)) and the dynamics is intermittent Eldredge (1971); Gould (1977); Gould and Eldredge (1977); Eldredge and Gould (1988). In the spirit of the traditional approach of statistical mechanics it is interesting to consider models, defined at a microscopic level, which are able to reproduce the large scale temporal and taxonomic structures.

In the present paper we consider a model of individuals identified solely by their genome. The model was introduced in Ref. Christensen et al. (2002) where we also presented a discussion of the qualitative behavior of the model. We combine ecology with evolution by considering interacting individuals which can multiply (sexually or asexually) subject, potentially, to mutations. The size of the total population fluctuates, the average being controlled by the amount of available resources. From these three minimal ingredients emerge segregation in genome space, to be interpreted as the appearance of species, and a complex intermittent dynamics, to be interpreted as extinction and creation events at the higher taxonomic levels. The entire taxonomic hierarchy is an emergent property of the dynamics at the microscopic level of individuals. We characterize the configurations generated in genotype space in terms of the species abundance curve, and find a good qualitative agreement with the functional form typically found for real ecosystems. The intermittent dynamics is characterized by the statistics of the duration of the quasi-stable epochs or in other words the waiting times between transitions. We find a broad distribution of durations and observe a gradual aging of the macro-dynamics. No stationary state is ever reached.

### I.1 Related Models

Many mathematical models of biological evolution have been developed according to the usual statistical mechanics agenda of generating the macroscopic complex behavior from simplistic microscopic definitions. An elegant review of this endeavor has recently been given by Drossel Drossel (2001). Here we limit ourselves to a discussion of similarities and differences between our model and related studies.

Let us first mention models which define the ecosystem in terms of individuals. Higgs and Derrida Higgs and Derrida (1992) studied speciation in a model consisting of a fixed number of individuals. Each individual is represented by a genome modeled as a string of zeros or ones, like in Eigen and collaborators’ seminal work on quasi-species Eigen et al. (1988). Higgs and Derrida demonstrated that a sexually reproducing population breaks up into distinct species when only individuals with a sufficiently similar genome sequence are allowed to produce offspring. This agrees with a large bulk of experimental work Hoster (1993). Gavrilets and collaborators Gravilets (1999); Gavrilets et al. (1998, 2000) have made use of similar models generalized in particular to be able to study geographical and temporal aspects of speciation. These studies differ from ours in assuming a fixed population size and by defining a fitness function for pairs of individuals which is constant if the Hamming distance between the genomes of two individuals is small enough, and zero otherwise. Our model allows the total size of the population to fluctuate and the fitness of pairs of individuals (or in the asexual case single individuals) depends on the composition of the population at a given instant in time.

It is also important to mention the fitness landscape approach first pioneered by Wright Wright (1932, 1988), who considered gene frequencies, and was brought to the attention of the statistical mechanics community mainly through Kauffman’s so called NK model Kauffman and Levine (1987); Kauffman (1993). The main focus of the NK model, and of the later co-evolutionary NKC model Kauffman (1993), is the study of epistatic interactions (the influence of one gene on another) by use of fitness functions. The main difference between our model and Kauffman’s models is that the fitness of an individual in our system depends on the frequencies by which other locations in genotype space are occupied.

Taylor and Higgs Taylor and Higgs (2000) have studied pleiotropy and epistasis (the influence of one gene on several traits and the influence of one gene on another) in a model that combines and generalizes aspects of the Higgs-Derrida model with the epistatic interactions of Kauffman’s models. Taylor and Higgs then derive a phenotypical fitness for the specific genotype. Kaneko and Yomo Kaneko and Yomo (2000) have also studied models in which the difference between phenotype and genotype is accounted for explicitly. In our model we make the drastic simplification not to distinguish between genotype and phenotype.

Other models consider species as the elementary building block; these models neglect the specifics of the dynamics arising from reproduction and mutations at the level of individuals. The simplest of these models is the Bak-Sneppen model Bak and Sneppen (1993). The model aims to demonstrate that co-evolutionary interactions are sufficient to produce intermittent dynamics which is then related to intermittency in the fossil record and to Eldredge and Gould’s concept of punctuated equilibrium Eldredge (1971); Gould (1977); Gould and Eldredge (1977); Eldredge and Gould (1988). Each species is characterized by a single number between zero and one, the fitness, and the total number of species is kept constant. The model has interesting statistical properties but is difficult to relate to biological evolution.

Species level models of more detail than the Bak-Sneppen model have been formulated recently by McKane, Alonso and Solé McKane et al. (2000) and by Drossel, Higgs and McKane Drossel et al. (2001). The emphasis in these models is on predator-prey interactions and food-webs and are generalizations of early work by May May (1972) and May and Anderson May and Anderson (1983). Our model is intended to include all types of interactions between individuals, e.g. antagonistic or collaborative relationships, in addition to predator-prey competitions. Another important difference is that we define our model at the level of individuals in order to be able to study the emergence of species, something not possible in a species based model.

Most models of biological evolution assume that the dynamics is in a statistically stationary state. One marked exception is the model considered by Sibani and collaborators Sibani et al. (1995); Sibani (1997); Sibani and Brandt (1998); Sibani and Pedersen (1999). This is an abstract species based model consisting of random walks in a rugged fitness landscape. The statistics of the jumps in this landscape are the same as the record statistics considered some time ago by Sibani and Littlewood Sibani and Littlewood (1993). The pace of the dynamics of the model gradually slows down as indicated by a logarithmically decreasing extinction rate. As we shall see below our individual based model also exhibits aging, a property found to be consistent with analysis of the fossil record Newman and Sibani (1999).

The paper is organized as follows. In the next section we define the model in detail. In Sec. III we discuss the modes of the models emergent dynamics. In Sec. IV we show how the configurations generated dynamically gradually become better adapted in a collective sense. Sec. V demonstrates that the ecologies generated in the model exhibit characteristics similar to those observed in real ecologies. In Sec. VI we discuss the models ability to address the issue concerning why sexual reproduction can compete evolutionary with asexual reproduction. Sec. VII contains an analysis of the error threshold. We briefly present in Sec. VIII a scan of the behavior of the model for a range of the control parameters and in Sec. IX we conclude and summarize.

## II Definition of Model

We describe here in detail the structure and dynamics of the model which we, with an allusion to Darwin’s notion of the Tangled Bank, called the Tangled Nature (TaNa) model to stress the model’s emphasis on ecological interactions Christensen et al. (2002).

### II.1 Interaction

We represent an individual by a vector $`𝐒^\alpha =(S_1^\alpha ,S_2^\alpha ,\mathrm{},S_L^\alpha )`$ in genotype space $`𝒮`$. This representation is frequently used, see e.g. Ref. Eigen et al. (1988); Kauffman (1993); Higgs and Derrida (1992); Gravilets (1999); Wagner et al. (1998). Here $`S_i^\alpha `$ may take the values $`\pm 1`$, i.e. $`𝐒^\alpha `$ denotes one of the corners of the $`L`$ dimensional hypercube (in the present paper we use $`L=20`$). The coordinates $`S_i^\alpha `$ may be interpreted as genes with two alleles, or a string of either pyrimidines or purines. We think of genotype space $`𝒮`$ as containing all possible ways of combining the genomic building blocks into genome sequences. Many sequences may not correspond to viable organisms. Whether this is the case or not is for the evolutionary dynamics to determine. All possible sequences are made available for evolution to select from.

Individuals are labeled by Greek letters $`\alpha ,\beta ,\mathrm{}=1,2,\mathrm{},N(t)`$. When we refer, without reference to a specific individual, to one of the $`2^L`$ positions in genome space, we use roman superscripts $`𝐒^a`$, $`𝐒^b`$, … with $`a,b,\mathrm{}=1,2,\mathrm{},2^L`$. Many different individuals $`𝐒^\alpha ,𝐒^\beta `$,…, may reside on the same position, say $`𝐒^a`$, in $`𝒮`$.

The ability of an individual $`\alpha `$ to reproduce is controlled by $`H(𝐒^\alpha ,t)`$:

$$H(𝐒^\alpha ,t)=\frac{1}{cN(t)}\underset{𝐒𝒮}{}J(𝐒^\alpha ,𝐒)n(𝐒,t)\mu N(t),$$
(1)
where $`c`$ is a control parameter (see below), $`N(t)`$ is the total number of individuals at time $`t`$, the sum is over the $`2^L`$ locations $`𝐒`$ in $`𝒮`$ and $`n(𝐒,t)`$ is the occupancy of position $`𝐒`$. Two positions $`𝐒^a`$ and $`𝐒^b`$ in genome space are coupled with the fixed random strength $`J^{ab}=J(𝐒^a,𝐒^b)`$ which can be either positive or negative or zero. The coupling is non-zero with probability $`\mathrm{\Theta }`$ (throughout the paper we use $`\mathrm{\Theta }=0.25`$), in which case we assume $`J^{ab}J^{ba}`$ to be a deterministic but erratic function of the two positions $`𝐒^a`$ and $`𝐒^b`$. We have checked that the specific details of the form of the distribution of the non-zero values of the function $`J(𝐒^a,𝐒^b)`$ are irrelevant. We choose accordingly a form mainly determined by its numerical efficiency. In the next subsection we describe the details of the specific procedure used. The distribution of the generated interaction strengths is shown in Fig. 3 below.

#### II.1.1 Generation of interaction matrix

The interaction between two locations in genotype space, $`𝐒^a`$ and $`𝐒^b`$ is generated as a product $`J(𝐒^a,𝐒^b)=\mathrm{\Theta }(𝐒^a,𝐒^b)I(𝐒^a,𝐒^b)`$. The first factor $`\mathrm{\Theta }(𝐒^a,𝐒^b)`$ is obtained by interpreting the sequences $`𝐒^a`$ and $`𝐒^b`$ as binary numbers (letting $`10`$) and perform the XOR operation on the binary pair to obtain a new integer. This integer is used as an index in a lookup list to obtain either a 0 or 1 as the value of $`\mathrm{\Theta }(𝐒^a,𝐒^b)`$ . In case 1 is returned, the element of the $`I(𝐒^a,𝐒^b)`$ matrix is obtained in a similar way. This time, however, two arrays are needed. Each auxiliary array is of length $`2^L`$ and now the arrays contain uniformly distributed random numbers drawn from the interval $`[1,+1]`$. The pair of arrays is necessary in order to reproduce the asymmetry of the $`I(𝐒^a,𝐒^b)`$ matrix. Two indices are generated from the $`𝐒^a`$ and $`𝐒^b`$. The first via the same XOR operation used to calculate the $`\mathrm{\Theta }(𝐒^a,𝐒^b)`$ matrix element, whereas the second is simply the integer representing $`𝐒^b`$. The strength of interaction is taken to be the product of the members of each array at the appropriate location. This ensures that the elements of the matrices are non-symmetric due to the second array index depending on the order of the operation. This procedure is numerically extremely efficient and deterministic, but has the side effect of generating a distribution of a slightly unusual form, see Fig. 3.

We stress that the coupling matrix $`J(𝐒^a,𝐒^b)`$ is meant to included all possible interactions between two individuals of a given genomic constitution. In our simplistic approach, a given genome is imagined to lead uniquely to a certain set of attributes (phenotype) of the individuals/organisms. The locations $`𝐒^a`$ and $`𝐒^b`$ represent blueprints for organisms that exist in potentia. The positions may very likely be unoccupied but, if we were to construct individuals according to the sequences $`𝐒^a`$ and $`𝐒^b`$ the two individuals would have some specific features. The relationship between an organism of design $`𝐒^a`$ and one of design $`𝐒^b`$ may be as predator and prey or parasitic, i.e. $`J^{ab}>0`$ and $`J^{ba}<0`$, but it can also be collaborative ($`J^{ab}>0`$ and $`J^{ba}>0`$) or antagonistic ($`J^{ab}<0`$ and $`J^{ba}<0`$), see Fig. 1. And certainly in some cases $`J^{ab}`$ may represent less direct couplings, e.g. some animals may not eat trees, nevertheless they breath the oxygen produced by the rain forest. In order to emphasis co-evolutionary aspects we have excluded “self-interaction” among individuals located at the same position $`𝐒`$ in genome space, i.e. $`J(𝐒,𝐒)=\mathrm{𝟎}`$ for all $`𝐒𝒮`$. It is important to mention that including self-interactions of the same order of strength as the $`J`$-couplings does not change the qualitative behavior of the model.

The conditions of the physical environment are simplistically described by the term $`\mu N(t)`$ in Eq. (1), where $`\mu `$ determines the average sustainable total population size. That is, the total carrying capacity of the environment. An increase in $`\mu `$ corresponds to harsher physical conditions. This is a simplification, though one should remember that what is often considered as the physical conditions, e.g. temperature or oxygen density, is to a degree determined by the activity of other organisms and is therefore really a part of the biotic conditions. Consider, for example, the environment experienced by the bacterial flora in the intestines. Here one type of bacteria live very much in an environment strongly influenced by the presence of other types of bacteria. In this sense some fluctuations in the environment may be thought of as included in the coupling matrix $`J(𝐒^a,𝐒^b)`$.

### II.2  Reproduction, mutations and annihilation

Asexual reproduction consists of one individual being replaced by two copies. Successful reproduction occurs for individuals $`𝐒^\alpha `$ with a probability per time unit given by

$$p_{off}(𝐒^\alpha ,t)=\frac{\mathrm{exp}[H(𝐒^\alpha ,t)]}{1+\mathrm{exp}[H(𝐒^\alpha ,t)]}[0,1].$$
(2)
In the case of sexual reproduction an individual $`𝐒^\alpha `$ is picked at random and paired with another randomly chosen individual $`𝐒^\beta `$ with Hamming distance $`d=\frac{1}{2}_{i=1}^L|S_i^\alpha S_i^\beta |d_{max}`$ (allowing at most $`d_{max}`$ pairs of genes to differ). The pair produces an offspring $`\gamma `$ with probability $`\sqrt{p_{off}(𝐒^\alpha ,t)p_{off}(𝐒^\beta ,t)}`$, where $`S_i^\gamma `$ is chosen at random from one of the two parent genes, either $`S_i^\alpha `$ or $`S_i^\beta `$. For $`d_{max}1`$ this procedure may be thought of as being similar to recombination. The maximum separation criterion has been studied by several authors, see e.g. Higgs and Derrida (1992); Gravilets (1999).

We allow for mutations in the following way: with probability $`p_{mut}`$ per gene we perform a change of sign $`S_i^\gamma S_i^\gamma `$, during the reproduction process.

For simplicity, an individual is removed from the system with a constant probability $`p_{kill}`$ per time step (we use $`p_{kill}=0.2`$). This procedure is implemented both for asexual and sexually reproducing individuals.

A time step consists of one annihilation attempt followed by one reproduction attempt. One generation consists of $`N(t)/p_{kill}`$ time steps, which is the average time taken to kill all currently living individuals.

Initially we place $`N(0)=500`$ individuals at randomly chosen positions. The results are independent of initial conditions. We obtain the same results if all individuals are located at the same position initially.

The present paper’s main focus is on the asexual mode of reproduction and results presented are for asexual individuals except otherwise stated. For completeness and for comparison we do, however, in Secs. VI and VIII consider sexual reproduction, simplicity which is defined as follows.

## III  Dynamical stability

Neglecting fluctuations in the occupancy $`n(𝐒,t)`$ the above dynamics is described by the following set of equations (one equation for each position in the genotype space):

$`n(𝐒,t+1)=n(𝐒,t)`$ (3)
$`+`$ $`\{p_{off}(𝐒,t)[2(1p_{mut})^L1]p_{kill}\}{\displaystyle \frac{n(𝐒,t)}{N(t)}}`$
$`+`$ $`2p_{mut}(1p_{mut})^{L1}{\displaystyle \underset{𝐒^{},𝐒}{}}p_{off}(𝐒^{},t){\displaystyle \frac{n(𝐒^{},t)}{N(t)}},`$
where the sum is over the nearest neighbors of $`𝐒`$. Stationary solutions require the system to find configurations in genotype space for which all positions satisfy the demand that either $`n(𝐒,t)=0`$ or if $`n(𝐒,t)0`$ (neglecting the mutational back flow represented by the last term in Eq. (3)), we must have

$$p_{off}=\frac{p_{kill}}{2(1p_{mut})^L1}p_{qESS}$$
(4)
The fitness $`p_{off}(𝐒^a,t)`$ of individuals at a position $`𝐒^a`$ depends on the occupancy $`n(𝐒^b,t)`$ of all the sites $`𝐒^b`$ with which site $`𝐒^a`$ is connected through couplings $`J^{ab}`$. Accordingly, a small perturbation in the occupancy at one position may be able to disturb the balance in Eq. (4) between $`p_{off}(𝐒,t)`$, $`p_{kill}`$ and $`p_{mut}`$ on connected sites. In this way an imbalance at one site can spread as a chain reaction through the system, possibly causing a global reconfiguration of the occupancy in genotype space.

We show in Fig. 2 the occupancy in genotype space plotted as a function of time for asexual reproduction. Periods of stable configurations are separated by fast transitions. We have called the stable periods “quasi-Evolutionary Stable Strategies” or q-ESS since they are reminiscent of the Evolutionary Stable Strategies (ESS) introduced by Maynard Smith Maynard Smith (1982).

It is interesting to investigate just how stable the q-ESS are. We have done this by applying different types of perturbations in the q-ESS. The result is that the q-ESS are very stable against global perturbations such as a brief or a lasting increase in control parameters $`\mu `$, $`c`$ or $`p_{kill}`$. Changes of up to 50% in these parameters, either permanently or for a period of 100 generations, only effect the total population size and is typically not able to kick the population out of its present q-ESS configuration in genotype space. In contrast, a similar perturbation of the mutation rate easily destabilizes the q-ESS configuration.

We stress that the segregation (or speciation) to be discussed below is an effect of different couplings between different positions $`𝐒^a`$ and $`𝐒^b`$. When we assume $`J(𝐒^a,𝐒^b)=J_0`$ independent of $`𝐒^a`$ and $`𝐒^b`$, the population is not concentrated around a subset of the positions in genotype space, instead the population is smeared out through the space in a diffuse manner. Self-interaction, however, can cause segregation in a rather trivial way. Namely, if we include a distribution of $`J(𝐒,𝐒)`$ values, segregation may occur even in the case where all interaction terms assume the same value: $`J(𝐒^a,𝐒^b)=J_0`$ for $`𝐒^a𝐒^b`$. However, this type of selection of configurations in genotype space is not very interesting since the sites to become occupied is determined by the arbitrarily assigned self-interactions $`J(𝐒,𝐒)`$ and not by the collective dynamical adaptation at play when $`J(𝐒,𝐒)=0`$ and $`J(𝐒^a,𝐒^b)`$ assumes a distribution of different $`J`$-values. In reality one will expect the selection of species to be caused by a mixture of self-interaction and interaction between different species. To decide which one is dominant might be difficult and will certainly be system specific.

There is a significant difference between the distribution of active couplings, $`p_{act}(J(𝐒^a,𝐒^b))`$, in the q-ESS and the distribution during the hectic transitions. We show in Fig. 3 the distribution from which the $`J_{bare}(𝐒^a,𝐒^b)`$ are sampled together with the distribution of couplings between occupied sites after a large number of generations. During the hectic phases there is clearly no noticeable difference between the “bare” distribution of the $`J(𝐒^a,𝐒^b)`$ and the distribution of active couplings, i.e. couplings between occupied positions. During the q-ESS we observe a slight bias towards positive $`J`$-values of the active couplings. This slight shift towards more positive couplings will, according to Eqs. 1 and 2, lead to an increased reproduction rate during the q-ESS. The manifestation of this difference between the hectic periods and the q-ESS is illustrated in Fig. 4. We see that the distribution of $`H`$-values in the q-ESS during both modes of the dynamics contains a narrow peak. In the q-ESS, the peak in $`p(H)`$ is separated from a strongly negative band of support. The values of $`H`$ in this band are so negative that the corresponding $`p_{off}(H)`$ are negligible (see the insert in Fig. 4). Genotype positions corresponding to this band consist of unfit positions next to highly occupied and very fit positions. The reason these positions are occupied at all is that they are supplied by mutations occurring on the neighboring fit positions. The conclusion of these considerations is that the dynamics during the q-ESS as well as during the hectic periods are controlled by the reproduction of individuals with $`H`$-values in the two respective peaks of $`p(H)`$.

The location of the peaks of $`p(H)`$ is determined in the following way. During the hectic periods the occupation of positions in genotype space is highly unstable and $`n(𝐒,t+1)`$ is only related to $`n(𝐒,t)`$ in an erratic way, the balance equation Eq. (3) is never fulfilled for nonzero $`n(𝐒,t)=0`$. The only constraint on $`p_{off}`$ during the hectic periods is accordingly that the total population remains constant on average which implies that on average $`p_{off}=p_{kill}`$. This explains why in Fig. 4 the peak in $`p(H)`$ during the hectic periods corresponds to a peak in $`p(p_{off})`$ centered at $`p_{kill}=0.2`$. The situation is different during the q-ESS. Here the occupation of the selected positions in genotype space remains approximately constant and Eq. (4) applies. Substituting the relevant values $`p_{kill}=0.2`$, $`p_{mut}=0.0125`$ and $`L=20`$ into Eq. (4) produces $`p_{off}=0.36`$ which explains the position of the peak in $`p(p_{off})`$ during the q-ESS.

## IV  Aging

For simplicity we concentrate again on the asexual model in this section. For large genome length $`L`$ the system is always in a transient. The time needed to reach the stationary state increases exponentially with $`L`$ and is therefore unreachable for any biologically relevant values of $`L`$.

### IV.1 Increasing q-ESS durations

The gradual change in the statistical measures of the model is seen directly as a slow increase with time of the average duration of the q-ESS. To demonstrate this we show in Fig. 5 the average number of transitions $`\mathrm{\Omega }_T(t)`$ between q-ESS within a time window of fixed size $`T`$ as a function of time $`t`$ measured in number of generations. It is clear that $`\mathrm{\Omega }_T(t)`$ decreases with increasing $`t`$, however it is very difficult to obtain sufficient statistics to be able to determine the functional dependence of $`\mathrm{\Omega }_T(t)`$ on $`t`$, though a very slow exponential $`t`$ dependence is suggested by Fig. 5. Despite these sampling difficulties, it is evident that the duration of the q-ESS, on average, increases with time. This corresponds to a decrease in the extinction rate, consistent with analysis of the fossil record Newman and Sibani (1999).

### IV.2 Increasing population size, diversity and complexity

The gradual growth of the duration of the stable q-ESS epochs indicates that the dynamics of the system is able to produce more stable or better adapted configurations in genotype space. It is difficult to test quantitatively the stability of the q-ESS with respect to perturbations. That the population is distributed in an increasingly more efficient manner in genotype space can be seen directly from the increase in the total population size $`N(t)`$ averaged over an ensemble of different realizations of the stochastic elements of the dynamics. Fig. 6 contains the average total population $`N(t)`$ together with the ensemble average of the diversity $`D(t)`$, where $`D(t)`$ is defined as the number of different occupied positions $`𝐒`$ in genotype space at time $`t`$. The average diversity also increases with time.

Let us briefly consider how the total population size can increase. We saw in Sec. III that essentially $`p_{off}`$ is narrowly distributed either about $`p_{kill}`$ or, as in the q-ESS, about the value $`p_{qESS}`$ in Eq. (4). The increase in $`N(t)`$ is therefore not an effect of a gradual increase in $`p_{off}`$. Simulations indeed confirm that the average offspring probability always remains constant over the entire run. In biological observations and experiments, for reasons of uniqueness, the reproduction rate is identified as fitness. In this sense the fitness of the individuals remains, on average, constant in the TaNa model, as presumably is also the case in biological macroevolution; though the microbial experiments by Lenski Lenski and Travisano (1994) demonstrate that the reproductive fitness can increase as a result of adaptation in microevolution.

The increase in the average population size $`N(t)`$ observed in the TaNa model is caused by the system ability to generate configurations that increase the interaction term in the weight function $`H(𝐒,t)`$ defined in Eq. (1). When the first term increases, the second term $`\mu N(t)`$ in Eq. (1) can increase as well, while the total $`H(𝐒,t)`$ remains on, average, fixed.

The increase in the interaction term of $`H(𝐒,t)`$ is achieved in several ways. Firstly, the population is spread out onto an increasing number $`D(t)`$ of different genotypes, as seen in Fig. 6. Moreover, the evolutionary dynamics tends to produce occupied sites which are interacting with an increased number of other occupied sites, i.e., the number of non-zero terms $`J(𝐒^\alpha ,𝐒)n(𝐒)`$ in $`H(𝐒^\alpha ,t)`$ in Eq. (1) grows as the system produces configurations that are able to benefit better from the possible mutual interactions represented by $`J(𝐒^\alpha ,𝐒)`$. The distribution of active interaction links is shown at an early and a much later time in Fig. 7. We have not been able to resolve a shift with time in the distribution of the values of the active interaction strengths $`J(𝐒^a,𝐒^b)`$.

The increase in the diversity and the number of active links connected to an occupied position can be interpreted as an increase in the complexity of the configurations produced by the evolutionary dynamics. Selection and adaptation operate at the level of the entire configuration in genotype space rather than at the level of individual genotypes. This highlights that the biological concept of fitness makes most sense when considered as a collective property of an ecology, rather than an observable characteristic of the individual species or individual members of a population.

### IV.3 Record Statistics

The observation in the previous section that the dynamics of the TaNa model leads to an increase in a number of measurable quantities, taken together with the intermittent nature of the dynamics, suggests that the transitions between consecutive q-ESS epochs correspond to record transitions. One can imagine that some characteristic measures of the collective level of adaptation of the configurations generated in genotype space achieves an ever increasing value as the system undergoes a transition from one q-ESS to the next.

Sibani and collaborators Sibani et al. (1995); Sibani (1997); Sibani and Brandt (1998); Sibani and Pedersen (1999); Sibani and Littlewood (1993) have studied record dynamics and shown that the probability for $`n`$ records in a sequence of $`t`$ independently drawn random numbers is Poisson distributed on a logarithmic time scale, or equivalently: that the logarithm of the ratio of the time between the $`k`$-th and the $`(k1)`$-th record, $`\tau _k=\mathrm{ln}(t_k/t_{k1})`$, is exponentially distributed: $`P(\tau >x)=\mathrm{exp}(\lambda x)`$. Sibani and coworkers Sibani et al. (1995); Sibani (1997); Sibani and Brandt (1998); Sibani and Pedersen (1999) have also demonstrated the relevance of record statistics to the dynamics of the Kauffman NK model Kauffman and Levine (1987); Kauffman (1993).

Accordingly, it is interesting to investigate if the aging observed in the TaNa model exhibits signs of record statistics. To do this, we study the distribution of the variable $`\tau _k=\mathrm{ln}(t_k/t_{k1})`$ where $`t_k`$ denotes the time at which the $`k`$-th transition between consecutive q-ESS epochs occurs. We show in Fig. 8 that $`\tau _k`$ is exponentially distributed for large values and algebraically distributed for small values of $`\tau _k`$.

The exponential tail in Fig. 8 suggests that the transition times in the TaNa model follow record statistics in the region of large $`\tau `$ values, corresponding to the regime of long q-ESS durations. The algebraic form of $`p(\tau )`$ for small $`\tau `$ values indicates significant correlations for transitions that occurred in rapid succession. The question is which quantity evolves according to record dynamics. We have so far not been able to identify a variable of the system which jumps monotonously to ever higher values at the transition times. In the evolution models studied by Sibani et al. Sibani et al. (1995); Sibani (1997); Sibani and Brandt (1998); Sibani and Pedersen (1999) the fitness increases through consecutive records. As mentioned above in the TaNa model the reproductive fitness $`p_{off}`$ remains, on average, constant. The increase in the average duration of the q-ESS (see Fig. 5) suggests that the stability of the configurations in genotype space gradually increases. To explore this, one should study the temporal behavior of the eigenvalue spectrum of the stability matrix of the effective evolution equations in (3). We expect that the number of unstable directions, on average, decreases with time, though for a given realization fluctuations probably prevent a strictly monotonous behavior.

## V  Species abundance function

The fundamental quantity to describe an ecology is the species abundance function Hubbell (2001). The species abundance $`W(\rho )`$ is the ratio $`W`$ of species which contains a ratio $`\rho `$ of the total population. We plot in Fig. 9 the species abundance function for the TaNa model during a q-ESS, in this case we use the term species to denote individual positions in genotype space. A large number of positions are occupied by a small number of individuals; the occupancy of these positions is never established for extended periods during the q-ESS. The robust species contain a reasonable number of individuals and are distributed according to the broad peak. The peak can be fitted by a log-normal curve in a way similar to observed species abundance functions, see e.g. Hubbell (2001). We note that comparable species abundances functions are found in the predator-prey model studied by McKane, Alonso and Solé McKane et al. (2000).

## VI  Competition between sexual and asexual reproduction

We now consider the competition within the TaNa model between sexual and asexual reproduction in simulations of mixed populations. This is done by adding an extra reproduction determining gene to the genome (making $`L=21`$ in this section). This extra gene does not explicitly enter in $`H(𝐒,t)`$ but dictates an individual’s reproductive mode. Mutations to this gene occur during reproduction in the normal way. Obviously, the simplistic nature of the TaNa model excludes the model from realistically representing all biological features of the difference between sexual and asexual reproduction. Perhaps the most essential difference between sexual and asexual reproduction is the reshuffling of genes caused by the crossing over and recombination involved in sexual reproduction Burt (2000). It is possible that this effect can be captured by our definition of sexual reproduction for $`d_{max}1`$. To make the two reproductive modes as similar as possible in all aspects, except for the mixing of parent genes in the sexual case, we redefine in this section slightly the reproduction procedures in the following way. The only difference compared with the definitions in Sec. II.2 is that an asexual reproducing individual produces one new individual and we leave the parent in the system. In sexual reproduction the reproduction probability is assumed to be $`p_{off}`$ of the first of the two selected parents.

According to Weismann, sexual reproduction is more efficient to adapt to changed conditions because recombination produces a larger variation of types upon which Natural Selection is able to act, see Burt’s excelent review Burt (2000). It follows that in a harsher environment, sexual reproduction will be superior to asexual; a phenomenon which is observed Burt (2000). The TaNa model in its present form is not entirely able to account for these facts. It is informative, nevertheless, to investigate the behavior of the model and to understand why it fails in this particular respect.

In Fig. 10 we show the total size of the population together with the sizes of the sexual and asexually reproducing subpopulations. We see that the asexual population is always largest and even more so during the q-ESS.

In Fig. 11 we show the ratio of the average size $`N_a`$ of the asexual reproducing population and the average size $`N_s`$ of the sexually reproducing population. One notices that the asexual population is always more numerous, but less so for large $`d_{max}`$ and low $`\mu `$. That is the region where recombination is most important and the environment (modeled by $`\mu `$) least friendly. So in this respect the TaNa model exhibit a trend in the right direction in the sense of Weismann, see Burt’s review Burt (2000).

The reason the sexually reproducing population is unable to out-compete the asexual in the TaNa model is easy to understand. The offspring produced by two sexually reproducing individuals is likely to end up at a position in genotype space different from the parent positions. This effect makes it difficult to maintain the occupancy of the parent positions, even if these positions are highly fit. It is obviously this variation that is expected to make sexual reproduction more efficient in searching genotype space. However, in its present form, the TaNa model allows very large variations in the weight function $`H`$ in Eq. (1), even for positions very close by in genotype space. Whether this is reasonable or not depends on the interpretation of the genomic sequence $`𝐒`$. If we think of each “gene” $`S_i`$ as representing, in an averaged way, a section of length $`1/L`$ of the entire genome, then a separation by Hamming distance one is very significant for $`L`$ as small as 20; and it is reasonable that two positions in genotype space of this or larger separation may have very different weight functions $`H`$. On the other hand, this interpretation means that two individuals being different by a single Hamming distance differ by a major fraction of the entire genome and should therefore not belong to the same species and accordingly not be able to produce any (viable) offspring. We conclude that a more realistic study of the competition between sexual and asexual reproduction by use of the TaNa model should be possible if a more slowly varying weight function is used. Future studies will focus on this problem.

## VII  The Error Threshold

At sufficiently large mutation rates $`p_{mut}`$, offspring are so different from their parents that the occupation in genotype space rapidly moves from one position to the next. When this happens, it becomes impossible to establish the q-ESS seen in Fig. 2 and consequently the entire simulation consists of one hectic period. The change from the behavior depicted in Fig. 2, where the hectic periods are of much shorter duration than the q-ESS, to the behavior where the q-ESS are absent, occurs over a very narrow region of $`p_{mut}`$ values. Considering first large values of $`p_{mut}`$ we gradually decrease $`p_{mut}`$ in the simulations and we identify the threshold value, $`p_{th}`$ of $`p_{mut}`$ at which q-ESS are observed as the error threshold Eigen et al. (1988); Drossel (2001). In Fig. 12 we plot the simulated value of $`p_{th}`$ for different values of the width parameter $`c`$.

We can estimate the $`c`$ dependence of $`p_{th}`$ by the following argument. From Fig. 4 we know that the distribution of $`p_{off}`$ in the hectic periods is centered about $`p_{kill}`$ and in the q-ESS is centered about $`p_{qESS}`$ defined in Eq. (4). Changing the parameter $`c`$ will change the width of the distribution of the $`p_{off}`$ values (see Eqs. (1) and Eq. (2)). It will be possible to establish q-ESS in between the hectic periods if $`p_{kill}+\sigma _pp_{qESS}`$, where $`\sigma _p`$ is the half width of the peak in the distribution of $`p_{off}`$ in the hectic periods. We translate this argument to the distribution of the $`H`$ values and obtain the following estimate for $`p_{th}`$:

$$p_{th}=12^{1/L}[(1p_{kill})e^{\alpha /c}+1+p_{kill}]^{1/L}.$$
(5)
Here we have assumed that the width $`\sigma _p`$ of the peak in the distribution of $`H`$ values will be given by $`\sigma _p=\alpha /c`$ (see Eq.(1)) in which case $`\alpha `$ is a measure of the standard deviation of the factor in Eq. (1) multiplying $`1/c`$. We have used $`\alpha =0.07`$ to fit the simulation data in Fig. 12. This value is somewhat larger but of the right order of magnitude as the corresponding quantity measured during the simulation.

## VIII  Parameter dependence

For completeness, we present here the dependence on the parameters $`c`$ and $`\mu `$ in the Hamiltonian given in Eq. (1).

We show in Figs. 13 and 14 the averaged occupation measured as the ratio between the average number of individuals and the average number of occupied positions in genotype space for purely asexual and sexual reproducing populations respectively. As expected, the system is able to support the largest populations in the region of small $`\mu `$ parameter and broad distribution of coupling strength, i.e., small values of $`c`$. The sexual reproduction is most sensitive to a decrease in the carrying capacity (increase in $`\mu `$) or a decrease in the width of the range of possible $`J(𝐒^a,𝐒^b)`$ couplings (increase in $`c`$).

## IX  Discussion and Conclusion

The Tangled Nature model may be considered as a mathematical framework for the study of evolutionary ecology. The dynamics of the model is defined at the level of individuals, either as asexual or as sexually reproducing individuals. All ecological structure in the model arises through emergence. The model is able to generate many of the observed features of biological evolution starting from a basic implementation of the key assumptions in the Darwinian evolution paradigm.

The density of individuals in genotype space segregates corresponding to the emergence of distinct species. The interaction between individuals gives rise to a jerky or intermittent macro-dynamics, in which quasi-stable configurations (the q-ESS) in genotype space are abruptly replaced by new quasi-stable configurations. This mode of operation can be compared with the intermittent behavior observed in the fossil record and emphasized be Gould and Eldredge’s the term “punctuated equilibrium” Gould and Eldredge (1977). The TaNa model is always in a transient where the configurations generated as a result of adaptation to the co-evolutionary selective pressure gradually produce configurations, or ecologies, in genotype space which collectively exhibit a higher degree of adaptation, in the sense that the average lifetime of these q-ESS increases slowly. This behavior compares well with the observation that the fossil record indicates a decrease in the extinction rate Newman and Sibani (1999). The increase in the lifetime of the q-ESS is associated with an increase in the complexity (species diversity, number of active interactions) of successive configurations. This gradual aging together with the intermittent nature of the dynamics suggest that some characteristic of the evolving ecosystem might be undergoing record statistics in the sense of Sibani and co-workers Sibani et al. (1995); Sibani (1997); Sibani and Brandt (1998); Sibani and Pedersen (1999). So far we have unfortunately not been able to identify the appropriate effective variable which moves through the records but we expect this variable to be related to the stability matrix of the effective dynamical equation.

The species abundance distribution generated by the TaNa model encourage future studies of larger populations with longer genome sequences. This will enable a hierarchical study of the taxonomic organisation of the generated ecologies. Using distance criteria in genotype space one can study the clustering of individuals into species, of species into genera, etc. Future studies will also examine the phylogenetic structures in detail, especially during the radiation of species encountered in the transition periods between q-ESS. Using longer genome sequences and a more smoothly varying weight function, we expect the TaNa model to be able to illuminate the evolutionary competition between sexual and asexual reproduction.

###### Acknowledgements.

It is a pleasure to acknowledge very helpful discussions with Paolo Sibani. MH and SAC are supported by EPSRC studentships and KC gratefully acknowledges the financial support of EPSRC trough Grants No. GR/R44683/01 and GR/L95267/01. We thank Paul Anderson for reading the manuscript.