text:
# SPH Compressible Turbulence

## 1 Introduction

Many SPH simulations of astrophysical flows have regions where the SPH particles pass through each other. In a continuum simulation this streaming is unphysical because it implies that the velocity field is locally multi-valued. The problem can be cured by introducing sufficient viscosity, but this has undesirable side effects, for example excessive decay of shear. For this reason the version of SPH called XSPH was devised (Monaghan 1989, Monaghan 1992). The idea is simple. Each particle is moved with an smoothed velocity which is constructed by averaging the velocities of its neighbours. In this way the counterstreaming is prevented or greatly reduced The procedure is non-dissipative and conserves linear and angular momentum. The smoothing has the further advantage of reducing local disorder. Its defect is that energy is not conserved.

The way to conserve energy, and at the same time greatly enlarge the application of SPH is provided by the alpha model of turbulence (for early references see Holm, 1999). It is a minimalist model of turbulence which aims to take into account the finite resolution of any numerical simulation or experimental measurement by using an average velocity to move the fluid, while at the same time preserving the invariants of the flow. The latter is achieved by deriving the equations of motion from a Lagrangian which automatically provides the stress terms which are guessed in other models of turbulence (a good modern discussion of these models is given by Pope (2000)). The alpha equations preserve the circulation of the fluid (when this is an invariant of the original equations) as well as the additive integrals of energy and momentum.

A feature of the continuum alpha model is that the presence of velocity gradient terms in the force greatly reduces the deposit of energy in the high wave number modes. In the SPH version the force on a particle contains non linear terms which depend on the velocities of neighbouring particles. We conjecture that these terms drive the energy back and forth from the low and high wave numbers. A related example is the Fermi-Pasta-Ulam problem (Fermi et al. 1955, Segr$`\stackrel{`}{e}`$ 1965) where non linear force terms lead not to statistical equilibrium (which would mean most energy at the high wave numbers), but to a flow of energy back and forth amongst a few of the modes with low wave number. A similar back and forth flow occurs in three dimensional turbulence (Piomelli et al. 1991, Jimenez 1994 , Woodward et al. 2001, Pope 2000) although this is often obscured by describing just the nett flux which is from low to high wave numbers. In fact, at any time, different parts of a turbulent system can have energy flow either up or down the wave number scale.

The advantage of the alpha model is that, if the non dissipative system already reduces the energy in the high wave number modes, then it is possible to introduce weak dissipation to complete the task of removing energy. Other turbulence models try to do this by using a standard form of Navier Stokes equations with a dissipation strong enough to remove the energy at high wave numbers.

In this paper we derive the SPH alpha equations from a particle Lagrangian using the continuum constant density equations as a guide. We show that the system conserves a discrete form of the Kelvin circulation theorem in addition to the additive integrals of energy and momentum. The equations are simple to implement and the extra computational cost is negligible. Because the turbulence model only requires minimal dissipation, it should be possible to reduce the normal SPH viscosity significantly in regions of the flow away from shocks. By combining this with the idea of varying the viscosity coefficient for each particle separately (Morris and Monaghan 1997) the viscosity could be reduced by a factor $`100`$ away from shocks.

The plan of this paper is as follows. We derive a standard set of equations for SPH from a Lagrangian with a resolution length which is consistent with the density. The momentum equation is found to take the same form as that obtained by Springel and Hernquist (2001) using Lagrange multipliers to constrain the smoothing length. However, the thermal energy equation differs from any of their thermal energy equations. We show that the equations of motion lead to a discrete Kelvin circulation theorem in addition to the additive integrals of energy and momentum. We then define the XSPH velocity averaging and show that it agrees with the continuum alpha model. Guided by the alpha model we propose a discrete SPH Lagrangian. The Lagrangian equations of motion agree in the continuum limit with the equations of the alpha model. The conservation of circulation takes exactly the same form as in the case when the velocities are not averaged. Finally we discuss the dissipation and a scaling argument which suggests that the energy spectrum will fall off more rapidly than the Kolmogorov spectrum at short length scales. The extensive numerical work required to confirm the conjectures and the scaling argument will be presented elsewhere.

## 2 The Comparison equations

### 2.1 The momentum equation

In this section we derive a standard set of SPH equations in the absence of velocity averaging to provide a comparison with the alpha model equations. This standard set of SPH equations will be obtained for a Lagrangian allowing fully for the variation of resolution length $`h`$ with density. The equations are similar to those of Springel and Hernquist (2001), and Bonet (1999, 2001) however the energy equation differs from those considered by Springel and Hernquist (2001). The formulation can be easily extended to the relativistic case using the analysis of Monaghan and Price (2001).

The Lagrangian for compressible, non dissipative flow is (Eckart 1960)

$$L=\rho \left(\frac{1}{2}v^2u(\rho ,s)\right)\mathrm{𝐝𝐫},$$
(2.1)
where $`u(\rho ,s)`$ is the thermal energy per unit mass which is a function of density $`\rho `$ and entropy $`s`$. The SPH form of (2.1) with self gravity included is

$$L=\underset{b}{}m_b\left(\frac{1}{2}v_b^2u(\rho _b,s_b)+\frac{G}{2}\underset{k}{}\frac{m_km_b}{(|𝐫_b𝐫_k|)}\right),$$
(2.2)
where

$$\frac{d𝐫_a}{dt}=𝐯_a.$$
(2.3)
In these equations $`m_b`$ is the mass, $`𝐯_b`$ the velocity, $`𝐫_b`$ the position, $`s_b`$ the entropy, and $`\rho _b`$ the density of SPH particle $`b`$. In practice the gravitational term is smoothed to remove the singularity when $`𝐫_b=𝐫_k`$.

In the SPH formulation the density of particle $`a`$ can be written

$$\rho _a=\underset{b}{}m_bW_{ab}(h_a),$$
(2.4)
where $`h_a`$ is assumed to be a specified function of $`\rho _a`$ which we denote by $`H(\rho _a)`$ or $`H_a`$. In many astrophysical calculations $`H_a(1/\rho _a^{1/d})`$ where the number of dimensions is $`d`$, but a more general function could be used. For example, to prevent arbitarily large $`h`$ when $`\rho `$ becomes very small we could choose

$$H_a=\frac{A}{1+B\rho _a^{1/d}},$$

where $`A`$ and $`B`$ are constants. Furthermore, while the usual practice is to estimate $`\rho _a`$ at a given time using the value of $`h_a`$ from a previous time, it would be possible and more consistent to calculate $`\rho _a`$ self consistently from (2.4) as suggested recently by Bonet (2001). Equation (2.4) is a function defining the scalar quantity $`\rho _a`$ and it can be solved by any standard root solving algorithm using either the value at the previous time step as a starting value, or by estimating the value using the rate of change of the density.

The equations of motion follow from varying the action keeping the entropy constant. From Lagrange’s equations for particle $`a`$

$$\frac{d}{dt}\left(\frac{L}{𝐯_a}\right)\frac{L}{𝐫_a}=0,$$
(2.5)
we find

$$\frac{d𝐯_a}{dt}=\underset{b}{}m_b\left(\frac{u}{\rho }\right)_s\frac{\rho _b}{𝐫_a}G\underset{b}{}\frac{m_b(𝐫_a𝐫_b)}{|𝐫_a𝐫_b|^3}.$$
(2.6)
From (2.4)

$$\mathrm{\Omega }_b\frac{\rho _b}{𝐫_a}=\underset{c}{}m_c_aW_{ac}(h_a)\delta _{ab}m_a_bW_{ab}(h_b),$$
(2.7)
where $`_a`$ denotes the gradient of $`W`$ taken with respect to the coordinates of particle $`a`$ keeping $`h`$ constant (we keep this convention throughout this paper), and

$$\mathrm{\Omega }_b=1H_b^{}\underset{c}{}m_c\frac{W_{bc}(h_b)}{h_b}.$$
(2.8)
where $`H^{}`$ denotes $`H_b/\rho _b`$. If the density variation is smooth then $`\mathrm{\Omega }=1+O(h^2)`$.

The first law of thermodynamics, when the entropy is constant, gives

$$\left(\frac{u}{\rho }\right)_s=\frac{P}{\rho ^2},$$
(2.9)
where $`P`$ is the pressure. The acceleration equation (2.6) with (2.7) and (2.9) then becomes

$$\frac{d𝐯_a}{dt}=\underset{b}{}m_b\left(\frac{P_a}{\mathrm{\Omega }_a\rho _a^2}_aW_{ab}(h_a)+\frac{P_b}{\mathrm{\Omega }_b\rho _b^2}_aW_{ab}(h_b)\right)G\underset{b}{}\frac{m_b(𝐫_a𝐫_b)}{|𝐫_a𝐫_b|^3}.$$
(2.10)
The equation of motion (2.10) is exactly the same as the equation of motion due to Springel and Hernquist (2001) who introduce constraints on $`h`$ with a Lagrange multiplier. If the fluid is in a container the container can be represented by boundary forces (Monaghan 1994, Monaghan and Kos 1999).

### 2.2 Conservation Laws

The symmetry of the Lagrangian leads immediately to the conservation laws. In particular the density (and therefore P when $`s`$ is constant) and the gravitational potential are independent of translations and rotations of the coordinate system. Accordingly, linear and angular momentum are conserved. If the system is viscous then the viscous term can be constructed so that the linear and angular momentum are still conserved (see $`\mathrm{\S }2.3)`$. For the constant entropy case the conservation of energy can be shown directly by noting that

$$\frac{d}{dt}\left(\underset{a}{}𝐯_a\frac{L}{𝐯_a}L\right)=0.$$
(2.11)
If the entropy remains constant during the motion the particle system is invariant to other transformations. Consider, for example a set of particles each with the same mass and entropy marking out a necklace. Imagine each particle in the loop being shifted to its neighbour’s position (in the same sense around the necklace) and given its neighbour’s velocity. Since the entropy is constant, nothing has changed, and the Lagrangian is therefore invariant to this transformation. A discussion of the resulting invariant quantity has been given by Monaghan and Price (2001). We repeat the key elements of the argument here because we intend to show (see $`\mathrm{\S }6`$) that despite moving the fluid with a smoothed velocity field the same invariant holds.

The change in $`L`$ can be approximated by

$$\delta L=\underset{j}{}\left(\frac{L}{𝐫_j}\delta 𝐫_j+\frac{L}{𝐯_j}\delta 𝐯_j\right),$$
(2.12)
where $`j`$ denotes the label of a particle on the loop. Because the system is invariant to the transformation we take $`\delta l=0.`$ The change in position and velocity are given by

$$\delta 𝐫_j=𝐫_{j+1}𝐫_j,$$
(2.13)
and

$$\delta 𝐯_j=𝐯_{j+1}𝐯_j.$$
(2.14)
Using Lagrange’s equations we can rewrite (2.12) in the form

$$\delta L=\frac{d}{dt}\left(\underset{j}{}\frac{L}{𝐯_j}(𝐫_{j+1}𝐫_j)\right)=0.$$
(2.15)
Recalling that the particle masses are assumed identical, we deduce that

$$\frac{d}{dt}\underset{j}{}𝐯_j(𝐫_{j+1}𝐫_j)=0,$$
(2.16)
so that

$$C=\underset{j}{}𝐯_j(𝐫_{j+1}𝐫_j),$$
(2.17)
is conserved to this approximation, for every loop. The conservation is only approximate because the change to the Lagrangian is discrete, and only approximated by the first order terms. However, if the particles are sufficiently close together (2.17) approximates the circulation theorem to arbitrary accuracy. A related argument was used by Feynman (1957) to establish from the invariance of the wave function that circulation should be quantised.

The system is also invariant to the particles shifting around the loop in the opposite sense. This gives an approximation to the circulation with the opposite sign to that above. If these two are combined (taking account of their signs so we subtract one from the other) we get

$$\frac{d}{dt}\underset{j}{}𝐯_j\frac{(𝐫_{j+1}𝐫_{j1})}{2}=0.$$
(2.18)
which is a better approximation to the circulation of the continuous fluid. The errors in the discrete approximation to the circulation theorem arise because of the higher order terms in the change to the Lagrangian. These errors become less as the spatial scale of variation of the velocity field becomes greater than the particle separation. If the velocity is constant then

$$C=\frac{1}{2}𝐯\underset{j}{}(𝐫_{j+1}𝐫_{j1})=0,$$
(2.19)
as expected. If the velocity is due to rigid rotation with angular velocity $`𝛀`$, then

$$C=\frac{1}{2}𝛀\underset{j}{}𝐫_j\times (𝐫_{j+1}𝐫_{j1}).$$
(2.20)
Taking the origin of the coordinate system within the loop for convenience, the area of the triangle made by $`𝐫_j`$ and $`𝐫_{j+1}`$ is $`\frac{1}{2}(𝐫_j\times 𝐫_{j+1})`$. The area is counted twice so that if the loop lies in a plane with normal $`𝐧`$

$$C=2𝛀𝐧A,$$
(2.21)
where $`A`$ is the area within the loop of particles. The vorticity is $`2𝛀`$ and (4.18) gives the same result as would be obtained by applying Stokes theorem to the circulation integral. Further examples are given by Monaghan and Price (2001).

Salmon (1988), following Bretherton (1970) also establishes the circulation law by appealing to the invariance to particle interchange. However, because his analysis is based on the continuum equations, it is more complicated than our derivation.

### 2.3 Comments on thermal energy

If the entropy is constant we can obtain the rate of change of thermal energy from the rate of change of density. The rate of change of density is obtained by differentiating (2.4) with respect to time. We find

$$\frac{d\rho _a}{dt}=\frac{1}{\mathrm{\Omega }}_a\underset{b}{}m_b𝐯_{ab}_aW_{ab}(h_a),$$
(2.22)
where $`𝐯_{ab}=(𝐯_a𝐯_b)`$. The rate of change of thermal energy when the entropy is constant is then given by

$$\frac{du_a}{dt}=\frac{P_a}{\rho _a^2}\frac{d\rho _a}{dt}=\frac{P_a}{\rho _a^2\mathrm{\Omega }_a}\underset{b}{}m_b𝐯_{ab}_aW_{ab}(h_a),$$
(2.23)
which disagrees with all of the thermal energy equations proposed by Springel and Hernquist (2001).

Our thermal energy equation can also be deduced by first forming the vector dot product of $`m_a𝐯_a`$ and (2.10), followed by summation over $`a`$. We get

$$\frac{dE_K}{dt}=\underset{a}{}\underset{b}{}m_bm_a𝐯_a\left(\frac{P_a}{\mathrm{\Omega }_a\rho _a^2}_aW_{ab}(h_a)+\frac{P_b}{\mathrm{\Omega }_b\rho _b^2}_aW_{ab}(h_b)\right)\frac{dE_G}{dt},$$
(2.24)
where $`E_K`$ is the kinetic energy and $`E_G`$ is the gravitational energy. If the labels $`a`$ and $`b`$ are interchanged in that part of the double summation involving $`P_b`$ we can write (2.24) as

$$\frac{dE_K}{dt}+\frac{dE_G}{dt}=\underset{a}{}m_a\frac{P_a}{\mathrm{\Omega }_a\rho _a^2}\underset{b}{}m_b𝐯_{ab}_aW_{ab}(h_a)$$
(2.25)
Conservation of total energy then requires that the right hand side of (2.14) is $`dE_{\mathrm{𝑡ℎ}}/dt`$. Thus

$$\frac{dE_{\mathrm{𝑡ℎ}}}{dt}=\underset{a}{}m_a\frac{du_a}{dt}=\underset{a}{}m_a\frac{P_a}{\mathrm{\Omega }_a\rho _a^2}\underset{b}{}m_b𝐯_{ab}_aW_{ab}(h_a)$$
(2.26)
from which (2.23) follows.

In SPH calculations viscous dissipation is included by adding

$$\mathrm{\Pi }_{ab}=\frac{\sigma _{ab}𝐯_{ab}𝐫_{ab}}{|𝐫_{ab}|},$$
(2.27)
to the pressure terms (Monaghan 1997). The quantity $`\sigma _{ab}`$ is a positive definite parameter, invariant to the interchange of $`a`$ and $`b`$. It is related to the actual viscosity coefficient or to an artificial viscosity coefficient dependent on the speed of sound (see Monaghan 1997 for the particular form given in (2.27), and Monaghan 1992 for an alternative $`\mathrm{\Pi }_{ab}`$).

Accordingly, we add the term

$$\frac{1}{2}\underset{b}{}m_b\mathrm{\Pi }_{ab}\overline{_aW}_{ab},$$
(2.28)
to the right hand side of (2.10) where

$$\overline{_aW}_{ab}=\frac{1}{2}(_aW_{ab}(h_a)+_aW_{ab}(h_b)).$$
(2.29)
Angular and linear momentum are still conserved because the viscous force is along the line joining the centres of the SPH particles. Following the argument leading to (2.26) we find

$$\frac{du_a}{dt}=\frac{P_a}{\rho _a^2\mathrm{\Omega }_a}\underset{b}{}m_b𝐯_{ab}_aW_{ab}(h_a)\frac{1}{2}\underset{b}{}m_b\mathrm{\Pi }_{ab}𝐯_{ab}\overline{_aW}_{ab}.$$
(2.30)
We can write $`\overline{_aW}_{ab}`$ in the form $`𝐫_{ab}F_{ab}`$, where $`F_{ab}0`$. It then follows that the viscous contribution to $`du/dt`$ is $`0`$. It is also easy to show that when the fluid is viscous the entropy change is $`0`$.

There is another approach to the effect of dissipation which is closely analogous to that used in the simulation of liquids at normal temperatures and pressures. In these cases the entropy change can usually be neglected and the dominant effect of the dissipation is to reduce the kinetic energy to zero. This is the case for water initially sloshing in a laboratory tank. In this case we write (2.24), with the dissipation term added, in the form

$$\frac{dE^{}}{dt}=\frac{d}{dt}(E_k+E_G+E_{\mathrm{𝑡ℎ}})=\frac{1}{2}\underset{b}{}m_b\mathrm{\Pi }_{ab}𝐯_{ab}\overline{_aW}_{ab},$$
(2.31)
where only the pressure terms are included in the rate of change of thermal energy. Because the right hand side of (2.31)is negative definite the SPH equations guarantee that dissipation will cause $`E^{}`$ to decrease as it should. For example, an oscillating polytrope simulated keeping the entropy of every particle constant, will evolve to a state where $`E^{}`$ is a minimum. When it reaches this state $`E_K`$ will be zero.

## 3 Averaging the velocity field

In the standard SPH calculations the the position of particle $`a`$ is found by integrating

$$\frac{d𝐫_a}{dt}=𝐯_a.$$
(3.1)
However, it is possible for particles to stream through each other unless the viscosity is large. For that reason the XSPH version of SPH was suggested (Monaghan 1989) where the particles were moved with with a smoothed velocity, but the acceleration equation was unchanged. The smoothed velocity $`\widehat{𝐯}_a`$ is defined by an average over the velocities of neighbouring particles according to

$$\widehat{𝐯}_a=𝐯_a+ϵ\underset{b}{}\frac{m_b}{\overline{\rho }_{ab}}(𝐯_b𝐯_a)K(a,b),$$
(3.2)
where $`ϵ`$ is a parameter that is typically chosen to be 0.5, $`\overline{\rho }_{ab}`$ denotes an average density, for example

$$\frac{1}{\overline{\rho }_{ab}}=\frac{1}{2}\left(\frac{1}{\rho _a}+\frac{1}{\rho _a}\right),$$
(3.3)
and $`K(a,b)=K(b,a)`$ is a smoothing kernel. The kernel $`K(a,b)`$ does not have to be the same as the kernel used in other SPH equations. In this paper we assume the kernel has the form

$$K(a,b)=\frac{1}{2}\left(W_{ab}(h_a)+W_{ab}(h_b)\right).$$
(3.4)
The SPH particles are moved according to

$$\frac{d𝐫_a}{dt}=\widehat{𝐯}_a.$$
(3.5)
The smoothed velocity has a number of useful properties. First we note that

$$\underset{a}{}m_a𝐯_a=\underset{a}{}m_a\widehat{𝐯}_a,$$
(3.6)
which shows that the linear momentum is the same whether computed with $`\widehat{𝐯}`$ or $`𝐯`$. Second the centre of mass moves with the $`𝐯`$ because

$$\frac{d}{dt}\underset{a}{}m_a𝐫_a=\underset{a}{}m_a𝐯_a.$$
(3.7)
Third, the angular momentum is unaffected because

$$\underset{a}{}m_a𝐫_a\times \frac{d𝐯_a}{dt}=\frac{d}{dt}\underset{a}{}m_a𝐫_a\times 𝐯_a.$$
(3.8)
Finally we note that the motion can be exactly reversed by reversing $`𝐯`$ and $`\widehat{𝐯}`$.

The continuum limit of (3.2) is

$$\widehat{𝐯}_a=𝐯_a+\frac{ϵ}{2\rho _a}(\rho (𝐫)+\rho _a)(𝐯(𝐫^{})𝐯(𝐫_𝐚))K(𝐫^{},𝐫_a)\mathrm{𝐝𝐫}^{},$$
(3.9)
where

$$K(𝐫^{},𝐫_𝐚)=\frac{1}{2}(W(𝐫𝐫_a,h_a)+W(𝐫𝐫_a,h(\rho (𝐫^{}))$$
(3.10)
Taylor expansion about $`𝐫_a`$ followed by integration gives the approximation

$$\widehat{v}_a^j=v_a^j+\frac{\alpha ^2}{\rho _a}(\rho v^j)+\frac{1}{2}(\rho v^j)\frac{H_a}{\rho _a}\frac{\alpha ^2}{h_a},$$
(3.11)
where

$$\alpha ^2=\frac{ϵ}{2}(x^{}x_a)^2W(𝐫^{}𝐫_a,h_a)\mathrm{𝐝𝐫}^{},$$
(3.12)
defines the $`\alpha `$ in the alpha model although, in our analysis, $`\alpha `$ may be different for each particle and, because $`h_a`$ can vary with time, so can $`\alpha `$. If $`\rho `$ is constant (3.11) reduces to

$$\widehat{v}_a^j=v_a^j+\alpha ^2(^2v^j)_a,$$
(3.13)
as in the alpha model.

For typical SPH kernels $`\alpha ^2h^2`$. Accordingly, provided the scale of variation of the velocity field is much greater than $`h`$, the approximation (3.11) will be satisfactory. However, our original XSPH smoothing (3.2) applies regardless of the scale of the velocity field variation. In particular we expect that it can be used even when shocks occur.

The results of this section show that the XSPH smoothing has a number of attractive features. It prevents SPH particles passing through each other, it retains the conservation of linear and angular momentum, and it agrees with the alpha model when the density is constant. We now show that the equations of motion can be obtained from a Lagrangian which is a function $`𝐫`$ and $`\widehat{𝐯}`$.

## 4 Lagrange’s equations for XSPH

The analysis of Holm (1999) suggests that the Lagrangian will take the form

$$L=\underset{b}{}m_b\left(\frac{1}{2}\widehat{𝐯}_b𝐯_bu(\rho _b,s_b)\right),$$
(4.1)
where the canonical coordinates for particle $`b`$ are $`𝐫_b`$ with rate of change $`\widehat{𝐯}_b`$. In the interests of simplicity we have not included gravitational forces. These can be easily included in the final equations.

The kinetic energy part of this Lagrangian can be rewritten by inverting (3.2) to give

$$𝐯_a=\widehat{𝐯}_aϵ\underset{b}{}\frac{m_b}{\overline{\rho }_{ab}}(\widehat{𝐯}_b\widehat{𝐯}_a)W_{ab},$$
(4.2)
where the errors in the inversion are of $`O(h^4)`$. We then get

$$\frac{1}{2}\underset{b}{}m_b\widehat{𝐯}_b𝐯_b=\frac{1}{2}\underset{b}{}m_b\widehat{𝐯}_b\widehat{𝐯}_b+\frac{ϵ}{4}\underset{k}{}\underset{b}{}\frac{m_bm_k}{\overline{\rho }_{bk}}(\widehat{𝐯}_b\widehat{𝐯}_k)^2W_{kb}$$
(4.3)
which is positive definite.

If the sums over $`k`$ are approximated by integrals, and we assume the density is constant, Taylor series expansion of the velocity difference in the integrand, followed by integration, shows that

$$\underset{k}{}\frac{m_k}{\overline{\rho }_{ab}}(\widehat{𝐯}_b\widehat{𝐯}_a)^2W_{kb}=(\widehat{v}_b^i\widehat{v}_b^i)(x_bx)^2W(𝐫_b𝐫^{})𝑑\tau ^{},$$
(4.4)
where $`\widehat{v}_a^i`$ denotes the $`i`$ th component of $`\widehat{𝐯}_a`$, and there is an implied summation over $`i`$. Substituting this result into (4.3) gives

$$\frac{1}{2}\underset{b}{}m_b(\widehat{v}_b^2+\alpha ^2\widehat{v}_b^j.\widehat{v}_b^j).$$
(4.5)
In the absence of dissipation this quantity is an invariant for an incompressible fluid (Holm 1999) and is the alpha equivalent of the kinetic energy. Note particularly the $`v`$ terms which play a key role in the distribution of energy. The SPH equivalent is the sum involving the square of the velocity differences. This sum only involves neighbour particles, and is small if they have similar velocities. It becomes large only when the velocities of the neighbouring particles differ in sign or magnitude or both. This occurs when there is substantial change in $`v`$ on the scale of $`h`$. Such a change is equivalent to the existence of high wave number modes.

Returning to Lagrange’s equations the canonical momentum for particle $`a`$ is

$$\frac{L}{\widehat{𝐯}_a}=m_a\widehat{𝐯}_a+\frac{ϵ}{2}\underset{k}{}\underset{b}{}\frac{m_bm_k}{\overline{\rho }_{bk}}(\widehat{𝐯}_b\widehat{𝐯}_k)(\delta _{ba}\delta _{ka})W_{kb},$$
(4.6)
which reduces to

$`{\displaystyle \frac{L}{\widehat{𝐯}_a}}`$ $`=`$ $`m_a\widehat{𝐯}_aϵ{\displaystyle \underset{k}{}}{\displaystyle \frac{m_am_k}{\overline{\rho }_{ak}}}(\widehat{𝐯}_k\widehat{𝐯}_a)W_{ka},`$ (4.7)
$`=`$ $`m_a𝐯_a`$ (4.8)
so that the canonical momentum of the Lagrangian defined with $`\widehat{𝐯}`$ is the normal momentum constructed with $`𝐯`$. This dual relation between $`𝐯`$ and $`\widehat{𝐯}`$ reminds us of the relationship between Eulerian and Lagrangian velocities noted by Holm (1999).

The remaining term required for the equation of motion of particle $`a`$ is $`L/𝐫_a`$. The details of this calculation are given in the Appendix 1. We find

$$\frac{d𝐯_a}{dt}=\underset{b}{}m_b\left(\left(\frac{P_a}{\mathrm{\Omega }_a\rho _a^2}\frac{ϵ}{4}\left(\frac{v_{ab}^2}{\overline{\rho }_{ab}}\frac{\zeta _a}{\mathrm{\Omega }_a\rho _a^2}+\frac{2\nu _aH_a^{}}{\rho _a\mathrm{\Omega }_a}\right)\right)_aW_{ab}(h_a)+(ab)\right),$$
(4.9)
where $`(ab)`$ denotes terms that are identical to the previous but $`a`$ and $`b`$ are interchanged except in $`_a`$, and the other quantities, where not previously defined, have the following meaning

$$v_{ab}^2=(𝐯_a𝐯_b)^2,$$

$$\zeta _a=\frac{1}{2}\underset{k}{}m_kv_{ak}^2(W_{ak}(h_a)+W_{ak}(h_k)),$$

and

$$\nu _a=\underset{k}{}m_k\frac{v_{ak}^2}{\overline{\rho }_{ak}}\frac{W_{ak}(h_a)}{h_a}.$$

The equations required for the SPH alpha model are the acceleration equation (5.9), and the density equation either in the summation form (4.2) or the continuity equation

$$\frac{d\rho _a}{dt}=\frac{1}{\mathrm{\Omega }_a}\underset{b}{}m_b\widehat{𝐯}_{ab}_aW_{ab},$$
(4.10)
together with $`d𝐫_a/dt=\widehat{𝐯}_a`$. The derivative following the motion is given by

$$\frac{d}{dt}=\frac{}{t}+\widehat{𝐯}.$$
(4.11)
These beautifully simple equations comprise the SPH alpha model. They contain all the desirable features of Holm’s (1999) incompressible flow equations but generalised to compressible flow. If the fluid is self gravitating then the gravity terms can be added as in $`\mathrm{\S }2`$. These new equations differ from those in $`\mathrm{\S }2`$, for example (2.10), because the velocity averaging leads to velocity dependent terms in the force and the particles are moved with the smoothed velocity. Otherwise the equations are very similar.

In order to compare our results with that of the incompressible alpha model we assume $`h`$ and $`\rho `$ are constant. In this case $`\mathrm{\Omega }=1`$ and the term involving $`\nu `$ vanishes because $`h`$ is zero. The equation of motion becomes

$$\frac{d𝐯_a}{dt}=\underset{b}{}m_b\left(\left(\frac{P_a}{\rho _a^2}\frac{ϵ}{4}\left(\frac{v_{ab}^2}{\overline{\rho }_{ab}}\frac{\zeta _a}{\rho _a^2}\right)\right)_aW_{ab}(h_a)+(ab)\right).$$
(4.12)
If the sums in the SPH equations are converted to integrals the dominant terms give the continuum alpha model (for the details see Appendix 2). In particular, we recover Holm’s acceleration equation (Holm 1999, equation (143))

$$\frac{v^j}{t}+(\widehat{𝐯})v^j+v^{\mathrm{}}\frac{\widehat{v}^{\mathrm{}}}{x^j}=\frac{1}{\rho }\frac{}{x^j}\left(P\frac{1}{2}\widehat{𝐯}\widehat{𝐯}\frac{\alpha ^2}{2}\frac{\widehat{v}^k}{x^{\mathrm{}}}\frac{\widehat{v}^k}{x^{\mathrm{}}}\right).$$
(4.13)
We have therefore achieved our aim of constructing a Lagrangian based turbulence model which reduces to the alpha model in the continuum limit. In the following sections we will investigate some of the properties of the SPH model.

## 5 Conservation laws

In the absence of boundaries and external forces the Lagrangian is invariant to translations and rotations. Linear and angular momentum are therefore conserved. Because the Lagrangian has no explicit time dependence there is an invariant

$`E`$ $`=`$ $`{\displaystyle \underset{a}{}}\widehat{𝐯}_a{\displaystyle \frac{L}{\widehat{𝐯}_a}}L,`$ (5.1)
$`=`$ $`{\displaystyle \underset{a}{}}m_a\left({\displaystyle \frac{1}{2}}\widehat{𝐯}_a𝐯_a+u(\rho ,s)\right).`$ (5.2)
which can be appropriately called the energy. Using (5.3) we can write

$$E=\underset{b}{}m_b\left(\frac{1}{2}\widehat{𝐯}_b\widehat{𝐯}_b+\frac{ϵ}{4}\underset{k}{}\frac{m_k}{\overline{\rho }_{bk}}(\widehat{𝐯}_b\widehat{𝐯}_k)^2W_{kb}+u(\rho _b,s_b)\right).$$
(5.3)
If the entropy is constant then the system is also invariant to the necklace transformation considered earlier ($`\mathrm{\S }2.2`$. In the present case the particles, assumed to have the same mass and entropy, are shifted to the adjacent position on the loop and given the $`\widehat{𝐯}`$ of the new position.

The change in $`L`$ can be written

$$\delta L=\underset{j}{}\left(\frac{L}{𝐫_j}\delta 𝐫_j+\frac{L}{\widehat{𝐯}_j}\delta \widehat{𝐯}_j\right),$$
(5.4)
where $`j`$ denotes the label of a particle on the loop. The change in position and velocity are given by

$$\delta 𝐫_j=𝐫_{j+1}𝐫_j,$$
(5.5)
and

$$\delta \widehat{𝐯}_j=\widehat{𝐯}_{j+1}\widehat{𝐯}_j.$$
(5.6)
Using Lagrange’s equations we can rewrite (5.4) in the form

$$\frac{d}{dt}\left(\underset{j}{}\frac{L}{\widehat{𝐯}_j}(𝐫_{j+1}𝐫_j)\right)=0,$$
(5.7)
and recalling that the particle masses are identical, we deduce

$$\frac{d}{dt}\underset{j}{}𝐯_j(𝐫_{j+1}𝐫_j)=0.$$
(5.8)
which is identical to the discrete circulation theorem (2.16). The reader will note that this result depends on the fact that the canonical momentum of particle $`j`$ is $`m_j𝐯_j`$.

## 6  Is there a steady state ?

The SPH alpha model, in the absence of dissipation, is a conservative mechanical system derived from a Lagrangian. If the system is at rest, for example in a uniform periodic box, then small disturbances would propagate as sound waves. In the linear approximation the waves are not coupled so that energy put into a mode stays there. If the motion becomes non linear the waves are coupled and any is shared amongst the modes. If equilibrium statistical mechanics was valid for this system each mode would have the same energy, and the distribution of energy between wave number $`k`$ and $`k+dk`$ would be proportional to $`k^2`$ in three dimensions. However, it is known from other cases, for example the Fermi-Past-Ulam study of coupled non linear oscillators (Fermi et. al 1955, Segr$`\stackrel{`}{e}`$), that a non linear non-dissipative system may not evolve to statistical equilibrium, but instead moves energy back and forth between a relatively small number of modes.

In the present case we have non-linear coupling through both the velocity terms in the equation of motion, and the density variations with velocity. The results of Chen et. al (1999) for the continuum alpha model suggest that in wave number space the non-linear terms act to turn the spectrum down for large wave numbers. It is as if velocity terms in the equation of motion acted as a potential preventing a piling up of energy at the high wave numbers. This suggests that the non-dissipative model transfers energy back and forth between high and low wave numbers.

The importance of this for turbulence is the following. It is known that a direct numerical simulation without dissipation leads to a piling up of energy in the high wave number modes. If dissipation is included, but it is very weak at the resolution of the numerical simulation, then the piling up will still take place, though eventually, if there is no input of energy it will decay. If there is an input of energy at the low wave numbers, and it is sufficiently large, then the piling up will continue because it is arriving faster than it can leave. The idea of sub-grid models is to prevent this piling up by including extra dissipation. The alpha model has the ability, without dissipation, to reduce the piling up of energy at high wave numbers. For this to work it is necessary for the system to transfer energy back and forth between high and low wave numbers. As noted earlier, in three dimensional turbulence there is transfer both ways (Piomelli et al. 1991, Jimenez 1994, Pope 2000, Woodward et al. 2001) and, at any time, parts of the domain may have transfer up, and parts may have transfer down. The downward transfer is the residual between two much larger oppositely directed fluxes. The alpha model appears to mimic these features. How this occurs requires a numerical investigation of the non-dissipative system and this will be given elsewhere.

If energy is put into the system then it is still necessary to introduce dissipation but it can be very much weaker than the dissipation used in the usual sub-grid models.

## 7 Dissipation

The previous analysis concerns the non-dissipative fluid. If it is correct that the non-dissipative SPH model transfers energy between low and high wave numbers to prevent the piling up of energy at high wave numbers, then it should be possible to add a small amount of dissipation to remove energy at the resolution scale ($`h`$ in the SPH calculation) so that, with an energy input at low wave numbers, the system reproduces the Kolmogorov spectrum over all but the high wave number modes.

The dissipation term for SPH is based on $`\mathrm{\Pi }_{ab}`$ (see (2.27)). In the present case where the particles are moved with $`\widehat{𝐯}`$ we replace the velocity $`𝐯`$ in $`\mathrm{\Pi }_{ab}`$ with $`\widehat{𝐯}`$, and substitute the term in the right hand side of (4.9) to obtain the dissipative equations. For convenience we define

$$L_K=\underset{a}{}m_a\frac{1}{2}\widehat{𝐯}_a𝐯_a,$$

then (4.9) can be written

$$\frac{d}{dt}\left(\frac{L_K}{\widehat{𝐯}_a}\right)\frac{L_K}{𝐫_a}+m_a\underset{b}{}m_b\left(\frac{P_a}{\mathrm{\Omega }_a\rho _a^2}_aW_{ab}(h_a)+(ab)\right)=m_a\underset{b}{}m_b\mathrm{\Pi }_{ab}\overline{_aW}_{ab}.$$
(7.1)
We can deduce the thermal energy equation by requiring the total energy to be constant. Using the same arguments as in $`\mathrm{\S }2.3`$, but now taking the dot product of the acceleration equation with $`\widehat{𝐯}_a`$, we get the following equation for the thermal energy of particle

$$\frac{du_a}{dt}=\frac{P_a}{\rho _a^2\mathrm{\Omega }_a}\underset{b}{}m_b\widehat{𝐯}_{ab}_aW_{ab}(h_a)\frac{1}{2}\underset{b}{}m_b\mathrm{\Pi }_{ab}\widehat{𝐯}_{ab}\overline{_aW}_{ab}.$$
(7.2)
The term involving the dissipation is positive definite so that the SPH dissipation increases the thermal energy as it should. It is easy to deduce that the dissipation also guarantees that the entropy increases for each particle.

In the same way $`E^{}`$, the quasi energy appropriate for systems computed assuming the entropy of each particle is constant (see $`\mathrm{\S }2.3`$), changes with time according to

$$\frac{dE^{}}{dt}=\underset{a}{}\underset{b}{}m_am_b\mathrm{\Pi }_{ab}\widehat{𝐯}_a\overline{_aW}_{ab},$$
(7.3)
so that, as in the case without velocity averaging, the system will evolve to a state with $`E^{}`$ a minimum.

In order to establish these results it has been necessary to use the smoothed velocity $`\widehat{𝐯}`$ instead of $`𝐯`$ in the viscosity term $`\mathrm{\Pi }_{ab}`$. We can therefore expect that the SPH viscosity will be smaller than usual in those parts of the flow where the velocity is smooth. In SPH calculations the viscosity term $`\mathrm{\Pi }`$ has a coefficient (also denoted by $`\alpha `$) which can be assigned a value $`1`$. If desired, each SPH particle can have its own $`\alpha `$ which can then be determined by the dynamics (Morris and Monaghan 1997) in such a way that its value is $`0.1`$ away from shocks. If the smoothed velocity is used in $`\mathrm{\Pi }_{ab}`$, and the coefficient is varied as described by Morris and Monaghan (1997), then the coefficient away from shocks can be expected to be reduced and values of $`0.01`$ appear possible. Preliminary investigations of polytrope oscillations and Kelvin-Helmholtz instabilities show that this conjecture is true in these cases. A full investigation of these effects will be given elsewhere.

## 8 Scaling and Turbulence

The previous sections have laid out the formalism of the SPH alpha model of turbulence. If there is no dissipation we have argued that there is a back and forth flow of energy between low and high wave numbers. If there is an input of energy at low wave numbers, and the system has a small dissipation of the type considered above, then we conjecture that the SPH system will reproduce the Kolmogorov spectrum for the low wave number, energy containing modes. Support for this conjecture is provided by the continuum calculations of Chen et al. (1999).

Some features of this process can be guessed by scaling arguments (Kolmogorov and Oboukov as described by Landau and Lifshitz (1993) which assume (a) that the details of the dissipation are irrelevant to the low wave number modes and (b) the rate of loss of energy through the various length scales is constant.

We start with the energy

$$E=\underset{b}{}m_b\left(\frac{1}{2}\widehat{𝐯}_b\widehat{𝐯}_b+\frac{ϵ}{4}\underset{k}{}\frac{m_k}{\overline{\rho }_{bk}}(\widehat{𝐯}_b\widehat{𝐯}_k)^2\overline{W}_{kb}+u(\rho _b,s_b)\right).$$
(8.4)
and assume there is a specified rate of loss of energy from the large scale motions $``$ and this cascades through the various length scales until it is dissipated by either by our weak SPH viscosity or by molecular viscosity.

The rate of change of energy per unit mass $``$ determined by the first term in $`E`$ for velocity variations with characteristic length $`\mathrm{}`$ is

$$=\frac{v_{\mathrm{}}^3}{\mathrm{}}(1+\frac{\alpha ^2}{\mathrm{}^2},)$$
(8.5)
where the term involving $`\alpha ^2`$ can be estimated from the continuum limit of the energy which introduces the term $`\alpha ^2v^jv^j`$ (see (4.5)). Because $`\alpha `$ varies with the resolution length the following scaling formula should be interpreted in terms of the local conditions of the gas.

The characteristic velocity is

$$v_{\mathrm{}}=\left(\frac{\mathrm{}^3}{\mathrm{}^2+\alpha ^2}\right)^{\frac{1}{3}}.$$
(8.6)
If $`\mathrm{}\alpha `$ then this reduces to the usual Kolmogorov velocity scaling $`v_{\mathrm{}}=\mathrm{}^{1/3}^{1/3}`$. If $`\mathrm{}\alpha `$

$$v_{\mathrm{}}=\mathrm{}\left(\frac{}{\alpha ^2}\right)^{\frac{1}{3}},$$
(8.7)
The velocity therefore reduces more rapidly with $`\mathrm{}`$ in the smoothing region. The turnover time is given by

$$\tau =\frac{\mathrm{}}{v_{\mathrm{}}}=\left(\frac{\mathrm{}^2+\alpha ^2}{}\right)^{\frac{1}{3}}.$$
(8.8)
which is approximately constant when $`\mathrm{}\alpha `$. This shows that all length scales have essentially the same turnover time when the smoothing dominates. These effects on the velocity scale and the turnover time are not due to dissipation. They occur because of the quadratic velocity terms in the acceleration equation.

If the energy per unit mass in the wave number range $`k`$ to $`k+dk`$ is $`E(k)dk`$ then $`E(k)`$ has the dimensions of $`(length)^3/(time)^2`$. We can relate this to the energy dissipation and the wave number. For large length scales we expect that $`E(k)`$ will vary as $`^rk^n`$. We then find

$$E(k)=\frac{^{\frac{2}{3}}}{k^{\frac{5}{3}}}.$$
(8.9)
For the short length scales we relate $`E(k)`$ to that part of $``$ which depends on the characteristic velocity and length scales. Referring to (8.5) this is

$$\eta =\frac{v_{\mathrm{}}^3}{\mathrm{}^3}.$$
(8.10)
We then find

$$E(k)=\frac{\eta ^{\frac{2}{3}}}{k^3}.$$
(8.11)
The energy spectrum therefore steepens for small length scales. These approximate results have been confirmed for the continuum case by numerical experiments (Chen et al. 1999) but a similar extensive study is required to confirm these arguments for the SPH alpha model.

## 9  Conclusions

The results of this paper can be summarised as follows:

1. An averaged or smoothed velocity can be constructed for each particle such that, if each particle moves with this velocity, the linear and angular momentum are unchanged from the values with the unsmoothed velocity. This part of our analysis is the same as in the XSPH algorithm.
2. The smoothing of the velocity agrees with that used in the alpha model of turbulence in the continuum limit with the density constant.
3. A particle Lagrangian can be constructed using as canonical variables the smoothed velocity $`\widehat{𝐯}`$ and $`𝐫`$. This Lagrangian has the same structure as the continuum Lagrangian in the alpha model.
4. The SPH Lagrangian leads to equations which conserve energy, linear and angular momentum and satisfy the same discrete Kelvin circulation theorem as for SPH with no smoothing of the velocity.
5. In the continuum limit the new SPH equations agree with the equations of the alpha turbulence model when the density is constant.

In addition we conjecture that:

1. The non dissipative case should prevent energy piling up at the highest wave numbers allowed by the resolution and it does so by transferring energy back and forth between low and high wave numbers. This conjecture is supported by scaling arguments, by the fact that this occurs in laboratory turbulence, and by analogies to dynamical systems like the Fermi-Past-Ulam problem.
2. Weak dissipation should be adequate to reproduce the Kolmogorov spectrum for isotropic homogeneous turbulence. This conjecture is supported by the calculations for the continuum case by Chen et al. (1999).
3. The artificial viscosity used in SPH calculations can be very much smaller than that currently used in smooth parts of the flow provided each particle has its own viscosity coefficient. This conjecture is supported by preliminary calculations of the oscillations of a polytrope and the growth of Kelvin Helmholtz instabilities.

This paper has concentrated on the dynamics of the SPH system with and without viscous dissipation. We have had in mind those turbulent processes which arise from mechanical sources such as Kelvin Helmhotz instabilities in jets. We have not considered the more complex case of turbulence produced by thermal effects as in turbulent convection. In particular we have not considered how the smoothing of the velocity field affects the transport of heat. The formulation of the alpha model gives us two velocity fields and from these it appears possible to model the increased diffusivity expected from turbulence. How this should be done, however, is not clear.

Another important issue in turbulence dynamics is mixing (Warhaft 2000). The problem is often counter-intuitive, with many seemingly laminar flows leading to significant mixing (Aref et al. 1989). The natural approach to mixing is the Lagrangian formulation and it is possible that the SPH formulation provides an efficient practical means of determining the extent of mixing both in turbulent and in nearly laminar flow. This is another area of research which is opened up by the theoretical developments in this paper.

## 10  Acknowledgment

Numerous, pleasant and thought provoking interactions with Darryl Holm have greatly improved this paper.

## 11  References

1. Aref, H., Jones, S.W., Mofina, S., and Zadwadzki, I. Physica D, 37, 423, (1989)
2. Bonet, J., and Lok, T-S,L. Comp. Meth. App. Mech. Engrg180, 97, (1999)
3. Bonet, J. *Meshless Methods*, Bonn (2001).
4. Bretherton, F. J. Fluid Mech., 44, 19, (1970).
5. Chen, S., Holm, D. D., Margolin, L. G., and Zhang, R. Physica D, 133, 66 , (1999)
6. Eckart, C. Phys. Fluids, 3, 421, (1960)
7. Fermi,E., Pasta, J., and Ulam, S. Studies of Non Linear Problems. LA-1940 (1955)
8. Feynman, R. Prog. Low Temp. Phys., 1, 17, (1957)
9. Ghosal, S. AIAA J., 37, 425, (1999).
10. Holm, D. D. Physica D,133, 215, (1999)
11. Jim$`\stackrel{`}{e}`$nez, J. Fluid Physics. Lecture Notes on Summer Schools. Ed. M.G. Verlarde and C.I. Christov. Publ. World Scientific. (1995).
12. Landau, L.D., and Lifshitz, E.M. Fluid Mechanics Vol 6, Course of Theoretical Physics. (1993). Publ. Pergamon.
13. Monaghan, J.J. Jour. Computat. Phys., 64, 2, (1989)
14. Monaghan, J.J. Ann. Rev. Astron. Ap., 30, 543, (1992)
15. Monaghan, J.J.Jour. Computat. Phys., 110, 399, (1994)
16. Monaghan, J.J., and Kos, A., J. Waterways, Ports, Coastal and Ocean Eng., 125, 145, (1999)
17. Monaghan, J. J., and Price, D.L. Mon. Not. Roy. Astr. Soc.328, 381, (2001)
18. Morris, J.P., and Monaghan, J.J. J. Computat. Phys.,136, 41, (1997).
19. Piomelli, U., Cabot, W.H., Moin, P., and Lee, S. Phys. Fluids A 3,1766, (1991).
20. Pope, S.B. *Turbulent Flows*. Publ. Cambridge University Press. (2000).
21. Woodward, P., Porter. D. H., Sytine, I., Anderson, S. E., Mirin, A. A., Curtis, B. C., Cohen, R. H., Dannevik, W. P., Dimits, A. M., Eliason, D. E., Winkler, K-H., and Hodson S. W.  *Computational Fluid Dynamics*, 4th UNAM Supercomputing Conference. Ed. E. Ramos and G. Cisneros. Publ. World Scientific. (2001).
22. Salmon, R. Ann. Rev. Fluid Mech. , 20, 225, (1988)
23. E.Segr$`\stackrel{`}{e}`$ Ed.  *Enrico Fermi Collected Papers. Vol II*. page 266. (1965). Univ. of Chicago Press.
24. Springel, V., and Hernquist, L. Mon. Not. Roy. Astr. Soc (2001).
25. Warhaft, Z.Ann. Rev. Fluid Mech. , 32, 203, (2000).

## Appendix A The $`W/h`$ terms

The contribution to the Lagrangian from the thermal energy gives the usual terms involving the pressure. The velocity terms in the Lagrangian are

$$L_v=\frac{ϵ}{4}\underset{b}{}\underset{k}{}m_bm_k\frac{v_{bk}^2}{\overline{\rho }_{bk}}\overline{W}_{bk},$$
(A.1)
where

$$\overline{W}_{bk}=\frac{1}{2}(W_{bk}(h_k)+W_{bk}(h_b)).$$
(A.2)
By a change of label it is easy to show that $`\overline{\rho }_{bk}`$ can be replaced by $`\rho _b`$ in (A.1). Then

$$\frac{L_v}{𝐫_a}=\frac{ϵ}{4}\underset{b}{}\underset{k}{}m_bm_kv_{bk}^2\left(\overline{W}_{bk}\frac{}{𝐫_a}\left(\frac{1}{\rho _b}\right)+\frac{1}{\rho _b}\frac{\overline{W}_{bk}}{𝐫_a}\right).$$
(A.3)
Making use of the earlier expression for the density gradient (2.7), and some straightforward simplification, the first term in (A.3) can be written

$$\frac{m_aϵ}{4}\underset{b}{}\left(\frac{\zeta _a}{\rho _a^2\mathrm{\Omega }_a}_aW_{ab}(h_a)+\frac{\zeta _b}{\rho _b^2\mathrm{\Omega }_b}_aW_{ab}(h_a)\right),$$
(A.4)
where

$$\zeta _a=\underset{k}{}m_kv_{ak}^2\overline{W}_{ak}.$$
(A.5)
The second term in (A.3) is more complicated. We first consider the contribution from $`W_{bk}(h_b)`$ in $`\overline{W}_{bk}`$. We find

$$\frac{ϵ}{4}\underset{b}{}\underset{k}{}m_bm_k\frac{v_{bk}^2}{\rho _b}(_bW_{bk}(\delta _{ab}\delta _{ka})+\frac{W_{bk}}{h_b}H_b^{}\frac{\rho _b}{𝐫_a}.)$$
(A.6)
with a similar term from the $`W_{bk}(h_k)`$ in $`\overline{W}_{bk}`$. Combining these terms the contribution not involving $`W/h`$ is

$$\frac{m_aϵ}{4}\underset{b}{}m_b\frac{v_{ab}^2}{\overline{\rho }_{ab}}\left(_aW_{ab}(h_a)+_aW_{ab}(h_b)\right).$$
(A.7)
The $`W/h`$ terms give the terms involving $`H^{}`$ in (4.9).

## Appendix B  The continuum acceleration equation

From (4.9) the first term we consider is

$$\frac{ϵ}{2}\underset{b}{}m_b\frac{(\widehat{𝐯}_b\widehat{𝐯}_a)^2}{\overline{\rho }_{ab}}_aW_{ab}.$$
(B.1)
To simplify the analysis we only consider the two dimensional case and assume the density is constant. Without loss of generality we consider only the $`x`$ component.

The continuum limit of the summation is

$$\frac{ϵ}{2}[𝐯(𝐫)𝐯(𝐫^{})]^2\frac{W}{x^{}}𝑑x^{}𝑑y^{}.$$
(B.2)
Expanding $`𝐯(𝐫^{})`$ about $`𝐫`$ to second order, and only retaining terms of odd parity in $`(x^{}x)`$ and even parity in $`(y^{}y)`$, (B.2) becomes

$$\frac{ϵ}{2}\left(\mathrm{\Delta }x^3\frac{^2v^x}{^2x^2}\mathrm{\Delta }x\mathrm{\Delta }y^2\frac{^2v^x}{^2y^2}2\mathrm{\Delta }x\mathrm{\Delta }y^2\frac{v^x}{y}\frac{^2v^x}{xy}+(v^xv^y)\right)\frac{W}{x^{}}𝑑x^{}𝑑y^{},$$
(B.3)
where $`v^x`$ and $`v^y`$ denote the $`x`$ and $`y`$ components of $`𝐯`$ and $`\mathrm{\Delta }x`$ and $`\mathrm{\Delta }y`$ denote $`(x^{}x)`$ and $`(y^{}y)`$ respectively. The notation $`v^xv^y`$ means include a set of terms identical to the previous but with $`v^x`$ replaced by $`v^y`$. After integration by parts and noting that

$$(xx^{})^2W𝑑x^{}𝑑y^{}=(yy^{})^2W𝑑x^{}𝑑y^{},$$
(B.4)
(B.2) becomes

$$\alpha ^2\left(3\frac{v^x}{x}\frac{^2v^x}{x^2}+\frac{v^x}{x}\frac{^2v^x}{y^2}+2\frac{v^x}{x}\frac{^2v^x}{xy}+(v^xv^y)\right).$$
(B.5)
Combining terms we can write the previous expression as

$$\alpha ^2\left(\frac{v^x}{x}^2v^x+\frac{v^y}{x}^2v^y+\frac{}{x}(v^xv^x+v^yv^y)\right).$$
(B.6)
or, on replacing $`v^x`$ and $`v^y`$ by $`v^1`$ and $`v^2`$ respectively

$$\alpha ^2(\frac{v^{\mathrm{}}}{x}^2v^{\mathrm{}}+\frac{}{x}(v^{\mathrm{}}v^{\mathrm{}}).$$
(B.7)
with summation on repeated indices.

The remaining contribution from (4.9) is

$$\frac{ϵ}{4}\underset{b}{}m_b\left(\frac{\zeta _a}{\rho _a^2\mathrm{\Omega }_A}+\frac{\zeta _b}{\rho _b^2\mathrm{\Omega }_b}\right)_aW_{ab}(h).$$
(B.8)
The continuum expression for $`\zeta `$ is

$$\zeta =\rho (𝐫^{})[𝐯(𝐫)𝐯(𝐫^{})]^2W(𝐫𝐫^{})\mathrm{𝐝𝐫}^{}.$$
(B.9)
Expanding the velocity difference in a Taylor series as before, we get the approximation

$$\zeta =2\alpha ^2v^{\mathrm{}}v^{\mathrm{}}.$$
(B.10)
Furthermore, noting $`\rho `$ is constant,

$$\underset{b}{}m_b\frac{\zeta _a}{\rho ^2}_aW_{ab}\frac{\zeta _a}{\rho }_𝐫W\mathrm{𝐝𝐫}^{}.$$
(B.11)
This integral vanishes because it has odd symmetry. The remaining term is

$$\underset{b}{}m_b\frac{\zeta _b}{\rho ^2}_aW_{ab}\frac{_a\zeta _a}{\rho }.$$
(B.12)
Multiplying (B.12) by $`ϵ/4`$ and combining with (12.7) gives

$$\alpha ^2(\frac{v^{\mathrm{}}}{x}^2v^{\mathrm{}}+\frac{1}{2}\frac{}{x}(v^{\mathrm{}}v^{\mathrm{}}).$$
(B.13)
which, apart from the pressure term, gives the continuum form of the SPH force/mass.

The relevant term in Holm’s equation (see (4.22)) is

$$v^{\mathrm{}}\frac{\widehat{v}^{\mathrm{}}}{x}+\frac{}{x}\left(\frac{1}{2}\widehat{𝐯}\widehat{𝐯}+\frac{\alpha ^2}{2}\frac{\widehat{v}^j}{x^{\mathrm{}}}\frac{\widehat{v}^j}{x^{\mathrm{}}}\right),$$
(B.14)
which can be written

$$(\widehat{v}^{\mathrm{}}v^{\mathrm{}})\frac{\widehat{v}^{\mathrm{}}}{x}+\frac{\alpha ^2}{2}\frac{}{x}\left(\frac{\widehat{v}^j}{x^{\mathrm{}}}\frac{\widehat{v}^j}{x^{\mathrm{}}}\right).$$
(B.15)
The expression for the smoothed velocity can be inverted to give $`𝐯`$ in the form

$$𝐯=\widehat{𝐯}\alpha ^2^2\widehat{𝐯},$$
(B.16)
when the density is constant. Substituting this expression into (B.15) we get

$$\alpha ^2\left(^2\widehat{v}^{\mathrm{}}\frac{\widehat{v}^{\mathrm{}}}{x}+\frac{1}{2}\frac{}{x}\left(\frac{\widehat{v}^j}{x^{\mathrm{}}}\frac{\widehat{v}^j}{x^{\mathrm{}}}\right)\right).$$
(B.17)
Expression (B.13) the continuum limit of the SPH term and (B.17), the continuum term in Holm’s equation, are therefore identical.


format: text

text:
# 1 Introduction

## 1 Introduction

In quenched chiral perturbation theory , it has been asserted that some quantities in QCD possess logarithmic dependence on the bare quark masses as the latter approach zero, similar to the unquenched case. In particular, for the pseudoscalar meson mass and the condensate, there are extra chiral logarithms due to the $`\eta ^{}`$ loop void of topological screening in the quenched approximation. Explicitly, the pion mass to one-loop order can be written as

$`m_\pi ^2=Cm_q\{1\delta [\text{ln}(Cm_q/\mathrm{\Lambda }_\chi ^2)+1]\}+Bm_q^2,`$ (1)
where $`m_q`$ denotes the bare ( $`u`$ and $`d`$ ) quark mass, $`\mathrm{\Lambda }_\chi `$ the chiral cutoff which can be taken to be $`2\sqrt{2}\pi f_\pi `$ ( $`f_\pi 132`$ MeV ), $`C`$ and $`B`$ are parameters, and $`\delta `$ the coefficient of the quenched chiral logarithm.

Theoretically, $`\delta `$ can be estimated to be

$`\delta ={\displaystyle \frac{m_\eta ^{}^2}{8\pi ^2f_\pi ^2N_f}},`$ (2)
where $`m_\eta ^{}`$ denotes the $`\eta ^{}`$ mass in quenched chiral perturbation theory, $`f_\pi `$ the pion decay constant, and $`N_f`$ the number of light quark flavors. For $`m_\eta ^{}=900`$ MeV, $`f_\pi =132`$ MeV, and $`N_f=3`$, (2) gives $`\delta 0.2`$. Evidently, if one can extract $`\delta `$ from the data of $`m_\pi ^2`$ in lattice QCD, then $`m_\eta ^{}`$ in quenched chiral perturbation theory can be determined by (2).

Presumably, the chiral logarithm in (1) can be observed unambiguously in quenched lattice QCD. However, it is rather difficult to disentangle a logarithm from a power series. So far, there is no compiling evidence to single out (1) among many functional forms which seemingly can fit the data of $`m_\pi ^2`$ very well ( e.g., $`m_\pi ^2=A+Cm_q+Bm_q^2`$ ).

For lattice QCD with Wilson-Dirac quark, the chiral symmetry is explicitly broken and the quark mass is additively renormalized, thus one does not expect that quenched chiral logarithm can be identified unambiguously at finite lattice spacing. In fact, if one tries to fit (1) to the data of $`m_\pi ^2`$ with Wilson-Dirac quark, one would obtain the coefficient of quenched chiral logarithm, $`\delta 0.060.12`$ , quite smaller than its theoretical expectation $`0.2`$. For lattice QCD with staggered quark, the coefficient of chiral logarithm also turns out to be rather small, $`\delta 0.06`$ , similar to the case of Wilson-Dirac quark. These discrepancies indicate that quenched chiral logarithm may not be properly reproduced in lattice QCD with these two lattice fermion schemes, at finite lattice spacing.

With the realization of exact chiral symmetry on the lattice, one expects that quenched chiral logarithm can be reproduced in lattice QCD with overlap Dirac quarks . In a recent study , the coefficient of quenched chiral logarithm $`\delta `$ is estimated to be in the range $`0.150.4`$ for the chiral cutoff $`\mathrm{\Lambda }_\chi `$ in the range $`0.61.0`$ GeV, by fitting (1) to their data of $`m_\pi ^2`$. At $`\mathrm{\Lambda }_\chi a=2\sqrt{2}\pi f_\pi a=1.17`$ GeV, their data seems to give $`\delta 0.41(9)`$, which is too large comparing with the theoretical estimate, $`\delta 0.180.20`$. Thus, it is not clear whether quenched chiral logarithm has been unambiguously identified.

Besides from the data of $`m_\pi ^2`$, one can also obtain $`\delta `$ (2) by extracting $`m_\eta ^{}`$ from the propagator of the disconnected hairpin diagram. However, to compute the propagator of hairpin is very tedious. Fortunately, due to the chiral symmetry of $`D_c`$ (6) in the quark propagator $`(D_c+m_q)^1`$, only the zero modes of $`D_c`$ can contribute to the hairpin diagram, regardless of the bare quark mass $`m_q`$. Thus one can derive an exact relation (50) between the $`\eta ^{}`$ mass in quenched chiral perturbation theory and the index susceptibility of any Ginsparg-Wilson lattice Dirac operator, without computing the hairpin diagram explicitly. Explicitly, it reads as

$`m_\eta ^{}^2={\displaystyle \frac{4N_f}{f_\pi ^2}}{\displaystyle \frac{(n_+n_{})^2}{V}}`$ (3)
where $`N_f`$ denotes the number of light quark flavors, $`V`$ the space-time volume, and $`(n_+n_+)^2/V`$ the index susceptibility of any Ginsparg-Wilson lattice Dirac operator in the quenched approximation. Then (2) and (3) together gives

$`\delta ={\displaystyle \frac{1}{2\pi ^2(f_\pi a)^4}}{\displaystyle \frac{(n_+n_{})^2}{N_s}}`$ (4)
where $`N_s`$ is the total number of sites of the lattice. A salient feature of (4) is that $`\delta `$ can be determined at finite lattice spacing $`a`$, by measuring the index ( susceptibility ) of the overlap Dirac operator, and with $`f_\pi a`$ extracted from the pion propagator. Further, we suspect that the value of $`\delta `$ (4) is scale invariant at least for a range of lattice spacings including the continuum limit $`a0`$.

Now it is clear that, in order to confirm the presence of quenched chiral logarithm in lattice QCD, one needs to check whether the coefficient $`\delta `$ determined from the data of $`m_\pi ^2`$ agrees with that (4) obtained from the index susceptibility. Further, they should be close to the theoretical expectation $`\delta 0.2`$. This is a requirement for the consistency of the theory, since the quenched chiral logarithm in $`m_\pi ^2`$ (1) is due to the $`\eta ^{}`$ loop coupling to the pion propagator through the mass term ( in the chiral lagrangian ), thus $`\delta `$ ( $`m_\eta ^{}`$ ) must be the same in both cases. Otherwise, a consistent picture of the theory has not been established, and one may infer that something must have gone wrong, either in the formulation or in the implementation of the theory. We regard this consistency requirement as a basic criterion for lattice QCD ( with any fermion scheme ) to realize QCD chiral dynamics in continuum.

In this paper, we examine whether lattice QCD with overlap Dirac quark can satisfy this basic criterion.

## 2 Pion mass and decay constant

In this section, we set up our notations for computing quenched quark propagator as well as meson propagators in lattice QCD with exact chiral symmetry. Then we describe the practical implementation for overlap Dirac quark in our computations. Finally, we present our results of $`m_\pi ^2a^2`$ and $`f_\pi a`$ for 12 bare quark masses, and determine the coefficient of quenched chiral logarithm, $`\delta =0.203\pm 0.014`$.

### 2.1 Some general formulas

First, we recall some general formulas in lattice QCD with exact chiral symmetry. For any massless lattice Dirac operator $`D`$ satisfying the Ginsparg-Wilson relation

$`D\gamma _5+\gamma _5D=2raD\gamma _5D,`$ (5)
it can be written in terms of a chirally symmetric Dirac operator $`D_c`$ ,

$`D=D_c(1+raD_c)^1.`$ (6)
Then the bare quark mass is naturally added to the $`D_c`$ in the numerator of (6),

$`D(m_q)=(D_c+m_q)(1+raD_c)^1.`$ (7)
Thus the quark propagator is

$`(D_c+m_q)^1=(1rm_qa)^1[D^1(m_q)ra],`$ (8)
which is chirally symmetric in the massless limit. This exact chiral symmetry is the crucial feature one would like to preserve for any quark coupling to physical hadrons .

Then the pion propagator can be written as

$`M(\stackrel{}{x},t;\stackrel{}{0},0)=\text{tr}\{\gamma _5(D_c+m_q)^1(0,x)\gamma _5(D_c+m_q)^1(x,0)\},`$ (9)
where the trace tr runs over Dirac and color space. Note that any one of the quark propagators in (9) can form the $`\eta ^{}`$ loop which leads to the quenched chiral logarithm in the pion mass as well as the condensate.

With the pion propagators, one can compute its time correlation function

$`G(t)={\displaystyle \underset{\stackrel{}{x}}{}}M(\stackrel{}{x},t;\stackrel{}{0},0),`$ (10)
which can then be fitted by the usual formula

$`G_\pi (t)={\displaystyle \frac{Z}{2m_\pi a}}[e^{m_\pi at}+e^{m_\pi a(Tt)}],`$ (11)
to extract pion mass $`m_\pi a`$, and pion decay constant,

$`f_\pi a=2m_qa{\displaystyle \frac{\sqrt{Z}}{m_\pi ^2a^2}}.`$ (12)
Besides the exact chiral symmetry, we also require that $`iD_c`$ is hermitian, similar to the Dirac operator $`i\gamma _\mu (_\mu +iA_\mu )`$ in continuum. Then we have $`D_c^{}=\gamma _5D_c\gamma _5`$. As a consequence of hermiticity and chiral symmetry, both quark propagators in (9) can be expressed in terms of that propagating from $`(\stackrel{}{0},0)`$ to $`(\stackrel{}{x},t)`$, i.e.,

$`M(\stackrel{}{x},t;\stackrel{}{0},0)`$ $`=`$ $`\text{tr}\{[(D_c+m_q)^1(0,x)]^{}(D_c+m_q)^1(x,0)\}`$ (13)
$`=`$ $`{\displaystyle \underset{a,b=1}{\overset{3}{}}}{\displaystyle \underset{\alpha ,\beta =1}{\overset{4}{}}}\left[(D_c+m_q)_{}^{1}{}_{a\alpha }{}^{b\beta }(x,0)\right]^{}(D_c+m_q)_{}^{1}{}_{a\alpha }{}^{b\beta }(x,0)`$
where $`a,b`$ denote color indices, and $`\alpha ,\beta `$ Dirac indices. In other words, to obtain the time correlation function (10) of pion, one only needs to compute 12 columns of the matrix $`(D_c+m_q)^1`$.

### 2.2 Practical implementation for overlap Dirac quark

The massless overlap Dirac operator  reads as

$`D=m_0a^1\left(\text{1I}+\gamma _5{\displaystyle \frac{H_w}{\sqrt{H_w^2}}}\right)`$ (14)
where $`H_w`$ denotes the hermitian Wilson-Dirac operator with a negative parameter $`m_0`$,

$`H_w=\gamma _5D_w=\gamma _5(m_0+\gamma _\mu t_\mu +W),`$ (15)
$`\gamma _\mu t_\mu `$ the naive fermion operator, and $`W`$ the Wilson term. Then $`D`$ (14) satisfies the Ginsparg-Wilson relation (5) with $`r=1/(2m_0)`$. Note that the value of $`m_0`$ has to be chosen properly such that the overlap Dirac operator can capture the topology of the gauge configuration, as well as behaving like a massless Dirac fermion even in a topologically trivial gauge background. In other words, on a finite lattice, the proper range of $`m_0`$ for any gauge configuration is smaller than that $`(0<m_0<2)`$ in the free fermion limit. However, one usually does not need to finely tune $`m_0`$ except for rough gauge configurations at strong couplings . In this paper, we fix $`m_0=1.3`$ for our computations.

Now the difficult task is to compute the quark propagator with exact chiral symmetry (8), which amounts to computing

$`D^1(m_q)=D^{}(m_q)\{D(m_q)D^{}(m_q)\}^1,`$
where the positive hermitian operator $`D(m_q)D^{}(m_q)`$ can be written as

$`D(m_q)D^{}(m_q)=\left[m_q^2+\left(m_0^2{\displaystyle \frac{m_q^2}{4}}\right)\left(2+\gamma _5{\displaystyle \frac{H_w}{\sqrt{H_w^2}}}+{\displaystyle \frac{H_w}{\sqrt{H_w^2}}}\gamma _5\right)\right].`$ (16)
Then the quark propagators from the reference point $`(\stackrel{}{0},0)`$ to all sites $`(\stackrel{}{x},t)`$ can be obtained by solving 12 column vectors $`Y_{c\alpha }(\stackrel{}{x},t)`$ ( where $`c`$ and $`\alpha `$ are color and Dirac indices of the quark field at $`(\stackrel{}{0},0)`$ ) in the following linear system for 12 point source vectors on the r.h.s.,

$`D(m_q)D^{}(m_q)Y_{c\alpha }(\stackrel{}{x},t)=\text{1I}_{c\alpha },`$ (17)
where $`D(m_q)D^{}(m_q)`$ is expressed by (16). Due to the rather small size of the physical memory in the present generation of computers, the viable ways to solve the linear system (17) are iterative methods, in which the conjugate gradient algorithm is the most optimal for positive definite matrices. Starting from an inital guess of $`Y^{(0)}`$, we multiply $`DD^{}`$ to $`Y^{(0)}`$, and update $`Y^{(0)}`$ to $`Y^{(1)}`$ according to the conjugate gradient algorithm, then iterate this process until $`Y^{(n)}`$ converges to the solution with desired accuracy $`ϵ`$, i.e.,

$`D(m_q)D^{}(m_q)Y_{c\alpha }^{(n)}(\stackrel{}{x},t)\text{1I}_{c\alpha }<ϵ.`$ (18)
In this paper, we set $`ϵ=10^{11}`$.

Now the updating process involves the multiplication of the matrix

$`{\displaystyle \frac{1}{\sqrt{H_w^2}}}`$ (19)
to the vector $`Y^{(i)}`$. However, (19) does not have analytic closed form. This is the major challenge for lattice QCD with overlap Dirac quark. A way to proceed is to express (19) in terms of a rational approximation

$`{\displaystyle \frac{1}{\sqrt{H_w^2}}}{\displaystyle \underset{l=1}{\overset{k}{}}}{\displaystyle \frac{b_l}{H_w^2+d_l}},`$ (20)
where $`b_l`$ and $`d_l`$ are some positive definite coefficients. Then the matrix-vector multiplication

$`{\displaystyle \frac{1}{\sqrt{H_w^2}}}Y^{(i)}{\displaystyle \underset{l=1}{\overset{k}{}}}{\displaystyle \frac{b_l}{H_w^2+d_l}}Y^{(i)}={\displaystyle \underset{l=1}{\overset{k}{}}}b_lZ_l`$
can be evaluated by invoking another conjugate gradient process to the linear systems

$`(H_w^2+d_l)Z_l=Y^{(i)},l=1,\mathrm{},k.`$ (21)
Instead of solving each $`Z_l`$ individually, one can use multi-shift CG algorithm , and obtain all $`Z_l`$ altogether, with only a small fraction of the total time what one had computed each $`Z_l`$ separately. Evidently, one can also apply multi-shift CG algorithm to (17) to obtain several quark propagators with different bare quark masses.

Now the computation of overlap Dirac quark propagator involves two nested conjugate gradient loops : the so-called inner CG loop (21), and the outer CG loop (17). The inner CG loop is the price what one pays for preserving the exact chiral symmetry at finite lattice spacing.

The computational cost can be saved more than $`50\%`$ by observing that $`D(m_q)D^{}(m_q)`$ (16) commutes with $`\gamma _5`$. Thus one can choose chiral sources on the r.h.s. of (17),

$`(\text{1I}_{c\alpha })^{b\beta }(\stackrel{}{x},t)=\delta _{cb}\delta _{\alpha \beta }p(\stackrel{}{x},t)`$
where $`p(\stackrel{}{x},t)=\delta _{\stackrel{}{x},\stackrel{}{0}}\delta _{t,0}`$ for a point source. Then we have

$`\gamma _5\text{1I}_{c\alpha }`$ $`=`$ $`+\text{1I}_{c\alpha },\alpha =1,2,`$
$`\gamma _5\text{1I}_{c\alpha }`$ $`=`$ $`\text{1I}_{c\alpha },\alpha =3,4,`$
and

$`\gamma _5Y_{c\alpha }`$ $`=`$ $`+Y_{c\alpha },\alpha =1,2,`$
$`\gamma _5Y_{c\alpha }`$ $`=`$ $`Y_{c\alpha },\alpha =3,4.`$
Now one of the two matrix-vector multiplications in the outer CG loop can be eliminated, since the l.h.s. of (17) becomes

$`D(m_q)D^{}(m_q)Y_{c\alpha }=\left\{m_q^2+\left(2m_0^2{\displaystyle \frac{m_q^2}{2}}\right)\left[1+{\displaystyle \frac{(\gamma _5\pm 1)}{2}}H_w{\displaystyle \frac{1}{\sqrt{H_w^2}}}\right]\right\}Y_{c\alpha }`$ (22)
where ‘+’ for $`\alpha =1,2`$, and ‘-’ for $`\alpha =3,4`$. Further, due to the projector $`(\gamma _5\pm 1)/2`$, half of the matrix-vector multiplication in

$`H_w{\displaystyle \frac{1}{\sqrt{H_w^2}}}Y_{c\alpha }`$
can be saved.

The next crucial question is how to fix the values of coefficients $`b_l`$ and $`d_l`$ in (20) such that the inverse square root function can be approximated optimally. Fortunately, this problem has been solved by Zolotarev  in 1877, using Jacobian elliptic functions. A comparative study of Zolotarev’s optimal rational approximation versus other schemes for the overlap Dirac operator has been reported in Ref. .

For the inverse square root function

$`{\displaystyle \frac{1}{\sqrt{x}}},x[1,x_{max}],`$ (23)
the Zolotarev optimal rational approximation is

$`f(x)=d_0{\displaystyle \frac{\underset{l=1}{\overset{k1}{}}(x+c_{2l})}{_{l=1}^k(x+c_{2l1})}}`$ (24)
where

$`c_l={\displaystyle \frac{\text{sn}^2(lK/2k;\kappa )}{1\text{sn}^2(lK/2k;\kappa )}},\kappa =\sqrt{11/x_{max}},`$ (25)
the Jacobian elliptic function $`\text{sn}(u,\kappa )=\eta `$ is defined by the elliptic integral

$`u(\eta )={\displaystyle _0^\eta }{\displaystyle \frac{dt}{\sqrt{(1t^2)(1\kappa ^2t^2)}}},`$ (26)
$`K=u(1)`$ is the complete elliptic integral, and $`d_0`$ is uniquely determined by the condition

$`\text{max}[1\sqrt{x}f(x)]|_{x[1,x_{max}]}=\text{min}[1\sqrt{x}f(x)]|_{x[1,x_{max}]}`$
For our purpose, it suffices to fix $`d_0`$ such that $`f(1)=1`$.

Since the eigenvalues of $`H_w`$ are bounded,

$`\lambda _{min}|\lambda (H_w)|\lambda _{max},`$ (27)
we can rescale $`H_w`$,

$`h_wH_w/\lambda _{min},`$ (28)
such that the eigenvalues of $`h_w^2`$ fall inside the interval $`[1,(\lambda _{max}/\lambda _{min})^2]`$.

After obtaining the partial fraction of (24), we have the optimal rational approximation to the inverse square root of $`H_w^2`$

$`{\displaystyle \frac{1}{\sqrt{H_w^2}}}{\displaystyle \frac{1}{\lambda _{min}}}{\displaystyle \underset{l=1}{\overset{k}{}}}{\displaystyle \frac{b_l}{h_w^2+d_l}}`$ (29)
where

$`d_l`$ $`=`$ $`c_{2l1},`$ (30)
$`b_l`$ $`=`$ $`d_0{\displaystyle \frac{\underset{i=1}{\overset{k1}{}}(c_{2i}c_{2l1})}{_{i=1,il}^k(c_{2i1}c_{2l1})}},`$ (31)
and the parameter $`\kappa `$ in the Jacobian elliptic function $`\text{sn}(lK/2k;\kappa )`$ is

$`\kappa =\sqrt{1(\lambda _{min}/\lambda _{max})^2}.`$ (32)
Note that, for a given order $`k`$, the smaller the interval $`[1,(\lambda _{max}/\lambda _{min})^2]`$, the more accurate the rational approximation is. Thus, it is advantageous to narrow the interval $`[\lambda _{min},\lambda _{max}]`$ by projecting out the eigenmodes at both ends. Since the spectrum of $`H_w^2`$ is very dense near the upper bound, we only need to project out one or two eigenmodes near $`\lambda _{max}`$, just to get the precise value of $`\lambda _{max}`$. On the other hand, we usually project out 20 or more low-lying eigenmodes of $`H_w^2`$, such that $`(\lambda _{max}/\lambda _{min})^2`$ is always less than 2500. Then we project the inner CG loop to the complement of the vector space spanned by these eigenmodes. This not only gives a better rational approximation for the inverse square root of $`H_w^2`$, but also reduces the number of iterations of the inner CG loop, which is proportional to the condition number $`(\lambda _{max}/\lambda _{min})^2`$.

We use Arnoldi algorithm to compute a selected subset of eigenmodes of $`H_w`$, which correspond to high and low lying ones of $`H_w^2`$. Denoting these eigenmodes by

$`H_wu_j=\lambda _ju_j,j=1,\mathrm{},n,`$ (33)
we project the inner CG loop (21) to the complement of the vector space spanned by these eigenmodes ( i.e., multiplying both sides of (21) by the projector $`P=1_{j=1}^nu_ju_j^{}`$, and use $`PH_w=H_wP`$ ),

$`(h_w^2+d_l)\overline{Z}_l=\overline{Y}^{(i)}(1{\displaystyle \underset{j=1}{\overset{n}{}}}u_ju_j^{})Y^{(i)}.`$ (34)
Note that $`\lambda _{min}`$ and $`\lambda _{max}`$ in Eqs. (27)-(32) now refer to upper and lower bounds of $`|\lambda (\overline{H}_w)|`$, where

$`\overline{H}_w=PH_w=H_w{\displaystyle \underset{j=1}{\overset{n}{}}}\lambda _ju_ju_j^{}.`$ (35)
Then the matrix-vector multiplication at each step of outer CG loop can be expressed in terms of the projected eigenmodes (33) plus the solution obtained from the inner CG loop (34) in the complementary vector space, i.e.,

$`H_w{\displaystyle \frac{1}{\sqrt{H_w^2}}}Y^{(i)}{\displaystyle \frac{1}{\lambda _{min}}}\overline{H}_w{\displaystyle \underset{l=1}{\overset{k}{}}}b_l\overline{Z}_l+{\displaystyle \underset{j=1}{\overset{n}{}}}{\displaystyle \frac{\lambda _j}{\sqrt{\lambda _j^2}}}u_ju_j^{}Y^{(i)}S^{(i)}.`$ (36)
The accuracy of $`S^{(i)}`$ can be measured by the deviation

$`\sigma _i={\displaystyle \frac{|S_{}^{(i)}{}_{}{}^{}S^{(i)}Y_{}^{(i)}{}_{}{}^{}Y^{(i)}|}{Y_{}^{(i)}{}_{}{}^{}Y^{(i)}}},`$ (37)
which would be zero if (36) is exact.

In this paper, we project out 20 lowest-lying and 4 largest eigenmodes of $`H_w^2`$, and use 20 terms ( $`k=20`$ ) in the Zolotarev optimal rational approximation, and set the convergence criterion for inner CG loop to $`10^{12}`$, then $`\sigma _i`$ is always less than $`10^{11}`$ for any iteration of the outer CG loop.

### 2.3 Results

With the Wilson $`SU(3)`$ gauge action at $`\beta =5.8`$, we generate 100 gauge configurations on the $`8^3\times 24`$ lattice, using Creutz-Cabibbo-Marinari heat bath algorithm . For each configuration, we compute quenched quark propagators for 12 bare quark masses. Then the pion propagator (13) and its time correlation functions (10) are obtained, and the latter are fitted by (11) to yield the pion mass $`m_\pi a`$ and decay constant $`f_\pi a`$.

In Figs. 1 and 2, we plot the pion mass square $`(m_\pi a)^2`$ and decay constant $`f_\pi a`$ versus the bare quark mass $`m_qa`$ respectively. The data of $`f_\pi a`$ can be fitted by

$`f_\pi a=0.0984(3)+0.1635(15)(m_qa)`$ (38)
with $`\chi ^2/d.o.f.=0.03`$. Now we take the value $`f_\pi a=0.0984(3)`$ at $`m_qa=0`$ equal to the experimental value of pion decay constant $`f_\pi =132`$ MeV times the lattice spacing $`a`$,

$`f_\pi a=0.0984(3)=132\text{ MeV }\times a,`$ (39)
then this gives an estimate of the lattice spacing

$`a=0.147(1)\text{ fm }.`$ (40)
Thus the lattice size is $`(1.2\text{ fm })^3\times 3.5\text{ fm }`$. Since the smallest pion mass is $`418`$ MeV, the lattice size is $`(2.5)^3\times 7.5`$, in units of the Compton wavelength ( $`0.47`$ fm ) of the lowest-mass pion.

Fixing the cutoff $`\mathrm{\Lambda }_\chi =2\sqrt{2}\pi f_\pi `$ as in quenched chiral perturbation theory ( i.e., $`\mathrm{\Lambda }_\chi a=2\sqrt{2}\pi f_\pi a=0.8745`$ ), we fit (1) to our data of $`(m_\pi a)^2`$, and obtain

$`\delta `$ $`=`$ $`0.2034(140),`$ (41)
$`Ca`$ $`=`$ $`1.1932(182),`$ (42)
$`B`$ $`=`$ $`1.1518(556),`$ (43)
with $`\chi ^2/d.o.f.=0.03`$.

Even though one may not easily detect the quenched chiral logarithm in Fig. 1, it can be unveiled by plotting $`(m_\pi a)^2/(m_qa)`$ versus $`m_qa`$, as shown in Fig. 3. Further, we can subtract the quadratic term $`Bm_q^2`$ \[ with the value of $`B`$ given in (43) \], and check its dependence on $`\mathrm{ln}(m_q)`$ explicitly. In Fig 4, we plot $`(m_\pi a)^2/(m_qa)B(m_qa)`$ versus $`\mathrm{log}(m_qa)`$. The presence of quenched chiral logarithm is evident.

In passing, we note that a finite lattice must impose a lower bound for the pion mass, as well as the corresponding one for the bare quark mass. If one decreases the bare quark mass beyond its lower bound, then the resulting pion mass would be flattening ( or increasing ) rather than decreasing. This is essentially due to the finite volume effects of the zero modes of the overlap Dirac operator in the pion propagator, which is proportional to $`|n_+n_{}|/(m_q^2a^2N_s)`$ as $`m_qa0`$. Thus, only in the infinite volume limit ( $`N_s\mathrm{}`$ ), one can obtain zero pion mass with zero bare quark mass. Nevertheless, for a finite lattice, the finite volume effects of the zero modes can be suppressed if the bare quark mass $`m_qa\sqrt{|n_+n_{}|/N_s}`$. For a lattice of size $`8^3\times 24`$ at $`\beta =5.8`$, we find that $`m_qa0.06`$ is sufficient to suppress the finite volume effects of the zero modes. This can be verified by comparing the results (41)-(43) to those obtained in larger volumes, e.g., $`12^3\times 24`$ at $`\beta =5.8`$ ( which is about $`3.4`$ times of the present volume ). Another plausible check is to see whether the $`\delta `$ extracted from the pion mass agrees with that obtained from the index susceptibility, Eq.(4), since the index susceptibility and the $`f_\pi a`$ ( see Fig. 2 ) both presumably suffer the least from the finite volume effects. The latter check is performed in Section 3, and the agreement between these two $`\delta ^{}s`$ suggests that the finite volume effects may be under control.

## 3 Index susceptibility and the $`\eta ^{}`$ mass

In this section, we derive a relationship between the index susceptibility of any Ginsparg-Wilson lattice Dirac operator and the $`\eta ^{}`$ mass as depicted in quenched chiral perturbation theory. Then we measure the index of overlap Dirac operator for those 100 gauge configurations used for computing quark propagators in Section 2, and obtain the index susceptibility, the $`\eta ^{}`$ mass, and an estimate of $`\delta =0.197\pm 0.027`$, which is in good agreement with the $`\delta `$ determined from the pion mass in Section 2.

### 3.1 A formula due to exact chiral symmetry

The propagator of ( flavor-singlet ) $`\eta ^{}`$ can be written as

$`M_\eta ^{}(x,y)`$ $`=`$ $`\text{tr}\{\gamma _5(D_c+m_q)^1(y,x)\gamma _5(D_c+m_q)^1(x,y)\}`$ (44)
$`\text{tr}\{\gamma _5(D_c+m_q)^1(y,y)\}\text{tr}\{\gamma _5(D_c+m_q)^1(x,x)\}`$
where the second term on r.h.s. is the commonly called hairpin diagram. Note that the first term is just the usual pion propagator (9), while the second term can be regarded as pion propagator with the $`\eta ^{}`$ mass insertion due to the nontrivial interactions of gluons between two seemingly disconnected quark propagators.

In the quenched approximation, the expectation value of the hairpin diagram in momentum space can be written as

$`{\displaystyle \frac{1}{V}}{\displaystyle \underset{x,y}{}}e^{ip(xy)}\text{tr}\{\gamma _5(D_c+m_q)^1(y,y)\}\text{tr}\{\gamma _5(D_c+m_q)^1(x,x)\}`$ (45)
$`=`$ $`\sqrt{Z}{\displaystyle \frac{1}{p^2+m_\pi ^2}}m_\eta ^{}^2{\displaystyle \frac{1}{p^2+m_\pi ^2}}\sqrt{Z}`$
where the brackets $`\mathrm{}`$ denote quenched average over gauge configurations with weight $`e^{𝒜_g}`$ ( $`𝒜_g`$ : pure gauge action ), and

$`\sqrt{Z}={\displaystyle \frac{m_\pi ^2f_\pi }{2m_q}}`$
as defined in (11) and (12). Equation (45) can be regarded as the definition of the $`\eta ^{}`$ mass in quenched chiral perturbation theory.

At zero momentum $`p=0`$, (45) becomes

$`{\displaystyle \frac{1}{V}}{\displaystyle \underset{x,y}{}}\text{tr}\{\gamma _5(D_c+m_q)^1(y,y)\}\text{tr}\{\gamma _5(D_c+m_q)^1(x,x)\}={\displaystyle \frac{f_\pi ^2m_\eta ^{}^2}{4m_q^2}}`$ (46)
Now the expression inside the brackets $`\mathrm{}`$ can be evaluated exactly in terms of the index of $`D_c`$ ( $`D`$ ). Explicitly,

$`{\displaystyle \underset{x}{}}\text{tr}\{\gamma _5(D_c+m_q)^1(x,x)\}=\text{Tr}\{\gamma _5(D_c+m_q)^1\}`$ (47)
$`=`$ $`(1rm_qa)^1{\displaystyle \underset{\alpha }{}}\varphi _\alpha ^{}\gamma _5\varphi _\alpha \{[m_q+(1rm_qa)\lambda _\alpha ]^1ra\}`$
where (7) and (8) have been used, and the trace Tr is evaluated in terms of the eigenvalues $`\{\lambda _\alpha \}`$ and eigenvectors $`\{\varphi _\alpha \}`$ of the Ginsparg-Wilson lattice Dirac operator $`D`$ (6) with $`D_c`$ satisfying $`(iD_c)^{}=iD_c`$. Recalling the properties of the eigensystem of $`D`$ ,

$`\varphi _\alpha ^{}\gamma _5\varphi _\alpha =\{\begin{array}{cc}0\hfill & \text{ for }\lambda _\alpha \lambda _\alpha ^{}\text{ }\hfill \\ \pm 1\hfill & \text{ for }\lambda _\alpha =\lambda _\alpha ^{}\text{ }\hfill \end{array}`$
one immediately sees that only the real eigenmodes $`\lambda =0`$ and $`(ra)^1`$ can contribute to (47). Since the factor $`\{[m_q+(1rm_qa)\lambda _\alpha ]^1ra\}`$ vanishes for $`\lambda _\alpha =(ra)^1`$, thus only the zero modes contribute to (47), i.e.,

$`{\displaystyle \underset{x}{}}\text{tr}\{\gamma _5(D_c+m_q)^1(x,x)\}={\displaystyle \frac{n_+n_{}}{m_q}}.`$ (49)
Substituting (49) into l.h.s. of (46), one obtains

$`N_f{\displaystyle \frac{(n_+n_{})^2}{V}}={\displaystyle \frac{f_\pi ^2m_\eta ^{}^2}{4}},`$ (50)
where the factor $`N_f`$ denotes the number of light quark flavors, which accounts for the number of hairpin diagrams contributing to the $`\eta ^{}`$ mass. Since (50) is an exact result following from the definition of $`m_\eta ^{}`$ in (45), it offers the best way to obtain the $`\eta ^{}`$ mass in quenched chiral perturbation theory. One just measures the index ( susceptibility ) of the Ginsparg-Wilson lattice Dirac operator, without computing the hairpin diagram at all. This is one of the advantages of preserving exact chiral symmetry on the lattice. Note that (50) is independent of the bare quark mass $`m_q`$, thus the $`\eta ^{}`$ can remain massive even in the limit $`m_q0`$. This is the crucial feature of the $`\eta ^{}`$ to distinguish itself from the octet of approximate Goldstone bosons, the $`\pi ,K`$ and $`\eta `$.

Recently, (50) is derived  differently starting from the anomalous flavor-singlet axial Ward identity, which is formally equivalent to replacing the continuum topological charge density $`\rho `$ in the Witten-Veneziano formula

$`m_\eta ^{}^2={\displaystyle \frac{4N_f}{f_\pi ^2}}{\displaystyle d^4x\rho (x)\rho (0)},\rho ={\displaystyle \frac{1}{32\pi ^2}}ϵ_{\mu \nu \lambda \sigma }\text{tr}(F_{\mu \nu }F_{\lambda \sigma }),`$ (51)
with the axial anomaly $`q=\text{tr}[\gamma _5(1raD)]`$ of the overlap Dirac operator, and the integration $`d^4x`$ with the summation over all sites $`_x`$. That is, the integral in (51) is formally transcribed as

$`{\displaystyle d^4x\rho (x)\rho (0)}a^4{\displaystyle \underset{x}{}}q(x)q(0).`$ (52)
Since (51) is derived in the large $`N_c`$ limit ( which is valid to lowest order in $`1/N_c`$ ), a relevant question is what value of $`N_c`$ in Eq. (51) is supposed to be valid. Further, even though the r.h.s. of (52) is well-defined on a finite lattice, some subtleties  could emerge as the 2-point function $`q(x)q(0)`$ ( $`x0`$ ) becomes non-positive in the continuum limit, since $`\rho (x)`$ is odd under time reflections. Thus, a real $`m_\eta ^{}`$ would require a divergent contact term $`\rho ^2(0)`$. How to fix this contact term is the basic problem of Witten-Veneziano formula (51) in the continuum. On the other hand, for a finite lattice, one can use the well-known relation $`_xq(x)=n_+n_{}`$ to turn the r.h.s. of (52) into the index susceptibility,

$`a^4{\displaystyle \underset{x}{}}q(x)q(0)={\displaystyle \frac{1}{V}}{\displaystyle \underset{x}{}}{\displaystyle \underset{y}{}}q(x)q(y)={\displaystyle \frac{1}{V}}(n_+n_{})^2,`$ (53)
which is well-defined in the limit $`V\mathrm{}`$ and $`a0`$. Nevertheless, it is still interesting to see whether $`q(x)q(0)`$ ( $`x0`$ ) would become negative in the continuum limit, together with a divergent contact term.

It should be emphasized that our derivation of (50) is quite different from that of Ref. . We start from Eq. (45) which serves as the definition of the $`\eta ^{}`$ mass in quenched chiral perturbation theory. Then the subsequent steps leading to (50) are exact due to chiral symmetry, without any approximations or assumptions. This is in contrast to the approach used in Ref. , in which the anomalous flavor-singlet axial Ward identity is used, in the limit $`uN_f/N_c0`$. The subtleties arising from $`q(x)q(0)`$ in the continuum limit cannot be resolved unless one transcribes the expression in terms of the index susceptibility of a Ginsparg-Wilson lattice Dirac operator. Finally we note that a derivation of (50) similar to our approach is also presented in Ref.  for the overlap Dirac operator.

### 3.2 Measurement of the index of the overlap

The index of overlap Dirac operator (14) is

$`\text{index}(D)={\displaystyle \underset{x}{}}\text{tr}[\gamma _5(1arD(x,x))]={\displaystyle \frac{1}{2}}(h_{}h_+)`$ (54)
where $`h_+(h_{})`$ is the number of positive ( negative ) eigenvalues of the hermitian Wilson-Dirac operator $`H_w`$ (15). However, one does not need to obtain all eigenvalues of $`H_w`$ in order to know how many of them are positive or negative. The idea is simple. Since $`H_w`$ has equal number of positive and negative eigenvalues for $`m_00`$, then one can just focus at those low-lying ( near zero ) eigenmodes of $`H_w`$, and see whether any of them crosses zero from positive to negative, or vice versa, when $`m_0`$ is scanned from $`0`$ up to the value ( e.g., $`1.30`$ in this paper ) used in the definition of $`D`$. From the net number of crossings, one can obtain the index of $`D`$ (54). This is the spectral flow method used in Refs.  to obtain the index of overlap Dirac operator.

In Fig. 5, we plot the spectral flow of eight lowest-lying ( near zero ) eigenvalues of $`H_w(m_0)`$ in the interval $`0.9m_01.3`$, for one of the gauge configurations. Evidently, the flows are not as smooth as one may have expected. In this case, the net crossings is -6 ( from negative to positive ), so the index of the overlap Dirac operator is -6. In some cases ( see, for example, Fig 6 ), there are some intriguing eigenvalues lying very close to zero, thus for a coarse scan in $`m_0`$, it may not be so easy to determine whether they actually cross zero or not. These ambiguities can only be resolved by tracing them closely at a finer resolution in $`m_0`$.

The distribution of the indices of $`D`$ for these 100 gauge configurations is listed in Table 1. So the index susceptibility is

$`a^4\chi ={\displaystyle \frac{(n_+n_{})^2}{N_s}}=3.67(50)\times 10^4`$ (55)
where $`N_s`$ is the total number of sites. Now if we subtract the mean value $`n_+n_{}=0.19`$ from all indices, then the index susceptibility becomes

$`a^4\chi ={\displaystyle \frac{(n_+n_{})^2n_+n_{}^2}{N_s}}=3.64(50)\times 10^4`$ (56)
Now substitute the index susceptibility (56), the lattice spacing (40), $`f_\pi =132`$ MeV, and $`N_f=3`$, into the exact relation (50), we obtain the $`\eta ^{}`$ mass

$`m_\eta ^{}=(901\pm 64)\text{MeV},`$ (57)
which agrees with the theoretical estimate

$`\sqrt{𝐦_\eta ^{}^2+m_\eta ^22m_K^2}=853\text{ MeV },`$ (58)
with experimental values of meson masses<sup>1</sup><sup>1</sup>1Here we distinguish the physical $`\eta ^{}`$ mass ( $`𝐦_\eta ^{}`$ ), from the $`\eta ^{}`$ mass ( $`m_\eta ^{}`$ ) in queneched chiral perturbation theory.  : $`𝐦_\eta ^{}=958`$ MeV, $`m_\eta =547`$ MeV, and $`m_K=495`$ MeV.

Next we substitute (56) and $`f_\pi a=0.0984(3)`$ into (4), and get

$`\delta =0.197\pm 0.027,`$ (59)
which is in good agreement with the value (41) determined from the pion mass, as well as with the theoretical estimate $`\delta 0.2`$.

## 4 Concluding remarks

In this paper, we have examined the quenched chiral logarithm in lattice QCD with overlap Dirac quark. The coefficient of quenched chiral logarithm ( $`\delta =0.203\pm 0.014`$ ) extracted from the pion mass agrees very well with that ( $`\delta =0.197\pm 0.027`$ ) obtained from the index susceptibility of overlap Dirac operator. Further, they are in good agreement with the theoretical estimate $`\delta 0.2`$ in quenched chiral perturbation theory. This provides strong evidences that lattice QCD with overlap Dirac quark realizes quenched QCD chiral dynamics, as depicted by quenched chiral perturbation theory.

Our results on pion mass and decay constant rely very much on the viability of computing quark propagators to a high accuracy. We find that the accuracy in the implementation of the inverse square root of $`H_w^2`$ is the most crucial step in this two-level conjugate gradient paradigm. The Zolotarev optimal rational approximation together with projection of high and low-lying eigenmodes enables us to control the error of $`\text{sign}(H_w)Y`$ always less than $`10^{11}`$, at each iteration of the outer CG loop. The details of our implementation are described in Section 2.

Equation (50) is an exact result following from the definition of the $`\eta ^{}`$ mass in (45). It provides the best way to obtain the $`\eta ^{}`$ mass in quenched chiral perturbation theory, as well as to determine the coefficient of quenched chiral logarithm in terms of the index susceptibility of any topologically-proper Ginsparg-Wilson lattice Dirac operator. Equation (4) may even suggest that $`\delta `$ is scale invariant for a range of lattice spacings including the continuum limit $`a0`$. This also explains why the $`\delta `$ value we obtained at $`a=0.147`$ fm is so close to the theoretical estimate in the continuum.

At this point, it may be interesting to recall an example of lattice Dirac operator with exact chiral symmetry . Unlike the overlap Dirac operator, it does not have topological zero modes for any nontrivial gauge backgrounds. However, it reproduces correct axial anomaly in the continuum limit , at least for the trivial gauge sector. Since its index is zero, its index susceptibility must be zero, and it follows that the $`\eta ^{}`$ mass is zero, and the quenched chiral logarithm is absent in the pion mass as well as other related quantities. Nevertheless, it does not necessarily imply that this lattice Dirac operator could not be realized in nature. It seems that the real testing ground for this lattice Dirac operator is lattice QCD with dynamical quarks, in which all quenched pathologies are absent.

In passing, we briefly outline the essential features of our computing system. The platform is a home-made Linux PC cluster with 18 nodes, built with off-the-shelf components. Each node consists of one Pentium 4 processor ( 1.6 Ghz ) with one Gbyte of Rambus, one $`40`$ Gbyte hard disk, and a network card. The computational intensive parts ( matrix-vector operations ) of our program are written in SSE2 codes , which can double the speed of our earlier pure Fortran codes compiled by Intel Fortran Compiler 5.0 with maximum optimizations including the SSE2 option. The performance of our system is estimated to be around 20 Gflops.

This work was supported in part by the National Science Council, ROC, under the grant number NSC90-2112-M002-021, and also in part by NCTS.