text:
# 1 Introduction

## 1 Introduction

There has recently been an upsurge of interest in so-called scale-free networks, where the vertex order, or degree of connectivity per node $`k`$, follows a power-law distribution

$$P\left(k\right)k^\gamma $$
(1)
for $`k>k`$. This is in contrast to the exponential suppression expected from randomly wired networks. For a number of real-world networks like social networks, powergrids, citation networks, the World Wide Web and metabolic networks , the scale-free behaviour of Eq. (1) has been observed with exponents $`\gamma `$ in the range 2-4 . The nature of all examples above is that they originate from a growing process. Hence, it appears natural to develop a model of growing network, which was indeed done in .

The degree of connectivity is a measure of local properties of the networks. Less attention has been payed to more global properties for the above examples except for general features like network diameters. Global measures could be important discriminants of network properties, in particular when it comes to narrow in on potential underlying models.

In this paper, we analyze topological properties of networks, both of synthetic and real-world nature, by extracting spectral dimensions $`d_S`$ using a random walk procedure, from which the return-to-origin probability is estimated. The focus is on scale-fee networks; citation networks, metabolic networks and the growing network model  are investigated. For completeness and comparisons, we also extract the $`d_S`$ from simple three-dimensional regular lattice networks, random networks and small-world networks .

We find that the citation and metabolic networks are quite low-dimensional with $`d_S`$ around $`3`$, whereas for the growing network model $`d_S`$ is approximately $`7.5`$. The latter turns out to be more in parity with what characterizes random and small-world networks. As a consistency check, the dimension of regular lattice networks is also determined, which comes out as expected. One concludes that the spectral dimension offers a powerful and additional measure to $`\gamma `$ in Eq. (1) when it comes to characterizing network topologies. Furthermore, the diffusion process underlying the extraction of $`d_S`$ hints upon differences in the citation and metabolic network topologies; the former is tree-like whereas the latter is rich in loop structures.

This paper is organized as follows: In Sect. 2 we describe the method for extracting the spectral dimension and Sect. 3, 4 and 5 contain our investigations of the synthetic, citation and metabolic networks respectively. A summary can be found in Sect. 6.

## 2 Spectral Dimension

Our method to probe the topological properties of an interaction network, which is of more global nature than the degree of connectivity, is based upon the diffusion of a test particle in the metric space defined by the graph. In a continuous space with a fixed smooth metric, the diffusion equation has the form

$$\frac{}{t}K_g(ϵ,ϵ_0,t)=\mathrm{\Delta }_gK_g(ϵ,ϵ_0,t).$$
(2)
Here $`t`$ is the diffusion time, $`\mathrm{\Delta }_g`$ is the Laplace operator in the metric $`g`$ and $`K_g(ϵ,ϵ_0,t)`$ is the probability density to diffuse from $`ϵ_0`$ to $`ϵ`$ in time $`t`$. For small $`t`$ it is well known , that the average return probability has the following asymptotic expansion

$$K_g(0,0,t)t^{d_S/2}$$
(3)
with

$$K_g(ϵ,0,0)=\delta \left(ϵ\right)$$
(4)
and where the power $`t^{d_S/2}`$ reflects the dimension of the network.

In the spirit of these equations we extract the spectral dimension of the geometry defined by the interaction networks considered in this work. To this end, we use the transition matrix

$$C_{ij}=\frac{J_{ij}}{k_i}$$
(5)
with

$$J_{ij}=\{\begin{array}{cc}1\hfill & \text{if }i\text{ and }j\text{ are neighbours}\hfill \\ 0\hfill & \text{otherwise}\hfill \end{array}$$
(6)
as the discrete version of the Laplace operator $`\mathrm{\Delta }_g`$. The number $`k_i`$, counts the number of links connected to vertex $`i`$ – the vertex-order. In a simulation of the diffusion process defined by the transition matrix, the probability to return, $`P`$, after $`t`$ steps to the origin $`\delta _i`$ is then measured. More precisely, we choose a random subset $`\left\{\delta _i\right\},i=1\mathrm{}N`$ of vertices and extract

$$<P\left(t\right)>=\frac{1}{N}\underset{i}{\overset{N}{}}\delta _i\left|C_{i,j}^T\right|\delta _i.$$
(7)
We then fit the resulting distribution to the asymptotic form in eq. (3). For large $`t`$, Eq. (3) is dominated by the eigenvector of $`C_{ij}`$ with eigenvalue $`\lambda =1`$; the diffused particles reach an equilibrium distribution and Eq. (3) does not hold. Also, for too small $`t`$ the assumption of a smooth metric is not justified. Our extraction procedure is therefore somewhat more elaborate to account for these effects. A sliding window method is used, where the window is chosen such as to minimize the standard deviation per point used in the fit. In other words, we fit to the part of the distribution, which is closest to the functional form assumed in this procedure. Depending upon $`d_S`$, window sizes might vary from a few time points to the entire range.

## 3 Synthetic Networks

Next we extract $`d_S`$ for a few different types of synthetic networks with sizes up to around 30000 nodes in order to compare with typical sizes of citation networks. In addition, for all networks the average geodesic distance (diameter) $`R`$ is computed.

Lattice Networks. We use a simple square lattices with dimension $`d=3`$. For these it is only possible to return to the origin after an even number of steps. Therefore $`P\left(t\right)=0`$ for odd $`t`$. In terms of the transfer matrix eigenvectors this means that the symmetry of the square lattice leads to an eigenvector with $`\lambda `$ exactly $`1`$. The corresponding eigenstate does not decay for large $`t`$ and yields destructive interference with the $`\lambda =1`$ eigenstate after an odd number of steps. When extracting the spectral dimension $`d_S`$ from the data we therefore omit the odd time-points. As can be seen from Fig. 1, $`d_S`$ is slightly above the expected value $`3`$, the dimension of the lattice. We have repeated this experiment for a regular lattice including diagonals, hence allowing for closed loops with an odd number of steps. Although local properties of the graph are considerably changed since the number of links emerging from each vertex is increased from $`2d`$ to $`2d+d\left(d1\right)`$, we still extract the same value for $`d_S`$ within errors. This is a further justification of omitting the data at odd time points for the square lattice and emphasizes that the spectral dimension is sensitive to global geometric aspects rather than to local details.

Small-world Networks. We generate small-world networks  from the lattice network described above by rewiring the edges in the regular lattice to a randomly chosen vertex with probability $`p`$= $`0.01`$, $`0.05`$ and $`0.2`$ respectively. This generates three sets of models for each system size. Not surprisingly, these networks end up with spectral dimensions (see Fig. 1) in between the regular lattice and the homogeneous random networks described next.

Random Networks. Homogeneous random networks are generated as directed networks, where each vertex has two inputs from other vertices choosen at random. In the analysis the graph is considered undirected by ignoring the orientation on the edges connecting the vertices. Furthermore, the analysis is focused on the largest cluster; the largest connected part of the resulting graph. The spectral dimension $`d_S8.5`$ (see Fig. 1) is relatively large. This agrees qualitatively with the observation that a characteristic length-scale, e.g. average distance between vertices along links, grows slowly with the volume. It also signals a large dimensionality for this type of graph. It is interesting to note that again a distinct even-odd disparity emerges; the amplitude after an odd number of steps is almost zero. This indicates that the random network is dominated by a tree-like structure, the number of loops (at least of odd length) is negligibly small compared to the amplitude obtained from backtracking the same path.

Growing Network Model. In this model  one at each time step $`t`$ adds a node, which connects to $`m`$ existing nodes ($`i`$) with a probability $`\mathrm{\Pi }_i`$ given by

$$\mathrm{\Pi }_i=\frac{k_i^\beta }{_jk_j^\beta }$$
(8)
where $`k_i`$ is the connectivity of node $`i`$ and $`\beta `$ is a parameter. In  the $`\beta `$=1 case was investigated both analytically in the $`t\mathrm{}`$ case and numerically for finite $`t`$. With the approximations involved one arrived at $`\gamma =3`$ in the analytic approach and $`\gamma =2.9\pm 0.1`$ from simulations. For these networks, $`R`$ is approximately $`6`$ for large system sizes.

## 4 Citation Network

In a citation network the links and nodes correspond to the citations and the publications (citing and cited) respectively. We use the SPIRES data base  for our citation network studies. This data base is limited to high energy physics publications but is not confined to articles that have been published in refereed journals – citations to and from conference reports are also present. From the data we construct citation networks, which we treat as undirected, not fully connected graphs. By considering graphs generated by publications in a certain year or time span and the cited papers (regardless of publication year) we obtain a whole set of graphs. In these not fully connected graphs we focus on the largest connected cluster, the corresponding sizes are shown in Table 1. When computing the connectivity distribution of the nodes, we confirm the power law distribution (Cf. Eq. (1)) already observed in citation networks  with $`\gamma =2.7`$.

In Fig. 2 we plot the probability to return to the origin after $`t`$ time steps for a few typical networks. It is interesting to note that, within errors, the return probability for the 76 and 79 networks are in very good agreement indicating universal topological properties. This holds true also when comparing all of the other one-year networks with the exception of the $`75`$ and the $`87`$ networks (not shown), which cannot be fitted to Eq. (3). In these years the connectivity does probably not generate a homogeneous metric, because the geometry is dominated by two or more large clusters, which are interconnected by only a few links.

For the networks composed by several years we again observe a remarkable universality when comparing different time-periods. However, the distribution observed is considerably different from the one-year distributions. The reason presumably is a maturation of the citation networks leading to a tighter interconnection of the central cluster. In turn we observe larger dimension as can be seen from the steeper slope of the corresponding graph in Fig. 2. These issues for the citation networks have been subject to further investigations .

## 5 Metabolic Network

Another real-world interaction network is the metabolic network found in living cells. In these networks, substrates are treated as vertices, while chemical reactions connecting substrates and educts are treated as directed links. Recently it has been demonstrated that for this type of networks the connectivity distribution obeys a power-law behavior (Cf. Eq. (1)); these networks are scale-free with respect to the order distribution .

Using data from the EMP-Project  we constructed a network including all reactions found in the database without taking into account species and cell locations. As above, we neglect the orientation of the resulting graph in the analysis. The spectral dimension observed for this network is $`d_S2.8`$. This is a surprisingly small dimension taking into account the observation  that the average distance between vertices on this graph does not grow with the graph-size, implying a very large, possibly infinite dimension. Also, comparing $`d_S`$ with the corresponding value for the growing model scale-free network analyzed above, the difference is remarkable. This indicates that, although the metabolic network and the scale-free network are similar on a local (vertex order) scale, the more global topological properties are very different. Another important difference arises from the observation that the return probability for odd path-lengths do not vanish. This means that the metabolic networks consist of a very large number of closed loops. The graph has, in contrast to the random, growing model and citation networks, not a dominantly tree-like structure. We interprete this as an indication of a built-in stability of metabolic networks with respect to small modifications.

## 6 Summary and Outlook

In table 2 the results for $`d_S`$ are summarized for the different networks. Also shown here are the average (geodesic) distances $`R`$. For the regular lattice, the spectral dimension reproduces quite well the dimension of the underlying lattice. Depending upon the probability to rewire a connection ($`p`$), the small-world networks takes on values between the regular lattice ($`p=0`$) and the random network ($`p=1`$), which has the largest spectral dimension of all networks probed here. One might have guessed that the growing network has an even larger dimension – the average distance for a network with the same size is smaller. However, in contrast to the geodesic distance, which only counts the shortest paths, the spectral dimension takes into account all possible paths of a given length. This means that only a few links, in extreme cases even a single link, can considerably change the geodesic distance, while the spectral dimension probes larger parts of the geometry and is therefore only slightly affected. For example, consider the behavior of the average distance for small-world networks with different probabilities $`p`$ to rewire the underlying three dimensional regular structure (table 2). Even for $`p=0.01`$ with the resulting geometry still essentially a regular three dimensional lattice, the geodesic distance jumps to less than one half of the value for the regular case, while the spectral dimension is changed only moderately.

The spectral dimensions for the real-world networks, the citation and metabolic networks respectively, are strikingly similar in view of the much larger values observed for the network growing model and the random network. We interprete this as a sign of universality in the geometry of the real-world networks. The difference in $`d_S`$ when comparing to that of the network growing model is remarkable. As mentioned in the introduction, the latter reproduces quite nicely the power-law scale-free vertex order distribution observed for the real-world networks . However, the difference in $`d_S`$ observed here indicates that the geometries of real-world networks are not fully described by the network growing model.

Given the similarity in $`d_S`$ for the citation and metabolic networks, one might ask how similar they really are? An important difference can be observed by investigating the probability $`P_{odd}`$ to return to the origin after an odd number of steps. For the citation network, as for all the synthetic networks, one obtains $`P_{odd}0`$. This is in contrast to the metabolic network for which $`P_{odd}P_{even}`$. The case $`P_{odd}0`$ indicates a tree like geometrical structure with only few loops, where the return probability is dominated by the inverse of the forward paths. The regular lattice without diagonals seems to be a counter example. The return probability $`P_{odd}`$ is exactly zero, but the geometry is certainly not tree-like and there exist many non-backtracking loops. However, this lattice is a special case since its symmetry only allows for even step non-backtracking loops. It is extremely unlikely that such an effect plays a role in any of the other lattices.

Our results can be summarized as follows:

* The spectral dimension $`d_S`$, which reflects a diffusion process on networks, offers quite some promise in categorizing networks beyond what emerges from studying vertex order distributions.
* In particular, we find that $`d_S`$ are very similar for citation and metabolic networks. This may indicate a universal behaviour of real-world networks.
* With respect to $`d_S`$, the growing network model is significantly different from citation and metabolic networks.
* As a by-product when extracting $`d_S`$, we find that whereas citation networks have tree-like structures the metabolic networks appear to contain many loops.

In addition, lattice, random and small-world networks were probed using $`d_S`$ with results exhibiting internal consistency of the method.

When finishing the write-up of this work, a paper with somewhat similar scope and philosophy but with employing eigenvalue analysis techniques instead was released; I.J. Farkas et. al., ”Spectra of ”Real-World” Graphs: Beyond the Semi-Circle Law” \[cond-mat/0102335\].

Acknowledgements:

We are very much indebted to Louise Addis in the SPIRES-HEP data base group at the Stanford Linear Accelerator Center library for providing us with the SPIRES data base in a form that substantially facilitated this project. Alex Selkov from the EMP-Project supported us in acquiring the data for the metabolic network. We have also benefitted from discussions with Bo Söderberg regarding eigenvalue properties of network diffusion. This work was in part supported by the Swedish Foundation for Strategic Research and the Knut and Alice Wallenberg Foundation through the SWEGENE consortium.


format: text

text:
# Generalized Crossover in multiparameter Hamiltonians

## A  The Large $`𝐍`$ Limit

A rather explicit illustration of the mechanism described in the present paper is obtained by considering the tetragonal model in the limit of an infinite number of field components ($`N\mathrm{}`$). At variance with standard $`O(N)`$ vector models, the tetragonal model does not became trivial on this limit, because nontrivial contributions to all orders of $`\overline{v}`$ and $`\overline{w}`$ couplings are still present. However some semplifications occur which make our discussion, while stile quite general, formally much simpler.

In the large $`N`$ limit is possible to show that the RG functions take the following form:

$`\beta _{\overline{u}}(\overline{u},\overline{v},\overline{w})`$ $`=`$ $`A(\overline{v},\overline{w})\overline{u}B(\overline{v},\overline{w})\overline{u}^2`$ (A1)
$`\beta _{\overline{v}}(\overline{u},\overline{v},\overline{w})`$ $`=`$ $`\stackrel{~}{\beta }_{\overline{v}}(\overline{v},\overline{w})`$ (A2)
$`\beta _{\overline{w}}(\overline{u},\overline{v},\overline{w})`$ $`=`$ $`\stackrel{~}{\beta }_{\overline{w}}(\overline{v},\overline{w})`$ (A3)
$`\eta _\varphi (\overline{u},\overline{v},\overline{w})`$ $`=`$ $`\stackrel{~}{\eta }_\varphi (\overline{v},\overline{w})`$ (A4)
$`\eta _t(\overline{u},\overline{v},\overline{w})`$ $`=`$ $`\stackrel{~}{\eta }_t(\overline{v},\overline{w})+B(\overline{v},\overline{w})\overline{u}`$ (A5)
The system of equations $`\stackrel{~}{\beta }_{\overline{v}}(\overline{v},\overline{w})=0`$ and $`\stackrel{~}{\beta }_{\overline{w}}(\overline{v},\overline{w})=0`$ admits four sets of solutions $`(\overline{v}^{},\overline{w}^{})`$. For each set one finds two fixed points, corresponding to the values $`\overline{u}^{}=0`$ and $`\overline{u}^{}={\displaystyle \frac{A(\overline{v}^{},\overline{w}^{})}{B(\overline{v}^{},\overline{w}^{})}}`$.

Because of the above relationships, the differential equation satisfied by the function $`Z_t^1(\overline{u},\overline{v},\overline{w})`$ can be solved in the large $`N`$ limit by the Ansatz

$`Z_t^1(\overline{u},\overline{v},\overline{w})=\stackrel{~}{Z}_t^1(\overline{v},\overline{w})\left[1\overline{u}Y(\overline{v},\overline{w})\right]`$ (A6)
leading to the equations:

$`\left[\stackrel{~}{\beta }_{\overline{v}}(\overline{v},\overline{w}){\displaystyle \frac{}{\overline{v}}}+\stackrel{~}{\beta }_{\overline{w}}(\overline{v},\overline{w}){\displaystyle \frac{}{\overline{w}}}+\stackrel{~}{\eta }_t(\overline{v},\overline{w})\right]\stackrel{~}{Z}_t^1(\overline{v},\overline{w})=0`$ (A7)
$`\left[\stackrel{~}{\beta }_{\overline{v}}(\overline{v},\overline{w}){\displaystyle \frac{}{\overline{v}}}+\stackrel{~}{\beta }_{\overline{w}}(\overline{v},\overline{w}){\displaystyle \frac{}{\overline{w}}}+A(\overline{v},\overline{w})\right]Y(\overline{v},\overline{w})=`$ (A8)
$`B(\overline{v},\overline{w}).`$ (A9)
The first equation looks like eq.(6) and it is simply the restriction of the evolution to the $`u_0=0`$ plane; we can repeat our general arguments, obtaining in particular the RG trajectory that connects the $`I`$ and $`XY`$ fixed points. Notice however that, since the condition $`\stackrel{~}{Z}_t^1(\overline{v},\overline{w})=0`$ is independent of $`\overline{u}`$, it defines a surface in full parameter space, and the above discussion shows that the fixed points $`C`$ and $`T`$ must lay on this surface.

Once the functions $`\stackrel{~}{Z}_t^1`$ and $`Y(\overline{v},\overline{w})`$ have been determined, it is easy in the large $`N`$ limit to reconstruct the full $`Z_t^1=0`$ surface, which can be simply described by the above condition and by the function:

$`\overline{u}(\overline{v},\overline{w})={\displaystyle \frac{1}{Y(\overline{v},\overline{w})}}`$ (A10)
Notice that, as a consequence of eq.(A9), the function $`\overline{u}(\overline{v},\overline{w})`$ does not depend on the detailed form of the RG function $`\stackrel{~}{\eta }_t`$, as expected from our general arguments.

The intersections of the crossover surface with the planes $`\overline{v}=0`$ and $`\overline{w}=0`$ can now be found in a rather explicit form, by exploiting the above simplifications.

In terms of the generic variable $`z`$ we obtain the relevant equations:

$`\left[\stackrel{~}{\beta }(z){\displaystyle \frac{}{z}}+\stackrel{~}{\eta }_t(z)\right]\stackrel{~}{Z}_t^1(z)=0`$ (A11)
$`\left[\stackrel{~}{\beta }(z){\displaystyle \frac{}{z}}+A(z)\right]\left({\displaystyle \frac{1}{\overline{u}(z)}}\right)=B(z)`$ (A12)
It is straight-forward to solve the linear equations, obtaining

$`\stackrel{~}{Z}_t^1(z)=exp\left[{\displaystyle _0^z}{\displaystyle \frac{\stackrel{~}{\eta }_t(z^{})}{\stackrel{~}{\beta }(z^{})}}𝑑z^{}\right]`$ (A13)
$`\overline{u}(z)={\displaystyle \frac{X(z)}{{\displaystyle _0^z}{\displaystyle \frac{B(z^{})}{\stackrel{~}{\beta }(z^{})}}X(z^{})𝑑z^{}}}`$ (A14)
$`\text{where}X(z)=exp\left[{\displaystyle _0^z}{\displaystyle \frac{A(z^{})}{\stackrel{~}{\beta }(z^{})}}𝑑z^{}\right]`$ (A15)
It is easy to check that $`\overline{u}(z)`$ is a RG trajectory and that in the limits $`z0`$, $`zz^{}`$ ($`\stackrel{~}{\beta }(z^{})=0`$) we have $`\overline{u}(0)=\frac{A(0)}{B(0)}`$ and $`\overline{u}(z^{})=\frac{A(z^{})}{B(z^{})}`$ respectively, consistent with the boundary conditions at the fixed points.

The above expressions lend themselves to simple analytical integration in the one-loop approximation and to easy numerical integration in the more general case.

For the sake of illustration we computed explicitly the crossover lines on the $`(u,w)`$ plane. Fig.1 shows the results of our numerical integration of the equations, when resummed six-loop RG functions are employed.

The straight line connecting $`I`$ to $`C`$ is the intersection of the $`\overline{v}=0`$ plane with the $`\stackrel{~}{Z}_t^1=0`$ surface.

## B An effective one-loop model

It is possible to find a set of one-loop truncated RG equations which retains all the essential features of the above described model and lends itself to analytic integration, thus offering a rather explicit illustration of our results.

Our starting point is the set of equations:

$`\beta _u=u+a(v+w)u+u^2,`$ $`\beta _v=v+v^2`$ (B1)
$`\beta _w=w+2vw+2w^2,`$ $`\eta _t=ub(v+w)`$ (B2)
where some trivial rescaling of the couplings is understood.

The properties and symmetries of the fixed points are the same as in our general discussion, and in particular when $`1a2`$ also the stability properties are exactly reproduced. By direct integration one can find that:

$`\stackrel{~}{Z}_t^1(v,w)=(1v)^{\frac{b}{2}}(1v2w)^{\frac{b}{2}}`$ (B3)
and, being $`F`$ the standard hypergeometric function:

$`Y(v,w)=`$ (B4)
$`(1)^{\frac{a}{2}}w^{1a}(v^2+2vw)^{\frac{a}{2}1}[(1+{\displaystyle \frac{v}{w}})F({\displaystyle \frac{1}{2}},{\displaystyle \frac{a}{2}},{\displaystyle \frac{3}{2}},(1+{\displaystyle \frac{v}{w}})^2)`$ (B5)
$`(12v+{\displaystyle \frac{v}{w}}(1v))F({\displaystyle \frac{1}{2}},{\displaystyle \frac{a}{2}},{\displaystyle \frac{3}{2}},(12v+{\displaystyle \frac{v}{w}}(1v))^2)].`$ (B6)
One may easily check that the surface described by the function $`u(v,w)=Y^1(v,w)`$ and the condition $`\stackrel{~}{Z}_t^1(v,w)=0`$ enjoys all the properties discussed in the text.


text:
# 1 Coulomb contribution to pion dissociation into two jets.

1.  We study the diffractive dissociation of a pion on a nucleus into two jets. In this process a highly energetic pion interacts with a nucleus and produces two jets (di-jets) which in the lowest approximation are formed from the quark-antiquark ($`q\overline{q}`$) pair.

$$\pi Aq\overline{q}A$$
(1)
Early studies of this process can be found in Ref. . Recently it was measured in the E791 experiment  and was the subject of further theoretical investigations . The process is of interest for its potential to measure directly the pion distribution amplitude (i.e. the probability amplitude to find in a pion a valence $`q\overline{q}`$ Fock state with a certain momentum fraction) and to study the effects of colour transparency.

The main contribution to the diffractive process (1) is due to Pomeron exchange, or in QCD language, by the colour singlet gluon ladder. However, this process can also occur as result of the electromagnetic interactions between pion and target nucleus (Coulomb exchange). The strenght of the electromagnetic coupling is $`\alpha Z`$ ($`\alpha `$ is the electromagnetic fine coupling constant). For heavy nuclei $`\alpha Z`$ is not small and therefore one can ask about the size of the Coulomb contribution to (1). The answer to this question is the subject of the present paper.

2.  The cross section for this process is largest when the momentum transfer to the nucleus is smaller than the scale of the nucleus form factor $`\mathrm{\Lambda }60`$MeV (for Platinum, $`Z=78`$). Therefore we consider the process (1) in the case that large transverse momenta of the quark jet ($`q_1`$) and of the antiquark one ($`q_2`$) balance each other, $`|q_1+q_2||q_1|,|q_2|`$. The diagrams describing the Coulomb contribution to the process (1) are shown in Fig. 1. The wavy line denotes photon exchanged in the t-channel. The large transverse momenta of di-jets result from the hard gluon exchange denoted by the dashed line. They supply a hard scale to the process and therefore we neglect the pion mass. The incoming pion has momentum $`p_1,p_1^2=0`$. The nucleus mass is $`M`$, its momenta in the final and in the initial states are $`p_2`$ and $`p_3`$, respectively. The large cms energy squared of the process is given by the Mandelstam variable $`s=(p_1+p_2)^2`$, whereas the small momentum transfer squared equal to the photon virtuality $`k^2`$ is denoted as $`t=k^2=(p_2p_3)^2`$.

Let us also introduce the auxiliary Sudakov vector $`p_2^{}`$ such that $`(p_2^{})^2=0`$

$$p_2=p_2^{}+\frac{M^2}{\overline{s}}p_1,$$
(2)
where $`\overline{s}=2(p_1p_2^{})=sM^2`$. Then, the Sudakov decompositions of the on-shell $`q\overline{q}`$ momenta in the di-jets are

$$q_1=zp_1+\frac{q_1^2}{z\overline{s}}p_2^{}+q_1\text{and}q_2=\overline{z}p_1+\frac{q_2^2}{\overline{z}\overline{s}}p_2^{}+q_2.$$
(3)
Here, the variable $`z`$ describes the fraction of pion momentum carried by quark, $`0z1`$. For the corresponding fraction of the antiquark momentum we use the shorthand notation $`\overline{z}=1z`$. The decomposition of the photon momentum $`k`$ reads

$$k=\alpha p_1+\beta p_2^{}+k_{},$$
(4)
its parameters $`\alpha `$, $`\beta `$ and the transverse momentum $`k_{}`$ are expressed through jets variables

$$\beta =\frac{M_{2j}^2+(q_1+q_2)^2}{\overline{s}},\alpha =\frac{(q_1+q_2)^2+\frac{M^2\beta }{\overline{s}}}{\overline{s}},k_{}=q_1+q_2,$$
(5)
here $`M_{2j}^2`$ is the invariant mass of the di-jets

$$M_{2j}^2=\frac{q_1^2}{z}+\frac{q_2^2}{\overline{z}}(q_1+q_2)^2.$$
(6)
As already mentioned the main contribution to the process comes from the region of small momentum transfer (or photon virtuality $`k^2`$) $`t=k^2=k_{}^2+\alpha \beta \overline{s}`$. There $`k_{}<<q_1,q_2`$, which means that $`q_2q_1`$ and

$$M_{2j}^2=\frac{q_1^2}{z\overline{z}}.$$
(7)
The momentum transfer $`t`$ equals $`t=k^2(k_{}^2+k_{min}^2)`$, and according to Eq.(5) its minimal value is

$$k_{min}^2=\frac{M^2M_{2j}^4}{\overline{s}^2}=\frac{M_{2j}^4}{4E^2},$$
(8)
the last equality being valid in the rest frame of the nucleus in which the pion energy equals $`E`$.

3.  Our aim is to calculate the leading asymptotics of the scattering amplitude $`M`$ in powers of $`1/q_1^2`$ (at the leading twist level). It is given as the convolution of the hard scattering amplitude $`T_H(u,\mu _F^2)`$ with the pion light-cone distribution amplitude $`\varphi _\pi (u,\mu _F^2)`$

$$M=\underset{0}{\overset{1}{}}𝑑u\varphi _\pi (u,\mu _F^2)T_H(u,\mu _F^2)$$
(9)
where the hard scattering amplitude $`T_H(u,\mu _F^2)`$ describes the production of a free $`q\overline{q}`$ pair (the di-jets) in collision of the $`t`$-channel photon with $`q`$ and $`\overline{q}`$ (the later having momenta $`up_1`$ and $`\overline{u}p_1`$, respectively), collinear to the pion momentum $`p_1`$. The factorization scale $`\mu _F`$ is of the order of the di-jets transverse momentum $`q_1`$. At lowest order in $`\alpha _s`$, which we consider in this paper, $`T_H(u,\mu _F^2)`$ does not depend on $`\mu _F`$.

Eq.(9) describes the factorization procedure in QCD which disentangles the contributions to $`M`$ coming from large and small distances. The soft part is described by the pion distribution amplitude , which is defined as follows

$$0|\overline{d}(x)\gamma _\mu \gamma _5u(x)|\pi ^+(p)_{x^20}=ip_\mu f_\pi \underset{0}{\overset{1}{}}𝑑ue^{i(2u1)(xp)}\varphi _\pi (u).$$
(10)
The constant $`f_\pi =131MeV`$ is known experimentally from $`\pi \mu \nu `$ decay,
$`0|\overline{d}(0)\gamma _\mu \gamma _5u(0)|\pi ^+(p)=ip_\mu f_\pi `$.

The hard part, i.e. the amplitude $`T_H(u)`$ in Eq.(9), is calculable in the perturbative QCD, it is given by four tree diagrams shown in Fig. 1.

The factorization procedure described by Eq. (9) is similar to that used for the hard exclusive processes . Its validity in the case of our quasi-exclusive process can be clarified as follows. Let us consider, for example, the diagram Fig.1(b). The large transverse momentum of quark jets flows along the lines $`ABC`$. Therefore their virtualities are much larger than those of the other quark lines $`DA`$ and $`DB`$. Thus, at leading twist, the quark lines $`DA`$ and $`DB`$ have to be considered as being on mass shell and this part of the diagram can be factorized out of the hard part given by the highly virtual quark and gluon propagators.

Calculating the hard amplitude $`T_H(u)`$ in the Feynman gauge we perfom the usual substitution in the nominator of the photon propagator

$$g^{\mu \nu }2\frac{p_1^\mu p_2^\nu }{\overline{s}},$$
(11)
where $`p_2^\nu `$ acts on the upper quark vertex. This brings an inaccuracy $`1/\overline{s}`$ which at high energies is very small.

Since the energy of a photon $`k`$ is small the nucleus can be considered as a ”scalar particle” with the electromagnetic formfactor $`F_{QED}(k^2)`$, so the photon-nucleus vertex leads to the factor (see Eq.(11))

$$(ieZ)(p_2+p_3)_\mu e^\mu F_{QED}(k^2)=(ieZ)(p_2+p_3)_\mu \frac{p_1^\mu }{\overline{s}}F_{QED}(k^2)=(ieZ)F_{QED}(k^2).$$
(12)
After calculating convolutions over colour and Dirac indices we obtain the contribution of diagram (b) in the form

$$M_{(b)}=i\underset{0}{\overset{1}{}}𝑑u\varphi _\pi (u)\frac{(t^at^a)_{ij}}{4N_c}\frac{2\overline{s}g^2(e_qeZ)f_\pi }{k^2}\frac{\overline{u}(q_1)\frac{\widehat{p}_2^{}}{\overline{s}}(\widehat{q}_1\widehat{k})\gamma _\mu \gamma _5\widehat{p}_1\gamma _\mu v(q_2)}{(q_1k)^2(p_1\overline{u}q_2)^2}F_{QED}(k^2).$$
(13)
The produced quark and antiquark which have the colours $`i,j`$ are described by the spinors $`\overline{u}(q_1),v(q_2)`$. $`g`$ is the strong coupling constant. The electric charge of the nucleus we denote as $`eZ`$, whereas $`e_q`$ is the electric charge of the quark.

The contributions of the other diagrams in Fig. 1 can be writen down in a similar way. One remark is now in order. As we stressed above the hard scattering amplitude can be calculated in leading twist accuracy by considering the pion splitting into on-shell quarks. Then, the sum of diagrams Fig.1(a) and Fig.1(b) as well as the sum of diagrams Fig.1(c) and Fig.1(d) are separately gauge invariant, i.e. both sums vanish when $`\widehat{p}_2^{^{}}`$ in the photon vertex is substituted by $`\widehat{k}`$. This property together with Eqs. (4,5) permits us to substitute in all quark-photon vertices

$$\widehat{p}_2^{^{}}\frac{\overline{s}}{M_{2j}^2t}\widehat{k}_{},$$
(14)
since the term proportional to $`p_1`$ in Eq.(4) gives in our kinematics a contribution suppressed by power of $`\overline{s}`$. The substitution (14), being the consequence of gauge invariance, implies that the coupling of the photon to the $`q\overline{q}`$ pair vanishes linearly when $`k_{}0`$.

The final result for the scattering amplitude corresponding to the sum of diagrams shown in Fig. 1 is given by the formula

$`M^{\pi ^++A2j+A}=i\delta _{ij}f_\pi {\displaystyle \frac{2^6}{3^2}}\pi ^2\alpha ZF_{QED}(k^2)\alpha _s(q_1){\displaystyle \frac{\overline{s}}{(k_{}^2+k_{min}^2)(q_1^2)^2}}`$
$`(e_u\overline{z}+e_dz)\overline{u}(q_1)\gamma ^5\left[z\widehat{k}_{}\widehat{q}_1+\overline{z}\widehat{q}_1\widehat{k}_{}\right]{\displaystyle \frac{\widehat{p}_2^{}}{\overline{s}}}v(q_2){\displaystyle \underset{0}{\overset{1}{}}}𝑑u{\displaystyle \frac{\varphi _\pi (u)}{u}}.`$ (15)
where $`\alpha _s=\frac{g^2}{4\pi }`$, and we used the symmetry property $`\varphi _\pi (u)=\varphi _\pi (\overline{u})`$. The behaviour of the pion distribution amplitude at the end-points is known: $`\varphi _\pi (u)u`$ for $`u0`$, and $`\varphi _\pi (u)\overline{u}`$ for $`\overline{u}0`$ . Thus the integral over the momentum fraction fraction $`u`$ is well defined what confirms that factorization holds for the process we discuss. Eq. (S0.Ex1) is the main result of our paper.

Let us emphasize that the integral of $`\varphi _\pi (u)`$ over $`u`$ in Eq. (S0.Ex1) generates only an overall factor. Therefore the dependence of the amplitude $`M`$ on $`z`$ is universal, i.e. it doesn’t depend on the shape of the pion distribution amplitude. It is interesting also to note that the amplitude (S0.Ex1) vanishes for $`z=e_u/(e_ue_d)=2/3`$. Due to the opposite signs of the electric charges of the pion constituents, $`e_u=2e/3,e_d=e/3`$, the contribution of the diagrams (a) and (b) in Fig. 1 cancels the one of diagrams (c) and (d).

The expression in the square bracket on the r.h.s. of Eq.(S0.Ex1) can be put in the form

$$\overline{u}(q_1)\gamma ^5[\mathrm{}]\frac{\widehat{p}_2^{}}{\overline{s}}v(q_2)=\sqrt{u\overline{u}}\chi _q[(k_{}q_1)+i([k_{}\times q_1]\stackrel{}{\sigma })(z\overline{z})]\chi _{\overline{q}}$$
(16)
where $`\chi _q,\chi _{\overline{q}}`$ are two dimensional spinors of quark and antiquark. Since the incoming pion fluctuates into a $`q\overline{q}`$ states with total helicity zero we have in principle to consider two amplitudes having different spin configurations of quark and antiquark as denoted by $`()`$ and $`()`$. However, according to Eq.(16) $`M()=M()^{}`$, so effectively we deal here with only one independent helicity amplitude.

The Colomb contribution to the process (1) was considered in Ref. . Our result (S0.Ex1) differs from the formula (78) of . According to the factorization formula (9) we express our result in terms of the pion distribution amplitude $`\varphi _\pi (u)`$. We predict the universal $`z`$ dependence of the scattering amplitude $`M`$, independent of the shape of $`\varphi _\pi (u)`$, moreover $`M`$ vanishes for $`z=2/3`$. Whereas according to Eq. (78) of  the amplitude is proportional to the lowest pion Fock state wave function, $`M\mathrm{\Psi }_\pi (z,q_1)`$ (in our notation). This disagreement is, in our opinion, related to the treatment of the hard gluon exchange in : it is not included in the hard scattering block as is usually done within the QCD factorization approach but it is attributed to the high transverse momentum tail of $`\mathrm{\Psi }_\pi (z,q_1)`$.

4.  Next, we estimate the cross section. For simplicity we use the asymptotic pion distribution amplitude $`\varphi _\pi =6u\overline{u}`$. Using the standard rules and performing the sum over the colours and over the helicities of produced quarks we obtain

$`(q_1^2)^4{\displaystyle \frac{d\sigma ^{\pi ^++A2j+A}}{d^2k_{}d^2q_1dz}}=f_\pi ^2{\displaystyle \frac{2^5}{3\pi }}\alpha ^2\alpha _s^2(q_1)Z^2(e_u\overline{z}+e_dz)^2`$ (17)
$`{\displaystyle \frac{1}{(k_{}^2+k_{min}^2)^2}}F_{QED}^2(k_{}^2)\left[(z\overline{z})^2\left(q_1^2k_{}^2(q_1k_{})^2\right)+(q_1k_{})^2\right].`$
For the electromagnetic formfactor of a nucleus (Platinum, $`Z=78`$) we use a naive dipole approximation

$$F_{QED}(k^2)=\frac{\mathrm{\Lambda }^2}{\mathrm{\Lambda }^2+k^2},\mathrm{\Lambda }=60\text{MeV}.$$
(18)
In Fig. 2 we present the differential cross section (17) integrated over $`k_{}`$ for Platinum and the pion energy $`E=500\text{GeV}`$. Due to the different charges of $`u`$ and $`d`$ quarks the cross section has an asymmetric shape under the intechange $`z\overline{z}`$. As already mentioned it is zero for $`z=2/3`$. The cross section vanishes also for $`z\mathrm{\hspace{0.17em}0}`$ and $`z\mathrm{\hspace{0.17em}1}`$, this suppression is caused by the electromagnetic formfactor. According to Eq. (8), in the vicinity of the end-points $`z(\overline{z})q_1^2/(2E\mathrm{\Lambda })`$, the minimal momentum transfer becomes much larger than the scale of the electromagnetic formfactor $`\mathrm{\Lambda }`$.

It is quite difficult to distinguish jets originating from $`u`$ and $`d`$ quarks (by identifying the leading hadrons), therefore we show in Fig. 3 the symmetrized cross section, $`d\sigma ^{sym}(z)=1/2(d\sigma (z)+d\sigma (\overline{z}))`$. The dip structure in the plot for the symmetrized cross section is due to above mentioned zero of the amplitude (S0.Ex1). One could think that the shapes of the plots in Fig. 3 resemble the form of the Chernyak-Zhitnitsky distribution amplitude , i.e. that $`d\sigma ^{sym}(\varphi _\pi ^{(CZ)}(z))^2`$. This similarity is nevertheless accidental since as we explained above the dependence of the amplitude $`M`$ on $`z`$ is universal and the shape of curves in Fig. 3 has other reasons.

The Coulomb effect studied in this paper is complementary to the contribution from Pomeron exchange. Although the QCD part of the cross section is large it scales like $`1/q_1^8`$, i.e. it decreases faster at large $`q_1`$ than QED contribution which behaves like $`1/q_1^6`$. The relative magnitude of the Coulomb part should also be enhanced in the region of small $`z`$ or $`\overline{z}`$. The QCD contribution to the cross section is proportional to $`\varphi _\pi ^2(z)F_{QED}^2(t)`$ which vanishes at $`z(\overline{z})0`$ as $`z^2(\overline{z}^2)`$. This behaviour should be compared with the QED contribution given by Eq. (17) which behaves, for $`z(\overline{z})0`$, as $`\text{const}F_{QED}^2(t)`$.

The comparison of our results on the Coulomb diffractive pion dissociation into two jets with theoretical predictions or with experimental data is difficult. On the one hand side the theoretical results existing in the literature  for the QCD contribution are controvercial. In a fortcomming paper  we try to solve this issue within the QCD factorization scheme. On the other hand, in the kinematical region covered by E791 experiment the Coulomb effect is difficult to observe due to the large Pomeron background. Also, one cannot compare our predictions for Coulomb contribution directly with E791 data since their absolute normalization is still not known.

Acknowledgments

We are grateful to Vladimir Braun, Leonid Frankfurt and Andreas Schäfer for discussions and critical remarks.

L.Sz. acknowledges support by DFG and the warm hospitality at Regensburg University. D. I. acknowledges the support of Alexander von Humboldt Stiftung.


text:
# The Navier-Stokes-alpha model of fluid turbulence

## 1 Introduction

The energy in a turbulent fluid flow cascades toward ever smaller scales until it reaches the dissipation scale, where it can be transformed into heat. This cascade – creating fluid motions at ever smaller scales – is a characteristic feature of turbulence. This feature is also the main difficulty in simulating turbulence numerically, because all numerical simulations will have finite resolution and will eventually be unable to keep up with the cascade all the way to the dissipation scale, especially for complex flows, e.g., near walls and interfaces.

The effects of subgrid-scale fluid motions occurring below the available resolution of numerical simulations must be modeled. One way of modeling these effects is simply to discard the energy that reaches such subgrid scales. This is clearly unacceptable, though, and many creative alternatives have been offered. A prominant example is the large eddy simulation (LES) approach, see, e.g., , . The LES approach is based on applying a spatial filter to the Navier-Stokes equations. The reduction of flow complexity and information content achieved in the LES approach depends on the characteristics of the filter that one uses, its type and width. In particular, the LES approach introduces a length scale into the description of fluid dynamics, namely, the width of the filter used. Note that the LES approach is conceptually different from the Reynolds Averaged Navier-Stokes (or, RANS) approach, which is based on statistical arguments and exact ensemble averages, rather than spatial and temporal filtering. After filtering, however, just as in the RANS approach, one faces the classic turbulence closure problem: How to model the filtered-out subgrid scales in terms of the remaining resolved fields? In practice, this problem is compounded by the requirement that the solution be simulated numerically, thereby introducing further approximations.

This paper begins by reviewing a modeling scheme — called here the Navier-Stokes-alpha model, or NS$`\alpha `$ model (also called the viscous Camassa-Holm equations in  \-  ) — that introduces an energy “penalty” inhibiting the creation of smaller and smaller excitations below a certain length scale (denoted alpha). This energy penalty results in a nonlinearly dispersive modification of the Navier-Stokes equations. The alpha-modification appears in the nonlinearity, it depends on length scale and we emphasize that it is dispersive, not dissipative. We shall re-derive the modified equations in Section 2 from the viewpoint of Kelvin’s circulation theorem. As we shall show in section 3, this modification causes the translational kinetic energy wavenumber spectrum of the NS$`\alpha `$ model to roll off rapidly below the length scale alpha as $`k^3`$ in three dimensions, instead of continuing to follow the slower Kolmogorov scaling law, $`k^{5/3}`$. This roll off shortens the inertial range of the NS$`\alpha `$ model and thus makes it more computable. This is the main new result of the paper, so we shall comment now on its main implications.

Since the energy spectrum rolls off faster, the wavenumber $`k=\kappa _\alpha `$ at which viscous dissipation takes over in the NS$`\alpha `$ model must be lower than for the original Navier-Stokes equations. Hence, for a given driving force and viscosity, the number of active degrees of freedom $`N_{dof}`$ for the NS$`\alpha `$ model must be smaller than for Navier-Stokes. In section 3 we give a heuristic estimate of the number $`N_{dof}`$ and compare it with the rigorous estimate derived in  for the fractal dimension $`D_{frac}`$ of the global attractor for the NS$`\alpha `$ model. Namely,

$$D_{frac}(N_{dof})^{3/2}\text{and}N_{dof}(L\kappa _\alpha )^3\frac{L}{\alpha }Re^{3/2},$$
(1.1)
where $`L`$ is the integral scale (or domain size), $`\kappa _\alpha `$ is the end of the NS$`\alpha `$ inertial range and $`Re=L^{4/3}ϵ_\alpha ^{1/3}/\nu `$ is the Reynolds number (with NS$`\alpha `$ energy dissipation rate $`ϵ_\alpha `$ and viscosity $`\nu `$). The corresponding number of degrees of freedom for a Navier-Stokes flow with the same parameters is

$$N_{dof}^{NS}(L/\mathrm{}_{Ko})^3Re^{9/4},$$
(1.2)
where $`\mathrm{}_{Ko}`$ denotes the Kolmogorov dissipation length scale. The implication of these estimates of degrees of freedom for numerical simulations that access a significant number of them using the NS$`\alpha `$ model would be an increase in computational speed relative to Navier-Stokes of

$$\left(\frac{N_{dof}^{NS}}{N_{dof}}\right)^{4/3}=\left(\frac{\alpha }{L}\right)^{4/3}Re.$$
(1.3)
Thus, if $`\alpha `$ tends to a constant value, say $`L/100`$, when the Reynolds number increases – as found in  \-  by comparing steady NS$`\alpha `$ solutions with experimental data for turbulent flows in pipes and channels – then one could expect to obtain a substantial increase in computability by using the NS$`\alpha `$ model at high Reynolds numbers. An early indication of the reliability of using these estimates to gain a relative increase in computational speed in direct numerical simulations (DNS) of homogeneous turbulence in a periodic domain is given in . There, a computational speed up was achieved by using the NS$`\alpha `$ model in DNS of turbulence in a periodic domain by a factor about equal to the fluctuation Reynolds number ($`250`$). In the case considered in , this factor happens to be about $`Re/100`$.

In section 4, we discuss the relation of the NS$`\alpha `$ model to similar equations which are derived in a different physical context, namely in the context of non-Newtonian fluids. The main difference between the NS$`\alpha `$ model and the non-Newtonian second-grade fluids which it resembles lies in the choice of dissipation. The NS$`\alpha `$ model uses the standard Navier-Stokes viscosity, while the second-grade fluid uses a weaker form of dissipation – namely, wavenumber independent damping of the fluid velocity. Discussions of the relative advantages of the two forms of dissipation (e.g., in terms of boundary data requirements and well-posedness) are beyond the scope of this review. However, we remark that the steady solutions of the NS$`\alpha `$ model with the standard Navier-Stokes viscosity were found in  \-  to agree with experimental mean velocity profile data for turbulent flows in pipes and channels. The corresponding steady solutions of the second-grade fluid with its weaker form of dissipation were not found to so agree with this experimental data. Finally, in section 5, we provide a brief guide to the recent literature for those who might be interested in the mathematical context in which the alpha model was originally derived, as well as in its potential applications in turbulence modeling. In regard to the latter, the NS$`\alpha `$ model was recently shown to be transformable to a generalized similarity model for large eddy simulation (LES) turbulence modeling .

## 2 Kelvin-filtered turbulence models

Although it was first derived from energetic considerations using the Euler-Poincaré variational framework in , , the Navier-Stokes-alpha model may be motivated by an equivalent argument based on Kelvin’s circulation theorem. The original Navier-Stokes (NS) equations are

$$\frac{𝐯}{t}+𝐯𝐯+p=\nu \mathrm{\Delta }𝐯+𝐟,\text{with}𝐯=0,$$
(2.1)
for a forcing $`𝐟`$ and constant kinematic viscosity $`\nu `$. These equations satisfy Kelvin’s circulation theorem,

$$\frac{d}{dt}_{\gamma (𝐯)}𝐯𝑑𝐱=_{\gamma (𝐯)}(\nu \mathrm{\Delta }𝐯+𝐟)𝑑𝐱,$$
(2.2)
for a fluid loop $`\gamma (𝐯)`$ that moves with velocity $`𝐯(𝐱,t)`$, the Eulerian fluid velocity.

### Kelvin-filtering the Navier-Stokes equations.

The equations for the NS$`\alpha `$ model emerge from a modification of the Kelvin circulation theorem (2.2) to integrate around a loop $`\gamma (𝐮)`$ that moves with a spatially filtered Eulerian fluid velocity given by $`𝐮=g𝐯`$, where $``$ denotes the convolution,

$$𝐮=g𝐯=g(𝐱𝐲)𝐯d^3y.$$
(2.3)
The “inverse” is denoted

$$𝐯=𝒪𝐮,$$
(2.4)
thereby defining an operator $`𝒪`$ whose Green’s function is the filter $`g`$ and which we shall assume is positive, symmetric, isotropic and time-independent. Under these assumptions the quantity (kinetic energy)

$$E=\frac{1}{2}𝐮𝐯d^3x=\frac{1}{2}𝐯g𝐯d^3x=\frac{1}{2}𝐮𝒪𝐮d^3x,$$
(2.5)
defines a norm.

We obtain a modification to the Navier-Stokes equations (2.1) by replacing in their Kelvin’s circulation theorem (2.2) the loop $`\gamma (𝐯)`$ with another loop $`\gamma (𝐮)`$ moving with the spatially filtered velocity, $`𝐮`$. Then we have,

$$\frac{d}{dt}_{\gamma (𝐮)}𝐯𝑑𝐱=_{\gamma (𝐮)}(\nu \mathrm{\Delta }𝐯+𝐟)𝑑𝐱.$$
(2.6)
After taking the time derivative inside the Kelvin loop integral moving with filtered velocity $`𝐮`$ and reconstructing the gradient of pressure, we find the Kelvin-filtered Navier-Stokes equation,

$$\frac{𝐯}{t}+𝐮𝐯+𝐮^\mathrm{T}𝐯+p=\nu \mathrm{\Delta }𝐯+𝐟,$$
(2.7)
with

$$𝐮=0,\text{and}𝐯=𝒪𝐮.$$
(2.8)
The velocity $`𝐮(𝐱,t)`$ is the spatially filtered Eulerian fluid velocity in equation (2.3). Note that continuity equation is now imposed as $`𝐮=(g𝐯)=0`$. The energy balance relation derived from the Navier-Stokes-alpha equations (2.7) is

$$\frac{d}{dt}\frac{1}{2}𝐮𝐯d^3x=𝐮𝐟d^3x\nu |𝒪^{1/2}𝐮|^2d^3x,$$
(2.9)
where for the moment we have dropped the boundary terms that appear upon integrating by parts. That is, for the moment we ignore the boundary effects and consider either the case of the whole space with solutions vanishing sufficiently rapidly at infinity, or the case of periodic boundary conditions.

### Similarities with previous work.

1. Except for the term $`(𝐮)𝐯`$, the Kelvin-filtered Navier-Stokes equation (2.7) is otherwise quite similar to Leray’s regularization of the Navier-Stokes equations proposed in 1934 . Extension of the Leray regularization to satisfy the Kelvin circulation theorem was cited as an outstanding problem in Gallavotti’s review . Looking at the equation of motion for the vorticity $`𝐪=\times 𝐯`$, reveals even more similarity with Leray’s regularization of the Navier-Stokes equations. (See the section about the vorticity below.)
2. At first glance, the Kelvin-filtered equations (2.7) in the absence of dissipation and forcing may seem reminiscent of a form of the Euler equations that was discussed in Kuz’min  and Oseledets  (see also Gama and Frisch ), namely,

$`{\displaystyle \frac{\gamma _i}{t}}+u^j{\displaystyle \frac{\gamma _i}{x^j}}`$ $`=`$ $`\gamma _j{\displaystyle \frac{u^j}{x^i}},`$
$`\gamma _i`$ $`=`$ $`u_i+{\displaystyle \frac{\varphi }{x^i}},`$
$`{\displaystyle \frac{u_i}{x^i}}`$ $`=`$ $`0.`$ (2.10)
In light of the third equation in this set, the second one is essentially the Hodge decomposition of the vector $`\gamma `$ with Euclidean components $`\gamma _i`$. Comparing these equations with Euler’s equations for the incompressible motion of an ideal fluid,

$$\frac{u_i}{t}+u^j\frac{u_i}{x^j}+p=0,\frac{u_i}{x^i}=0,$$
(2.11)
gives a relation between the pressure gradient and the “gauge function” $`\varphi `$ that is reminiscent of (but different from) Bernoulli’s law,

$$p=\left(\frac{1}{2}|𝐮|^2\frac{\varphi }{t}u^j\frac{\varphi }{x^j}\right).$$
(2.12)
This type of relationship also arises in Hamilton’s principle when one uses Clebsch variables , in which case $`\varphi `$ is a Lagrange multiplier that enforces the continuity equation. In contrast, equations (2.7) in the ideal unforced case give the Kelvin-filtered Euler model,

$`{\displaystyle \frac{v_i}{t}}+u^j{\displaystyle \frac{v_i}{x^j}}`$ $`=`$ $`v_j{\displaystyle \frac{u^j}{x^i}}{\displaystyle \frac{p}{x^i}},`$ (2.13)
$`v_i`$ $`=`$ $`𝒪u_i,`$
$`{\displaystyle \frac{u_i}{x^i}}`$ $`=`$ $`0.`$
The Kuz’min-Oseledets form of the Kelvin-filtered Euler model appears by replacing $`u_iv_i`$ only in the second of the three equations in (2).

### Conservation laws.

References to the interesting geometrical properties of the Euler$`\alpha `$ model equations – namely equations (2.13) when $`𝒪`$ is the Helmholtz operator – are cited in the last Section. At this point, we only comment that in the case of periodic boundary conditions, or the case of the whole space with solutions vanishing sufficiently rapidly at infinity, these nondissipative equations preserve the following two quadratic invariants, the kinetic energy,

$$E=\frac{1}{2}𝐮𝐯d^3x,$$
(2.14)
and the $`𝐯`$ helicity,

$$\mathrm{\Lambda }=𝐯\text{curl}𝐯d^3x.$$
(2.15)
The kinetic energy, $`E_\alpha `$, and the $`𝐯`$-helicity, $`\mathrm{\Lambda }`$, defined above, are also preserved in bounded domains provided appropriate boundary conditions are imposed. Boundary conditions sufficient for energy conservation when $`𝒪`$ is the Helmholtz operator $`1\alpha ^2\mathrm{\Delta }`$ are

$$\widehat{𝐧}\times \left(\widehat{𝐧}(𝐮+𝐮^T)\right)=0,$$
(2.16)
where the vector $`\widehat{𝐧}`$ is normal to the boundary and superscript $`()^T`$ denotes matrix transpose. Boundary conditions sufficient for $`𝐯`$-helicity conservation in general are $`\lambda \widehat{𝐧}𝐮=0`$ and $`p\widehat{𝐧}\text{curl}𝐯=0`$, as seen from the helicity equation,

$$\frac{\lambda }{t}=\text{div}(\lambda 𝐮+p\text{curl}𝐯),\text{where}\lambda 𝐯\text{curl}𝐯.$$
(2.17)
This equation is obtained by using only the motion equation in (2.13) and its curl. Therefore, conservation of helicity holds with these boundary conditions simply because of the form of the Kelvin-filtered motion equation in (2.13), independently of the relation between $`𝐯`$ and $`𝐮`$, and regardless of whether $`𝐮`$ is incompressible.

### Vortex transport and stretching.

Let $`𝐪=\times 𝐯`$ be the vorticity of the unfiltered Eulerian velocity. The curl of the Kelvin-filtered Navier-Stokes equation (2.7) gives the vortex transport and stretching equation,

$$\frac{𝐪}{t}+𝐮𝐪𝐪𝐮=\nu \mathrm{\Delta }𝐪+\times 𝐟.$$
(2.18)
We note that the coefficient $`𝐮`$ in the vortex stretching term is the gradient of the spatially filtered Eulerian velocity $`𝐮`$. Thus, Kelvin-filtering tempers the vortex stretching in the modified Navier-Stokes equations (2.7), while preserving the original form of the vortex dynamics. This tempered vorticity stretching is also reminiscent of Leray’s approach  for regularizing the Navier-Stokes equations.

Indeed, when $`𝒪`$ is the Helmholtz operator we take full advantage of this regularized vorticity stretching effect to prove in  the global existence and uniqueness of strong solutions for the three dimensional NS$`\alpha `$ model. In fact, when the forcing is in a certain Gevrey class of regularity (real analytic), then the NS$`\alpha `$ solutions are also Gevrey regular, so that the Fourier coefficients of the solution decay exponentially fast, with respect to the wavenumbers, at a rate that increases as $`\alpha ^2`$ increases. See ,  for the corresponding analysis of the Gevrey properties of the Navier-Stokes equations. This result suggests that the filtering, or smoothing of the velocity $`𝐮`$ due to the presence of alpha in the NS$`\alpha `$ model could enhance the decay of its wavenumber spectrum and produce an earlier observation of this exponentially falling tail when the forcing is sufficiently smooth, even when the viscosity is small.

### Specializing to the Navier-Stokes-alpha model.

The special case of the NS$`\alpha `$ model emerges from the Kelvin-filtered Navier-Stokes equation (2.7) when we choose the operator $`𝒪`$ to be the Helmholtz operator, thereby introducing a constant $`\alpha `$ that has dimensions of length,

$$𝒪=1\alpha ^2\mathrm{\Delta },\text{with}\alpha =const.$$
(2.19)
In this case, the filtered and unfiltered fluid velocities in equation (2.4) are related by

$$𝐯=(1\alpha ^2\mathrm{\Delta })𝐮.$$
(2.20)
The Navier-Stokes-alpha model is given by the Kelvin-filtered Navier-Stokes equation (2.7) with definition (2.20). The original derivation of the ideal Euler-alpha model (the $`\nu =0`$ case of the NS$`\alpha `$ model) obtained by using the Euler-Poincaré approach is given in , . The physical interpretations of $`𝐮`$ and $`𝐯`$ as the Eulerian and Lagrangian mean velocities are given in .

The corresponding kinetic energy norm (2.5) for the NS$`\alpha `$ model is given by

$$E_\alpha =\left[\frac{1}{2}|𝐮|^2+\frac{\alpha ^2}{2}|𝐮|^2\right]d^3x,$$
(2.21)
and we repeat that $`𝐮`$ is the spatially filtered Eulerian fluid velocity. This kinetic energy is the sum of a translational kinetic energy based on the filtered velocity $`𝐮`$, and a gradient-velocity kinetic energy, multiplied by $`\alpha ^2`$. Thus, by showing the global boundedness, in time, of the kinetic energy (2.21) one concludes that the coefficient $`𝐮`$ in the vortex stretching relation (2.18) is bounded in $`L^2`$ for the NS$`\alpha `$ model. The second term in the kinetic energy (2.21) imposes an energy penalty for creating small scales. The spatial integral of $`|𝐮|^2`$ in the second term has the same dimensions as filtered enstrophy (the spatial integral of $`|\times 𝐮|^2`$, the squared filtered vorticity). For a domain with boundary the spatial integral of $`|𝐮|^2`$ in the second term in the kinetic energy norm (2.21) is replaced by the integral of $`\mathrm{trace}(4𝐃𝐃)`$, where $`𝐃`$ is the strain rate tensor, $`𝐃=(1/2)(𝐮+𝐮^\mathrm{T})`$ in Euclidean coordinates . For the case we treat here, in Euclidean coordinates and in the absence of boundaries, all of these norms are equivalent for an incompressible flow.

## 3 Spectral scaling for the NS$`\alpha `$ model

### A preliminary scaling argument.

Scaling ideas originally due to Kraichnan  for the case of two-dimensional turbulence may be applied to estimate the effect of the second term in the NS$`\alpha `$ model energy (2.21) on the energy spectrum in the present case. For sufficiently small wavenumbers ($`k\alpha 1`$, but $`kL1`$) the first term in the energy $`E_\alpha `$ in (2.21) dominates and the second term may be neglected. In this wavenumber region, the standard Kolmogorov scaling argument for turbulence gives a $`k^{5/3}`$ spectrum by the usual dimensional argument for the inertial range,

$$k=[L]^1,ϵ_\alpha =[L]^2[T]^3,E_\alpha (k)=[L]^3[T]^2=ϵ_\alpha ^ak^b,\text{if}k\alpha 1,$$
(3.1)
so that $`a=2/3`$ and $`b=5/3`$. Conversely, for sufficiently large wavenumbers ($`k\alpha 1`$) and the second term dominates in $`E_\alpha `$. A preliminary scaling argument would then indicate a $`k^3`$ spectrum for $`E_\alpha `$ in the large wavenumber region, according to

$$k=[L]^1,\eta _\alpha =[T]^3,E_\alpha (k)=[L]^3[T]^2=\eta _\alpha ^ak^b,\text{if}k\alpha 1,$$
(3.2)
so that $`a=2/3`$ and $`b=3`$ in this case. Here $`\eta _\alpha `$ is the rate of dissipation of the $`|𝐮|^2`$ part of the kinetic energy (with the same dimensions as enstrophy). Therefore, in the wavenumber region near $`k\alpha 1`$ the kinetic energy spectrum may be expected to have a break in slope and roll off from $`k^{5/3}`$ to $`k^3`$ scaling .

The expectation from this preliminary scaling argument seems to be confirmed in direct numerical simulations of the NS$`\alpha `$ model, as shown in Figure 1 below taken from . Figure 1 shows the energy spectra resulting from three DNS of the NS$`\alpha `$ model with mesh sizes of $`256^3`$ for three cases: with $`\alpha /L=0`$ (the Navier-Stokes equations), $`1/32`$ and $`1/8`$ for the same viscosity $`\nu =0.001`$. The corresponding Taylor microscale Reynolds numbers $`R_\lambda `$ are reported as $`147,182`$ and $`279`$, respectively. The higher $`\alpha `$ (higher Reynolds number) flows are found numerically to have more compact energy spectra.

### A more refined argument.

A more refined argument for the spectral scaling of the NS$`\alpha `$ model shall now be given that depends explicitly upon properties of the nonlinearity. Following , let $`u_\kappa `$ denote the component of $`𝐮`$ formed by the Fourier modes of fluid velocity $`𝐮`$ with wavenumbers in $`[\kappa ,2\kappa )`$ and similarly let $`v_\kappa `$ denote the corresponding component of $`𝐯`$, where $`\kappa `$ is well beyond the active wavenumbers of the driving force. The NS$`\alpha `$ energy balance for this component is then

$$\frac{1}{2}\frac{d}{dt}(u_\kappa ,v_\kappa )+\nu (\mathrm{\Delta }u_\kappa ,v_\kappa )=T_\kappa T_{2\kappa },$$
(3.3)
where $`T_\kappa `$ and $`T_{2\kappa }`$ denote, respectively, the energy transfer rate at $`\kappa `$ from the low wavenumber modes $`u_<`$ to the high wavenumber modes $`u_\kappa +u_>`$, and at $`2\kappa `$ from $`u_<+u_{2\kappa }`$ to $`u_>`$, where $`u_<`$ and $`u_>`$ are defined as

$$u_<\underset{j<\kappa }{}u_j,\text{and}u_>\underset{j2\kappa }{}u_j.$$
(3.4)
In particular,

$$T_\kappa =(\stackrel{~}{B}(u_<,v_<),u_\kappa )+(\stackrel{~}{B}(u_\kappa +u_>,v_\kappa +v_>),u_<),$$
(3.5)
where the bilinear operator $`\stackrel{~}{B}(u,v)`$ is given by

$$\stackrel{~}{B}(u,v)=P_\sigma \left(u\times (\times v)\right),$$
(3.6)
and $`P_\sigma `$ is the $`L^2`$ orthogonal projection (Leray projection), see, e.g.,  for more mathematical details. For equation (3.3) to hold, it is essential that $`\stackrel{~}{B}(u,v)`$ has the property,

$$(\stackrel{~}{B}(u,v),u)=0.$$
(3.7)
Note: property (3.7) does not hold for the enstrophy transfer rate of the NS$`\alpha `$ model, so the NS$`\alpha `$ model does not possess a Kraichnan type inertial range due to an enstrophy cascade. Time-averaging equation (3.3) gives

$$\nu (\mathrm{\Delta }u_\kappa ,v_\kappa )=T_\kappa T_{2\kappa }.$$
(3.8)
Consequently, upon introducing the energy spectrum $`E_\alpha (\kappa ),`$ this time-averaged energy transfer equation implies

$$\nu \kappa ^3E_\alpha (\kappa )\nu _\kappa ^{2\kappa }\kappa ^2E_\alpha (\kappa )𝑑\kappa T_\kappa T_{2\kappa }.$$
(3.9)
Hence, as long as the energy dissipation rate is small compared to the energy transfer rate, i.e., provided

$$\nu \kappa ^3E_\alpha (\kappa )T_\kappa ,$$
(3.10)
we have the inertial range condition,

$$T_\kappa T_{2\kappa }.$$
(3.11)
The total energy dissipation rate is estimated as

$$ϵ_\alpha =\frac{\nu }{L^3}_{[0,L]^3}u(\mathrm{\Delta }v)d^3x.$$
(3.12)
Kraichnan  posits the following mechanism for the turbulent cascade: In the inertial range the eddies $`u_\kappa `$ transfer their energy to the eddies $`u_{2\kappa }`$ in the time $`\tau _\kappa `$ it takes to travel their length $`1/\kappa `$. Their average velocity being

$$U_\kappa \frac{1}{L^3}_{[0,L]^3}u_\kappa u_\kappa d^3x^{1/2}=\left(_\kappa ^{2\kappa }\frac{E_\alpha (\kappa )d\kappa }{(1+\alpha ^2\kappa ^2)}\right)^{1/2}\left(\frac{\kappa E_\alpha (\kappa )}{(1+\alpha ^2\kappa ^2)}\right)^{1/2},$$
(3.13)
the eddy energy exchange, or turnover time $`\tau _\kappa `$ will be

$$\tau _\kappa 1/(\kappa U_\kappa ).$$
(3.14)
Consequently, the total energy dissipation rate is related to the $`\kappa `$ spectral energy density by using (3.13) as

$$ϵ_\alpha \tau _\kappa ^1_\kappa ^{2\kappa }E_\alpha (k)𝑑k\kappa U_\kappa _\kappa ^{2\kappa }E_\alpha (k)𝑑k\frac{\kappa ^{5/2}}{(1+\alpha ^2\kappa ^2)^{1/2}}E_\alpha (\kappa )^{3/2},$$
(3.15)
which yields the following spectral scaling law for the NS$`\alpha `$ inertial range,

$$E_\alpha (\kappa )\frac{ϵ_\alpha ^{2/3}(1+\alpha ^2\kappa ^2)^{1/3}}{\kappa ^{5/3}}.$$
(3.16)
Thus, the total energy present in the NS$`\alpha `$ inertial range is actually enhanced by the presence of alpha. Notice, however, that in terms of the filtered velocity $`𝐮`$ alone, the translational kinetic energy spectrum is given by

$$\frac{E_\alpha (\kappa )}{(1+\alpha ^2\kappa ^2)}\frac{ϵ_\alpha ^{2/3}}{\kappa ^{5/3}}\frac{1}{(1+\alpha ^2\kappa ^2)^{2/3}}=\{\begin{array}{cc}\frac{ϵ_\alpha ^{2/3}}{\kappa ^{5/3}}\hfill & \text{if}\kappa \alpha 1,\hfill \\ & \\ \frac{ϵ_\alpha ^{2/3}}{\alpha ^{4/3}\kappa ^3}\hfill & \text{if}\kappa \alpha 1.\hfill \end{array}$$
(3.17)
Hence, the anticipated $`\kappa ^3`$ behavior appears in the spectrum for the $`\alpha `$filtered velocity as a consequence of equation (3.17) arising from Kraichnan’s argument  associating energy transfer rates and eddy turnover times for the filtered velocity. According to this argument, the presence of alpha reduces the energy associated with the higher wave numbers in the $`L^2`$ norm of the filtered velocity $`𝐮`$. This conclusion from (3.17) agrees with the trends shown in Figure 1 obtained from high resolution DNS studies of the NS$`\alpha `$ equation in . However, these numerical studies do not have sufficient dynamic range to confirm this conclusion entirely. Thus, uncertainty remains in claiming a numerical confirmation because the quantity $`E_\alpha (\kappa )/(1+\alpha ^2\kappa ^2)`$ in equation (3.17) is unaffected by the alpha-modification for $`\kappa \alpha 1`$ and may already be out of the inertial range and into the dissipation range for $`\kappa \alpha 1`$. Studies in progress using DNS of the limit $`\alpha \mathrm{}`$ NS$`\alpha `$ equation hope to clarify this point . Next we shall discuss the extent of the inertial range for the NS$`\alpha `$ model.

Equation (3.17) holds only in the inertial range, that is, provided, cf. (3.10),

$$\nu \kappa ^3E_\alpha (\kappa )ϵ_\alpha .$$
(3.18)
This inertial-range inequality may be expressed equivalently as $`\kappa \kappa _\alpha `$ where $`\kappa _\alpha `$ (the end of the NS$`\alpha `$ inertial range) is given in terms of the NS$`\alpha `$ Kolmogorov wavenumber $`\kappa _{\alpha ,Ko}`$ by using (3.15) and (3.17) to find

$$\kappa _\alpha ^4(1+\alpha ^2\kappa _\alpha ^2)\frac{ϵ_\alpha }{\nu ^3}=\kappa _{\alpha ,Ko}^4.$$
(3.19)
Thus, as expected from the numerical simulations of , the inertial range is shortened to $`\kappa <\kappa _\alpha `$ for the NS$`\alpha `$ model by its nonlinear dispersive filtering with lengthscale $`\alpha `$. For sufficiently large NS$`\alpha `$ Kolmogorov wavenumber $`\kappa _{\alpha ,Ko}`$ and with $`\alpha `$ fixed, the wavenumber $`\kappa _\alpha `$ at the end of the NS$`\alpha `$ inertial range is determined from (3.19) to be

$$\kappa _\alpha \left(\frac{1}{\alpha }\right)^{1/3}\kappa _{\alpha ,Ko}^{2/3}.$$
(3.20)
This is a relationship among the three progressively larger wavenumbers,

$$1/\alpha <\kappa _\alpha <\kappa _{\alpha ,Ko}.$$

Shortening the inertial range for the NS$`\alpha `$ model to $`\kappa <\kappa _\alpha `$ rather than $`\kappa <\kappa _{\alpha ,Ko}`$ implies fewer active degrees of freedom in the solution.

### Counting degrees of freedom.

If one expects turbulence to be “extensive” in the thermodynamic sense, then one may expect that the number of “active degrees of freedom” $`N_{dof}`$ for alpha-model turbulence should scale as

$$N_{dof}(L\kappa _\alpha )^3(L/\alpha )(L\kappa _{\alpha ,Ko})^2\frac{L}{\alpha }Re^{3/2},$$
(3.21)
where $`L`$ is the integral scale (or domain size), $`\kappa _\alpha `$ is the end of the NS$`\alpha `$ inertial range and $`Re=L^{4/3}ϵ_\alpha ^{1/3}/\nu `$ is the Reynolds number (with dissipation rate $`ϵ_\alpha `$ and viscosity $`\nu `$). The corresponding number of degrees of freedom for Navier-Stokes with the same parameters is

$$N_{dof}^{NS}(L\kappa _{\alpha ,Ko})^3Re^{9/4},$$
(3.22)
and one sees a possible trade-off in the relative Reynolds number scaling of the two models. Should these estimates of the number of degrees of freedom needed for numerical simulations using the NS$`\alpha `$ model relative to Navier-Stokes not be overly optimistic, the implication would be one factor of $`(N_{dof}^{NS}/N_{dof})^{1/3}`$ in relative increased computational speed gained by the NS$`\alpha `$ model for each spatial dimension and yet another factor for the accompanying lessened CFL time step restriction. Altogether, this would be a gain in speed of

$$\left(\frac{N_{dof}^{NS}}{N_{dof}}\right)^{4/3}=\left(\frac{\alpha }{L}\right)^{4/3}Re.$$
(3.23)
Since $`\alpha /L1`$ and $`Re1`$, the two factors in the last expression compete, but the Reynolds number should eventually win out, because $`Re`$ can keep increasing while the number $`\alpha /L`$ is expected to tend to a constant value, say $`\alpha /L=1/100`$, at high (but experimentally attainable) Reynolds numbers, at least for simple flow geometries. Empirical indications for this tendency were found in  \-  by comparing steady NS$`\alpha `$ solutions with experimental mean-velocity-profile data for turbulent flows in pipes and channels.

Thus, according to this scaling argument, a factor of $`10^4`$ in increased computational speed for resolved scales greater than $`\alpha `$ could occur by using the NS$`\alpha `$ model at the Reynolds number for which the ratio $`\kappa _{\alpha ,Ko}/\kappa _\alpha =10`$. An early indication of the feasibility of obtaining such factors in increased computational speed was realized in the direct numerical simulations of homogeneous turbulence reported in , in which $`\kappa _{\alpha ,Ko}/\kappa _\alpha 4`$ and the full factor of $`4^4=256`$ in computational speed was obtained using spectral methods in a periodic domain at little or no cost of accuracy in the statistics of the resolved scales.

### Related mathematical results.

The paper  shows that strong solutions of the NS$`\alpha `$ model exist globally, they are unique, and they lie on a global attractor whose fractal (Lyapunov) dimension is bounded above by

$$𝒟_{frac}N_{dof}^{3/2},$$
(3.24)
with $`N_{dof}`$ defined as in equation (3.21). This rigorous mathematical bound exceeds the expected value obtained above by heuristically counting degrees of freedom assuming “extensive” turbulence. Thus, it would imply a smaller increase in computational speed than that of (3.23). However, this rigorous bound may have room for improvement.

### Oboukov cascade rate.

We have established in equation (3.14) that the eddy turnover rate based on the translational velocity for the NS$`\alpha `$ model is

$$\tau _\kappa ^1\kappa U_\kappa \sqrt{\frac{\kappa ^3E_\alpha (k)}{1+\alpha ^2\kappa ^2}}\left(\frac{ϵ_\alpha \kappa ^2}{1+\alpha ^2\kappa ^2}\right)^{1/3}.$$
(3.25)
Therefore,

* for low wavenumbers ($`\kappa \alpha 1`$), information (or error) propagates between scales $`\kappa `$ and $`2\kappa `$ at a rate proportional to $`\kappa ^{2/3}`$, in agreement with the classical accelerated cascade of Oboukov  as cited, e.g., in ; while
* for high wavenumbers ($`\kappa \alpha 1`$), information propagates at a constant rate, independently of wavenumber for the NS$`\alpha `$ model (as occurs also in the case of 2D turbulence discussed by Leith and Kraichnan ).

Without an accelerated cascade at high wavenumber, the NS$`\alpha `$ model should tend to be more predictable than the original Navier-Stokes model. Of course, it makes sense that an averaged, or filtered description of fluid flow would propagate high wavenumber information and errors at a slower rate than an instantaneous, unfiltered description does (e.g., think of climate vs weather). The $`\kappa ^3`$ spectrum of the translational kinetic energy in the alpha models for $`\alpha \kappa 1`$ is consistent with this interpretation of reduced error propagation rate as measured in the $`L^2`$ norm of the filtered velocity.

## 4 Rheology of NS$`\alpha `$ turbulence: 2nd-grade fluids

We rewrite the NS-alpha equations (2.7) with $`𝐯=(1\alpha ^2\mathrm{\Delta })𝐮`$ and $`\alpha `$ constant in their equivalent constitutive form,

$$\frac{d𝐮}{dt}=𝐓,\text{where}𝐓=p𝐈+2\nu (1\alpha ^2\mathrm{\Delta })𝐃+2\alpha ^2𝐃^{^{^{}}},$$
(4.1)
with $`𝐮=0`$, strain rate $`𝐃=(1/2)(𝐮+𝐮^\mathrm{T})`$, vorticity tensor $`\mathrm{\Omega }=(1/2)(𝐮𝐮^\mathrm{T})`$, and co-rotational (Jaumann) derivative given by $`𝐃^{^{^{}}}=d𝐃/dt+𝐃\mathrm{\Omega }\mathrm{\Omega }𝐃`$.

In equation (4.1), one recognizes the constitutive relation for the NS$`\alpha `$ model as a variant of the rate-dependent incompressible homogeneous fluid of second grade , in which the dissipation, however, is modified by composition with the Helmholtz operator $`(1\alpha ^2\mathrm{\Delta })`$. Thus, the NS$`\alpha `$ model has nearly the same stress tensor as the second-grade fluid. However, it is not quite the same. The stress tensor for the NS$`\alpha `$ model uses Navier-Stokes viscous dissipation, instead of the weaker form of dissipation used for second-grade fluids that is independent of wavenumber. Note: despite first appearances, there is no hyperviscosity in the NS$`\alpha `$ model, only the standard Navier-Stokes viscosity.

Equations for the second grade fluid were treated recently in the mathematical literature . Also, in  local well-posedness results for some initial value problems for differential type fluids arrising in non-Newtonian fluids (second-grade and third-grade fluid) were obtained by transferring the problem from the Eulerian to the Lagrangian setting, thereby extending the Arnold  and Ebin-Marsden program  to the case of second-grade fluids. For the case of second-grade fluids, for $`k\alpha >1`$ the wavenumber spectrum should still roll over to $`k^3`$, provided the Kolmogorov/Oboukov argument still holds for the weaker dissipation.

The association of turbulence closure models with non-Newtonian fluids is natural. There is a tradition at least since Rivlin  of modeling turbulence by using continuum mechanics principles such as objectivity and material frame indifference (see also ). For example, this sort of approach is taken in deriving Reynolds stress algebraic equation models . Rate-dependent closure models of mean turbulence have also been obtained by the two-scale DIA approach  and by the renormalization group method .

Despite its similarity to the rheology of second-grade fluids, the alpha parameter in the NS$`\alpha `$ model is actually not a material parameter. Rather, it is a flow regime parameter. Since the NS$`\alpha `$ model describes mean quantities, it was proposed as a turbulence closure model and this ansatz was tested by comparing its steady solutions (using the standard Navier-Stokes viscosity) to mean velocity measurements in turbulent channel and pipe flows in -. These experimental tests show that alpha depends slightly on Reynolds number and varies slightly with distance from the wall, for the low to moderate Reynolds numbers available in channel flow. However, at the high to very high Reynolds numbers available in pipe flows, alpha becomes independent of Reynolds number and takes a small value – about one percent of the pipe diameter.

## 5 Guide to recent NS$`\alpha `$ model literature

In - the authors introduce random fluctuations into the description of the fluid parcel trajectories in the Lagrangian in Hamilton’s principle for ideal incompressible fluid dynamics. They then take its statistical average and use the Euler-Poincaré theory of , ,  to derive Eulerian closure equations for the corresponding averaged ideal fluid motions. The Euler-Poincaré equation that is used in deriving the ideal dynamics of the NS$`\alpha `$ model is equivalent in the Eulerian picture to the corresponding Euler-Lagrange equation for fluid parcel trajectories for Lagrangians that are invariant under the right-action of the diffeomorphism group. See , , , and references therein for more discussions of Euler-Poincaré equations. The Euler-Poincaré theory is applied for modeling fluctuation effects on 3D Lagrangian mean and Eulerian mean fluid motion in .

The NS$`\alpha `$ model equations (also sometimes called in the literature the viscous Camassa-Holm equations) are proposed in - as a turbulence closure approximation for the Navier-Stokes equations. The analytic form of the velocity profiles based on the steady NS$`\alpha `$ equations away from the viscous sublayer, but covering at least 95% of the channel, depends on two free parameters: the flux Reynolds number $`R`$, and the wall-stress Reynolds number $`R_0`$. (Due to measurement limitations, most experimental data are contained in this region.) The authors further reduce the parameter dependence to a single free parameter by assuming a certain drag law for the wall friction $`DR_0^2/R^2`$. For most of the channel the steady NS$`\alpha `$ solution is shown to be compatible with empirical and numerical velocity profiles in this subregion. The NS$`\alpha `$ steady velocity profiles agree well outside the viscous sublayer with data obtained from mean velocity measurements and simulations of turbulent channel and pipe flow over a wide range of Reynolds numbers, -.

### DNS and comparisons with LES.

In  direct numerical simulations (DNS) of the NS$`\alpha `$ equations are compared with Navier-Stokes DNS and interpreted as behaving like an LES (Large Eddy Simulation) model. That is, the NS$`\alpha `$ model is shown to produce an accurate dynamical description of the large scale features of turbulence driven at large scales, even at resolutions for which the fine scales are not resolved. The NS$`\alpha `$-model may at first appear not to be an LES model, since it has rate-dependence that is not admitted by LES models. However, the NS$`\alpha `$ model was recently shown to transform to a generalized LES similarity model in . This is a promising result for LES modeling, since the mathematical theorems available in  for existence, uniqueness and finite dimensional global attractor for the NS$`\alpha `$ model will now be transferable to the continuous formulations of this class of generalized LES similarity models. As a basis for numerical schemes, the NS$`\alpha `$ model also resembles vorticity methods, as observed in , , .

### Geodesic motion.

The completely integrable one-dimensional Camassa-Holm equation  is expressed on the real line as

$$\frac{v}{t}+u\frac{v}{x}+2v\frac{u}{x}=0,u(x,t)=\frac{1}{2}_{\mathrm{}}^{\mathrm{}}e^{|(xy|}v(y,t)𝑑y.$$
(5.1)
Thus, we have $`v=u^2u/x^2`$, cf. equations (2.7) with definition (2.20). This one-dimensional equation is formally the Euler-Poincaré equation for geodesic motion on the diffeomorphism group with respect to the metric given by the mean kinetic energy Lagrangian, which is right-invariant under the action of the diffeomorphism group. See ,  for detailed discussions, applications and references to Euler-Poincaré equations of this type for ideal fluids and plasmas. See ,  for the original derivation of the $`n`$-dimensional Camassa-Holm, or Euler$`\alpha `$ equation in Euclidean space. See ,  for discussions of its generalization to Riemannian manifolds, its existence and uniqueness on a finite time interval, and more about its relation to the theory of second grade fluids. Additional properties of the Euler$`\alpha `$ equations, such as smoothness of the geodesic spray (the Ebin-Marsden theorem) are also shown to hold in  and the limit of zero viscosity for the corresponding viscous equations is shown to be a regular limit, even in the presence of boundaries for either homogeneous (Dirichlet) boundary conditions, or for boundary conditions involving the second fundamental form of the boundary. Functional-analytic studies of the ideal Euler$`\alpha `$ model are made in , . The methods introduced and applied in , , while geometrical in nature, also address analytical issues and obtain results that were not previously available by traditional techniques for partial differential equations. For example, these methods produce local in time existence of $`C^{\mathrm{}}`$ viscosity independent solutions for the Euler$`\alpha `$ equations in $`n`$ dimensions (and in particular 3D) for a fluid container that can be an arbitrary Riemannian manifold with boundary .

### Mathematical estimates for strong solutions with dissipation.

Paper  provides the mathematical estimates that are needed to show that the solutions of the NS-$`\alpha `$ model exist globally, are unique and possess a global attractor with finite fractal dimension satisfying equation (3.24). The estimates in  do depend on the viscosity remaining positive.

## Acknowledgements

We are grateful to S. Chen, J. A. Domaradzki, L. Margolin, J. E. Marsden, E. Olson, T. Ratiu, S. Shkoller, J. Tribbia, B. Wingate and S. Wynne for many constructive comments and enlightening discussions. DDH is also indebted to G. Eyink, A. Leonard, K. Moffat and D. Thomson for insightful discussions during our time together (June 1999) at the Turbulence Programme of the Isaac Newton Institute for Mathematical Sciences at Cambridge. Research by CF was supported in part by the National Science Foundation grant DMS-9706903. Research by DDH was supported by the U.S. Department of Energy under contracts W-7405-ENG-36 and the Applied Mathematical Sciences Program KC-07-01-01. The work of EST was supported in part by the National Science Foundation grants DMS-9704632 and DMS-9706964.