text:
# Cosmological SPH simulations: The entropy equation

## 1 Introduction

Smoothed particle hydrodynamics (SPH) was introduced by Lucy (1977) and Gingold & Monaghan (1977) as an alternative to grid-based fluid solvers and has developed into a mature and popular simulation technique. The Lagrangian nature of SPH enables it to adjust to the large dynamic range posed by problems such as the formation of galaxies and large-scale structure in a manner that is difficult to match with non-adaptive Eulerian methods. It has been shown that SPH can produce accurate results in many situations (e.g. Steinmetz & Müller, 1993), albeit at the expense of an artificial viscosity which broadens shocks over several smoothing lengths.

The use of SPH to study the formation and evolution of galaxies was pioneered more than a decade ago by a number of researchers (Evrard, 1988, 1990; Hernquist, 1989; Hernquist & Katz, 1989; Navarro & Benz, 1991; Barnes & Hernquist, 1991; Katz & Gunn, 1991; Hiotelis & Voglis, 1991; Thomas & Couchman, 1992; Katz, Hernquist & Weinberg, 1992, among others). Together with subsequent applications, these works demonstrated the reliability of the approach and ultimately influenced our understanding of e.g. quasar absorption line systems, the intracluster medium, and galaxy interactions. While cosmological SPH simulations have also provided a successful instrument for exploring galaxy formation, it has recently become clear that subtle numerical properties of existing algorithms severely influence characteristics of model galaxies, such as their luminosities and colours.

These numerical complications follow once radiative cooling of the gas is included. Cooling is critically important in the conventional paradigm of galaxy formation, because it enables baryons to accumulate in the centres of halos, where they act as reservoir of cold, dense gas for forming stars. Thus, cosmological hydrodynamic simulations must incorporate these dissipative effects to determine the formation sites and stellar masses of galaxies. Ultimately, an appropriate treatment of star formation and associated feedback processes is also required. However, it is clear that reliable cooling rates, uncompromised by numerical effects, are essential if the simulations are to predict accurate galaxy properties.

Unfortunately, calculating the proper amount of gas that can radiatively cool is a difficult problem, especially in hierarchical cosmologies where the most abundant objects are of low mass and hence will invariably be relatively poorly resolved. Naively, it is expected that coarse resolution will tend to inhibit cooling because the radiative cooling rate per unit volume is a sensitive function of local density and high density peaks will be washed out by smoothing. Near the resolution limit of SPH, this effect will suppress cooling entirely, since halos will not harbour sufficient quantities of high density gas if they contain $`N_{\mathrm{sph}}`$ gas particles, where $`N_{\mathrm{sph}}`$ is the number of SPH smoothing neighbours. When the mass resolution is poor, gas particles can also be heated by two-body interactions with heavy dark matter particles (Steinmetz & White, 1997), which further reduces the effective cooling rates.

However, inadequate resolution can also act in the opposite sense and enhance the cooling rate in halos. This can arise from, e.g. ‘in-shock cooling’ (Hutchings & Thomas, 2000). In SPH, shock fronts are broadened and gas passing through them is heated only ‘slowly’ to post-shock temperatures. During this transition, gas can radiate energy, artificially depressing post-shock temperatures and reducing the cooling time of post-shock gas. This effect is particularly severe when the functional dependence of the cooling rate on temperature exhibits a peak between pre- and post-shock temperatures (Martel & Shapiro, 2001). The gas can then cool severely in a broadened shock, in extreme cases converting it incorrectly into an isothermal shock. Additional problems can develop in density discontinuities between the hot phase of a halo and the condensed cold gas that has already cooled. Owing to limited resolution, SPH estimates of the density for hot particles passing close to cold particles will be biased high if they sample the dense cold phase, leading to accelerated cooling (Pearce et al., 1999; Thacker et al., 2000; Croft et al., 2001).

Given these complications, and the highly nonlinear nature of radiative cooling, it is unclear whether or not cooling rates obtained by existing cosmological SPH codes are reliable, in particular with respect to the majority of halos which are resolved with relatively few particles. However, despite the fundamental importance of this issue, only limited attempts have been made to examine the systematic dependence of the results on numerical resolution or technique (among others Thacker et al., 2000). Abadi et al. (2000) found good agreement between SPH predictions and a variant of the self-similar cooling wave solution derived by Bertschinger (1989). While self-similar tests are useful, it is difficult to validate codes in typical cosmological settings from them, because in reality the cooling function is not self-similar, and the objects forming early on contain few particles. The latter is a particularly important restriction, because for sufficiently large particle numbers SPH yields the proper continuum limit.

In this study, we are especially interested in deviations from the correct behaviour when the resolution is coarse. To examine this issue, we compare entropy formulations of SPH (Hernquist, 1993) with variants of the standard algorithm in which the thermal energy equation is integrated. By appealing to a variational principle, we also derive a novel version of SPH which manifestly conserves both the energy and the entropy (under appropriate conditions) even when the smoothing lengths are fully adaptive. In order to investigate the numerical properties of these formulations of SPH in situations where cooling is important, we study collapsing and virializing gas spheres in isolation, and we compute a set of full cosmological simulations of structure formation.

We also explore the behaviour of SPH when feedback scenarios are invoked that deposit large quantities of energy into single particles, as in the method proposed recently by Springel & Hernquist (2001) for modelling starbursts and associated galactic winds and outflows. Benz & Thielemann (1990) have argued that in explosive outbursts, energy must be deposited in a smooth manner for SPH to produce reasonable results, and a similar concern has been raised by Klypin (2001, private communication). We show that singular point explosions do indeed present severe difficulties for some of the standard formulations of SPH, but also demonstrate that reasonable solutions are obtained when the entropy equation is integrated, thereby justifying the numerical validity of feedback schemes that make use of point-like energy injection.

The outline of our paper is as follows. In Section 2, we briefly review several common treatments of SPH. We then derive a new conservative formulation in terms of the entropy equation in Section 3. In Section 4, we start our comparisons by considering strong explosions in a uniform background. In Section 5, we analyse systematic differences between various formulations of SPH in a series of collapse simulations of gas spheres that undergo radiative cooling. We then extend this analysis to small cosmological simulations of structure formation in Section 6. Finally, we conclude and summarise in Section 7.

## 2 Entropy formulation of SPH

In the conventional implementation of SPH (e.g. Monaghan, 1992), the density of each particle is computed according to

$$\rho _i=\underset{j=1}{\overset{N}{}}m_jW(|𝒓_{ij}|,h_i),$$
(1)
where $`𝒓_{ij}𝒓_i𝒓_j`$, and $`W(r,h)`$ is the SPH kernel***We assume that the kernel drops to zero at $`r=h`$.. Commonly, variable smoothing lengths are employed so that the number of neighbours for each particle with $`|𝒓_{ij}|h_i`$ is maintained at a nearly fixed value $`N_{\mathrm{sph}}`$. This adaptive kernel estimation implicitly obeys the continuity equation. The thermal energy and momentum equations can then be integrated in time according to

$$\frac{\mathrm{d}u_i}{\mathrm{d}t}=\frac{_i}{\rho _i}+\frac{1}{2}\underset{j=1}{\overset{N}{}}m_j\left(\frac{P_i}{\rho _i^2}+\frac{P_j}{\rho _j^2}+\mathrm{\Pi }_{ij}\right)𝒗_{ij}_i\overline{W}_{ij},$$
(2)
and

$$\frac{\mathrm{d}𝒗_i}{\mathrm{d}t}=\underset{j=1}{\overset{N}{}}m_j\left(\frac{P_i}{\rho _i^2}+\frac{P_j}{\rho _j^2}+\mathrm{\Pi }_{ij}\right)_i\overline{W}_{ij},$$
(3)
respectively (Monaghan & Gingold, 1983; Monaghan, 1992; Hernquist, 1993; Springel et al., 2001). Here,

$$\overline{W}_{ij}=\frac{1}{2}\left[W(|𝒓_{ij}|,h_i)+W(|𝒓_{ij}|,h_j)\right]$$
(4)
is a symmetrised kernel (Hernquist & Katz, 1989), and $`\mathrm{\Pi }_{\mathrm{ij}}`$ denotes the artificial viscosity, for which we adopt a standard form as in Steinmetz (1996). Note that we have also included an emissivity per unit volume $`_i=(\rho _i,u_i)`$ to describe external sinks or sources of energy due to radiative cooling or heating.

Some authors prefer to employ $`\overline{W}_{ij}=W(|𝒓_{ij}|,[h_i+h_j]/2)`$ in the above equations, so that smoothing lengths are symmetrised rather than the kernels. It is also possible to apply a similar procedure to the estimation of the local density (Hernquist & Katz, 1989). However, of larger importance are the choices made for the symmetrisation of the pressure terms in the equations of motion and in the energy equation.

As a variant of pair-wise symmetrisation with arithmetic means, Hernquist & Katz (1989) suggested the use of geometric means instead, which appeared to provide better integration stability for the thermal energy in their numerical tests. In this case, the SPH equations read

$$\frac{\mathrm{d}u_i}{\mathrm{d}t}=\frac{_i}{\rho _i}+\frac{1}{2}\underset{j=1}{\overset{N}{}}m_j\left(2\frac{\sqrt{P_iP_j}}{\rho _i\rho _j}+\mathrm{\Pi }_{ij}\right)𝒗_{ij}_i\overline{W}_{ij},$$
(5)
and

$$\frac{\mathrm{d}𝒗_i}{\mathrm{d}t}=\underset{j=1}{\overset{N}{}}m_j\left(2\frac{\sqrt{P_iP_j}}{\rho _i\rho _j}+\mathrm{\Pi }_{ij}\right)_i\overline{W}_{ij}.$$
(6)
A number of SPH codes have employed this formulation (among others Hernquist & Katz, 1989; Katz et al., 1996; Davé et al., 1997; Carraro et al., 1998; Springel et al., 2001).

While both of the above formulations are pleasingly symmetric between particles $`i`$ and $`j`$, there is actually no compelling need to distribute the pressure work done by a particle pair equally between them. In fact, if one simply uses the SPH estimate for the local velocity divergence to derive the energy equation, one obtains

$$\frac{\mathrm{d}u_i}{\mathrm{d}t}=\frac{_i}{\rho _i}+\underset{j=1}{\overset{N}{}}m_j\left(\frac{P_i}{\rho _i^2}+\frac{1}{2}\mathrm{\Pi }_{ij}\right)𝒗_{ij}_i\overline{W}_{ij}.$$
(7)
This form of the thermal energy equation has frequently been used in SPH computations instead of a symmetrised form (Evrard, 1988; Rasio & Shapiro, 1991; Navarro & White, 1993; Steinmetz & Müller, 1993; Nelson & Papaloizou, 1993; Couchman et al., 1995; Hultman & Källander, 1997; Thacker et al., 2000). It conserves energy just as well, but as Couchman et al. (1995) pointed out, it produces less scatter in the entropy. Note that this formulation cannot give rise to negative temperatures, as can occur for other formulations of the energy equation in certain situations.

In principle, there are many other forms that the dynamical equations can take in SPH. It is clear, however, that differences between various implementations will in general be small in the limit of relatively smooth flows. However, this need not be true for extreme cases, such as point-explosions and cooling flows, as we discuss in what follows. In our tests, we are especially interested in how well different versions of SPH perform when physical conditions are unusual and numerical resolution is only poor or moderate, as is often the case in cosmological simulations of galaxy formation.

A somewhat more fundamental change in the numerical scheme is to formulate SPH in terms of dynamical equations for the entropy (Lucy, 1977; Benz & Hills, 1987; Hernquist, 1993), rather than the internal energy. We can characterise the specific entropy $`s`$ of a fluid element in terms of an entropic function $`A(s)`$, defined by

$$P=A(s)\rho ^\gamma ,$$
(8)
where $`\gamma `$ is the adiabatic index. Rather than following the evolution of the internal energy, we can integrate

$$\frac{\mathrm{d}A}{\mathrm{d}t}=\frac{\gamma 1}{\rho ^\gamma },$$
(9)
for an inviscid fluid. Note that a convective time derivative is used here so that $`A(s)`$ is conserved for each fluid element in an adiabatic flow ($`=0`$). In this approach, the temperature is inferred from

$$u=\frac{A(s)}{\gamma 1}\rho ^{\gamma 1}.$$
(10)
If shocks occur, the entropic function $`A(s)`$ can vary with time even in the absence of other sources or sinks of entropy. In order to allow for this possibility, an artificial viscosity needs to be introduced into SPH. For example, a suitable SPH discretisation of equation (9) is given by

$$\frac{\mathrm{d}A_i}{\mathrm{d}t}=\frac{\gamma 1}{\rho _i^\gamma }(\rho _i,u_i)+\frac{1}{2}\frac{\gamma 1}{\rho _i^{\gamma 1}}\underset{j=1}{\overset{N}{}}m_j\mathrm{\Pi }_{ij}𝒗_{ij}_i\overline{W}_{ij},$$
(11)
which shows that entropy is only generated by the artificial viscosity in shocks, and by external sources of heat, if they are present. In the entropy formulation of SPH, this expression is integrated in place of equation (2). Note that this approach provides tight control on sources of entropy and, in particular, it is possible to manifestly guarantee that the specific entropy of a particle can only grow in time (assuming for the moment that external sources of entropy are unimportant). This property can be exploited to obtain sharper shock profiles, or to devise shock detection algorithms.

In the continuum limit, employing equation (9) is equivalent to solving the gas dynamics by means of the equation for the internal energy. However, when the number of neighbours defining smoothed estimates in SPH is finite, as is always the case in practice, the two formulations exhibit differences, with one or the other yielding a better approximation to the continuum solution. In particular, Hernquist (1993) has shown that while integrating the internal energy in SPH results in good energy conservation, entropy is not conserved even for purely adiabatic flows. On the other hand, if the entropy is integrated, the total energy is not necessarily conserved. Hernquist (1993) also showed that these errors are primarily caused by the use of variable smoothing and the neglect of relevant terms in the dynamical equations. However, even for fixed smoothing lengths, simultaneous conservation of energy and entropy is not manifest in the above formulations of SPH, but is only guaranteed in the continuum limit; i.e. for a very fine sampling of the fluid.

There is perhaps a tendency to take violations of entropy conservation less seriously than those of total energy, which may explain why most SPH implementations have been made by integrating the thermal energy equation together with the equations of motion, resulting in good energy conservation. However, as we demonstrate below, this choice can lead to significant violations of entropy conservation in certain situations which are important for cosmological simulations of galaxy formation.

## 3 A fully conservative formulation of SPH

In principle, it is possible to construct Lagrangian treatments of hydrodynamics that behave better in their conservation properties, for example based on spatial Voronoi tesselations (Serrano & Espanol, 2001). Alternatively, for SPH, one can explicitly account for terms arising from the variation of the smoothing lengths; the so-called $`h`$ terms. Nelson & Papaloizou (1993, 1994) have shown that conservation of energy and entropy can be improved substantially if such $`h`$ terms are included (see also Serna et al., 1996). However, this approach leads to somewhat cumbersome forms for the dynamical equations and introduces noise into smoothed estimates owing to the dependence of the correction terms on the single most distant neighbour. Consequently, this version of SPH has not found widespread usage in astrophysical applications.

With this in mind, we now derive a new formulation of SPH which employs variable smoothing lengths and which conserves energy and entropy (when appropriate) by construction. To this end, consider the Lagrangian

$$L(𝒒,\dot{𝒒})=\frac{1}{2}\underset{i=1}{\overset{N}{}}m_i\dot{𝒓}_i^2\frac{1}{\gamma 1}\underset{i=1}{\overset{N}{}}m_iA_i\rho _i^{\gamma 1}$$
(12)
in the independent variables $`𝒒=(𝒓_1,\mathrm{},𝒓_N,h_1,\mathrm{},h_N)`$, where the thermal energy acts as the potential generating the motion of SPH particles. The densities $`\rho _i`$ are functions of $`𝒒`$, as defined by equation (1). The quantities $`A_i`$ are treated as constants; i.e. the flow is assumed to be strictly adiabatic for now.

We select smoothing lengths by requiring that a fixed mass is contained within a smoothing volume, viz. $`(4\pi /3)h_i^3\rho _i=M_{\mathrm{sph}}`$, where $`M_{\mathrm{sph}}=\overline{m}N_{\mathrm{sph}}`$ relates the mass $`M_{\mathrm{sph}}`$ to the typical number $`N_{\mathrm{sph}}`$ of smoothing neighbours for an average particle mass $`\overline{m}`$. These equations provide $`N`$ constraints

$$\varphi _i(𝒒)\frac{4\pi }{3}h_i^3\rho _iM_{\mathrm{sph}}=0$$
(13)
on the coordinates of the Lagrangian.

We now obtain the equations of motion from

$$\frac{\mathrm{d}}{\mathrm{d}t}\frac{L}{\dot{q}_i}\frac{L}{q_i}=\underset{j=1}{\overset{N}{}}\lambda _j\frac{\varphi _j}{q_i},$$
(14)
where $`N`$ Lagrange multipliers $`\lambda _i`$ have been introduced. The second half of these $`2N`$ equations gives the Lagrange multipliers as

$$\lambda _i=\frac{3}{4\pi }\frac{m_i}{h_i^3}\frac{P_i}{\rho _i^2}\left[1+\frac{3\rho _i}{h_i}\left(\frac{\rho _i}{h_i}\right)^1\right]^1.$$
(15)
With this result, the first half of the equations then yields

$$m_i\frac{\mathrm{d}𝒗_i}{\mathrm{d}t}=\underset{j=1}{\overset{N}{}}m_j\frac{P_j}{\rho _j^2}\left[1+\frac{h_j}{3\rho _j}\frac{\rho _j}{h_j}\right]^1_i\rho _j.$$
(16)
Using

$$_i\rho _j=m_i_iW_{ij}(h_j)+\delta _{ij}\underset{k=1}{\overset{N}{}}m_k_iW_{ki}(h_i),$$
(17)
we finally obtain the equations of motion

$$\frac{\mathrm{d}𝒗_i}{\mathrm{d}t}=\underset{j=1}{\overset{N}{}}m_j\left[f_i\frac{P_i}{\rho _i^2}_iW_{ij}(h_i)+f_j\frac{P_j}{\rho _j^2}_iW_{ij}(h_j)\right],$$
(18)
where the $`f_i`$ are defined by

$$f_i=\left[1+\frac{h_i}{3\rho _i}\frac{\rho _i}{h_i}\right]^1,$$
(19)
and the abbreviation $`W_{ij}(h)=W(|𝒓_i𝒓_j|,h)`$ has been used.

Note that because the potential (thermal) energy of the Lagrangian depends only on coordinate differences, the pairwise force in equation (18) is automatically anti-symmetric. Total energy, entropy, momentum, and angular momentum are therefore all manifestly conserved, provided that the smoothing lengths are adjusted locally to ensure constant mass resolution as defined by equation (13). Although not explicitly present, the $`h`$ terms are thus consistently included to all orders.

Moreover, the equations of motion (18) are remarkably similar to some prior implementations of SPH. In fact, our equations of motion reduce to those employed by Thomas & Couchman (1992) if one sets $`f_i=f_j=1`$. Hence, the modifications needed to realize the benefits noted above in existing SPH codes are minor, necessitating only a few additional computations. In particular, the quantities $`\frac{\rho _i}{h_i}`$ can be easily computed along with the densities themselves. A slightly more involved alteration is required for the algorithm which updates smoothing lengths, since it is now necessary to ensure that a ‘constant mass’ resides within a smoothing volume rather than a constant number of neighbours. The requisite changes are minimal, however, if an iterative bisection algorithm is used to determine the $`h_i`$ and $`\rho _i`$, as for example in the parallel version of GADGET (Springel et al., 2001).

In order to complete our derivation, we now incorporate an artificial viscosity to allow for the handling of shocks. We invoke a viscous force

$$\frac{\mathrm{d}𝒗_i}{\mathrm{d}t}|_{\mathrm{visc}.}=\underset{j=1}{\overset{N}{}}m_j\mathrm{\Pi }_{ij}_i\overline{W}_{ij},$$
(20)
just as in the standard formulation of SPH, and add it to the acceleration given by equation (18). The resulting dissipation of kinetic energy is exactly balanced by a corresponding increase in thermal energy if the entropy is evolved according to equation (11).

In the following sections, we describe detailed tests of various implementations of SPH. In this context, we distinguish between ‘energy’ formulations of SPH, where the thermal energy equation is integrated, and ‘entropy’ formulations, where the entropic function is integrated. We will consider three variants of the energy approach: a ‘standard’ one described by equations (2) and (3); one where the symmetrisation is done using a geometric mean; and one where the standard form for the equations of motion is combined with the asymmetric form of the energy equation. We contrast these schemes with a ‘standard’ entropy formulation, where equation (2) is replaced by the entropy equation (11), and we finally investigate the performance of our new SPH formulation given by equations (11), and (18) together with (20), which we refer to as the ‘conservative entropy approach’. The different combinations of SPH equations used are summarised in Table 1.

## 4 Point-like energy injection

In the feedback model of Springel & Hernquist (2001), starbursts are accompanied by the release of substantial amounts of thermal energy, injected into individual SPH particles which represent the hot ISM left behind by starbursts. These particles can be sufficiently energetic to entirely escape from galactic halos.

As we now discuss, the standard formulation of SPH is not well-suited for dealing with ‘delta-function’ energy distributions such as this. For definiteness, we consider the propagation of spherical Taylor-Sedov blast waves, which describe the gas dynamics resulting from the point-like injection of energy into a homogeneous medium of negligible pressure. This problem thus serves as a test case for the numerical behaviour of a code when starburst-‘explosions’ are produced. The analytic similarity solution for a strong explosion is well known (e.g. Landau & Lifshitz, 1966). After a time $`t`$, the blast wave propagates a distance

$$R(t)=\beta \left(\frac{Et^2}{\rho }\right)^{1/5},$$
(21)
where the constant $`\beta `$ depends on the adiabatic index $`\gamma `$ ($`\beta =1.15`$ for $`\gamma =5/3`$), $`E`$ is the explosion energy, and $`\rho `$ describes the initial density of the ambient gas. Directly at the spherical shock front, the gas density jumps to a maximum compression of $`\rho ^{}/\rho =(\gamma +1)/(\gamma 1)`$, with most of the mass inside the sphere being swept up into a thin radial shell. Behind the shock, the density rapidly declines and ultimately vanishes towards the explosion centre. Modelling this problem accurately in 3D is a challenge for any hydrodynamical code.

In the following, we describe explosions set up in periodic boxes of unit length per side and unit density, using a Cartesian grid to initialise particle distributions. At time $`t=0`$, we inject an explosion energy $`E=1`$ into a single particle, and compute the subsequent 3D-evolution using the parallel SPH-code GADGET (Springel et al., 2001).

In the ‘standard’ formulation of SPH (integration of the thermal energy equation and arithmetic mean for the pressure terms), a perhaps surprising phenomenon is observed: The direct neighbours of the explosion seed are quickly driven to unphysical negative temperatures. This is a result of the pair-wise symmetrisation of the pressure terms in the energy equation. First, the neighbouring particles of the explosion seed are pushed away. Once they gain velocity, the pressure work done by the particle pairs needs to be supplied by a corresponding decline of thermal energy. Due to the symmetrisation, the (large) rate of decline of thermal energy is equally distributed between the particles of each pair. For a strong explosion, this quickly consumes any thermal energy the neighbouring particles of the explosion seed might have had previously, and drives them towards negative temperatures. Formally, total energy would still be conserved in this case, but our code suppresses unphysical negative temperatures and enforces a minimum value of zero for the thermal energy, thereby leading to an unavoidable, strong violation of energy conservation in this case. Note that this behaviour is not a consequence of time integration errors, but is a result of the initial conditions and the particular set of equations used to evolve them. In our tests, we observe energy violations of up to $`25\%`$, generated entirely in the initial phase. After a short time, however, the energy distribution of the central particles becomes sufficiently ‘smooth’ that a successful explosion still develops. Of course, the initial ‘cooling’ of the neighbouring particles and the error in the explosion energy are substantial artifacts of the solution. These problems can be avoided completely if the asymmetric form of the thermal energy equation is used, in which case energy is conserved accurately for these point explosions.

If a geometric mean is instead used to symmetrise the pressure terms, it is clear that the maximum rate of decline of internal energy in particle pairs that involve the explosion seed becomes smaller. However, for the problem of a single point explosion, this damping leads to a completely unphysical solution: No explosion takes place. This is because the geometric mean (nearly) vanishes if the background is (nearly) pressureless.

However, integrating the entropy equation leads to a solution that is as physically well behaved as the one for the asymmetric form of the thermal energy equation. As before, the explosion seed pushes particles away, but because they are compressed, their temperatures increase, as they should, because their entropy is required to increase monotonically. Unphysical negative temperatures cannot result from the dynamics, regardless of how the energy is initially distributed. However, the total energy is not conserved as well with this approach. We observe a characteristic pattern of fluctuations in the total energy, as shown in Figure 1. There is a maximum deviation of $`4\%`$ occurring in the early phase of the expansion of the blast wave, but there is no long-term secular trend, indicating that the entropy-method conserves energy quite well unless the conditions are extreme. The initial fluctuation is related to the strong variations of the smoothing lengths near the beginning of the simulation. However, if our new conservative formulation of the entropy approach is used, the variations of the smoothing lengths are properly accounted for in the equations of motion and, so, energy is conserved at all times.

In Figure 2, we show an example of the time evolution of an explosion set up in a $`32^3`$ grid by injecting the energy into a single particle, and by using the conservative entropy-method for SPH. The explosion is not exactly spherical, but is slightly modulated by effects due to the geometric pattern of the Cartesian grid seen around the exploding particle. Note that the resolution available for following the self-similar evolution of the blast wave effectively grows with time. At a given time, the resolution can be assessed by the ratio of the radius that the blast wave has reached to the initial mean interparticle separation. For example, at time $`t=0.04`$, the effective resolution is 10.17 for the $`32^3`$ grid, with about 4405 particles contained inside the spherical shock front.

In Figure 3, we show radial density profiles for different $`32^3`$ simulations at a time $`t=0.04`$ after the explosion, computed with all five variants of SPH discussed above, but with the initial explosion energy smoothed according to the SPH kernel. Due to the more regular initial conditions, the spherical symmetry of the explosions is now nearly perfect. The results for integrating the entropy, or the thermal energy in its standard or asymmetric form, are similar, but notice that our new conservative entropy-method produces the smallest amount of scatter in the density estimates for particles in the blast wave. However, it is obvious that the variant of SPH where pressure terms are symmetrised with a geometric mean gives a considerably worse match to the analytic Sedov solution than the other methods. The shock wave appears to be delayed in time with this technique.

Finally, in Figure 4 we show radial density profiles for $`16^3`$, $`32^3`$, and $`64^3`$ simulations at a time $`t=0.02`$ after the explosion. In these tests, the explosion energy was deposited into a single particle, and the entropy equation was integrated. This resolution test makes it clear that it is difficult to resolve the internal structure of a thin shock wave in 3D accurately, and that it is challenging to obtain the expected maximum value of the density compression, albeit increasing the number of particles drastically improves the solution. It is also possible to improve resolution by modifying the smoothing algorithm of SPH so that it locally adjusts to properties of the flow (Shapiro et al., 1996; Owen et al., 1998). Note, however, that the propagation of the shock front in time is described quite accurately in all cases, as well as the density structure behind the blast wave. This is true even for the very low resolution of $`16^3`$ and illustrates the remarkable robustness of SPH in three dimensions. The $`32^3`$-result shown also compares very well to the one in Figure 3, where the energy had been deposited in a smoothed fashion.

From these experiments we conclude that formulations of SPH in terms of entropy or in terms of the asymmetric form of the energy equation are able to produce quite reasonable results for the explosive release of energy into a single particle. It has been argued (Benz & Thielemann, 1990; Owen et al., 1998) that SPH can treat a strong explosion only if the energy is initially deposited sufficiently smoothly; i.e. if the explosion energy is distributed according to the smoothing kernel into a group of particles. We think that this argument only really applies when the thermal energy is integrated in its symmetrised form, where neighbouring particles of a single explosion seed can suffer from an unphysical decline in temperature and entropy. While an initial smoothing of the explosion energy can cure this problem, certain starburst algorithms rely on the possibility of locally releasing a large amount of energy in a single particle. Note that using the geometric mean to symmetrise the pressure terms does increase the robustness of the thermal energy method, but yields noticeably worse or even unphysical solutions for the explosion problem.

## 5 Cooling within halos

Radiative cooling within halos is a highly non-linear process that can be severely influenced by numerical resolution. In order to straightforwardly examine such systematic effects, we consider a series of idealised test problems consisting of gas spheres that collapse and virialize under their own self-gravity. We parameterise the models by a ‘virial velocity’ $`V`$. Each gas sphere has a total mass $`M=V^3/(10GH_0)`$ and is initially at rest with a density profile $`\rho (r)1/r`$, a radius $`R=V/(5H_0)`$, and a thermal energy per unit mass of $`u=5.0\times 10^4V^2`$. After being released, these gas spheres collapse to their centres, bounce back, and virialize through a strong outgoing shock. As such, this problem is a version of the ‘Evrard’-collapse (Evrard, 1988), a common test problem for SPH, except that we allow the gas to cool radiatively like a plasma composed of a primordial mix of helium and hydrogen. The actual cooling rates are computed as in Katz et al. (1996). In order to arrive at a clean inner boundary condition to allow for well defined resolution studies, we set the gravitational softening to a fixed value of $`ϵ=1.5\times 10^3V/H_0`$ in all the tests, and we turn cold gas within a distance of less than $`0.25ϵ`$ from the origin and with a density higher than $`7.2\times 10^7\rho _{\mathrm{crit}}`$ into stationary sink particles that only gravitationally interact with the remainder of the gas. This ensures that once it can cool efficiently, the gas at the origin will essentially be fixed at a density $`7.2\times 10^7\rho _{\mathrm{crit}}`$, allowing a meaningful attempt to achieve convergence for the entire density profile as the mass resolution is increased. For the particle numbers employed here, the SPH smoothing lengths of particles in the sink region will be smaller than $`0.5ϵ`$; hence, they are no longer interacting with particles in the actual cooling flow once they are converted into sink particles.

While this test problem is highly idealised, it is nevertheless relevant for simulations of galaxy formation, because in that case the virialization of halos and the subsequent cooling of their gas happens in an analogous, yet dynamically more complicated manner. Note that we deliberately excluded dark matter from this test to highlight gas-dynamical effects without introducing complexities arising from collisionless dynamics. An interesting alternative test would be to examine self-similar cooling wave solutions (Bertschinger, 1989). While having an analytic solution is an advantage, we prefer to use a set-up that is closer to actual cosmological simulations with respect to particle number, to collapse scenario, and to the cooling function. This simplifies the interpretation of numerical effects.

Because the cooling function for primordial gas introduces characteristic temperatures, the evolution depends on the mass scale of the collapsing sphere. We here consider spheres of mass $`M=10^{15}h^1\mathrm{M}_{}`$, where $`h=0.7`$ parameterises the Hubble constant $`H_0=100h^1\mathrm{km}\mathrm{s}^1\mathrm{Mpc}^1`$. These gas spheres thus have the mass of a rich cluster of galaxies, and a virial temperature of $`5\times 10^7\mathrm{K}`$. Note that because these ’clusters’ consist of only gas, their cooling times will be shorter than those of real, dark-matter dominated clusters of comparable mass.

In Figure 5, we show the deposition rate of cold gas in the halo centre as a function of time, computed for mass resolutions ranging from 552 to 92096 particles, and for the five variants of SPH listed in Table 1. The deposition rate is given in units of total mass per unit time, where unit time is $`1/(10H_0)=0.98h^1\mathrm{Gyr}`$. A value of 0.5 thus means an accretion rate such that 50% of the gas cools within $`1.4\mathrm{Gyr}`$ for $`h=0.7`$.

From Figure 5, it is evident that at low resolution there are strong systematic differences between the energy and entropy formulations of SPH: When the thermal energy equation is integrated, cooling rates are substantially overestimated, while with the entropy equation, they are slightly underestimated compared to the ‘converged’ values measured at very high resolution. As the particle number is increased, the various treatments yield results that are in better and better agreement. In itself, this is reassuring, but it is unsettling that the cooling rate can be as dramatically overestimated as is seen in the 552 and 1472 particle runs which integrated the thermal energy equation. After all, many halos in cosmological simulations are resolved with far fewer particles. Therefore, all of these halos are likely subject to strong overcooling. Moreover, in hierarchical scenarios of structure formation, larger systems will inherit this problem if their progenitors suffered from overcooling. Such behaviour may well explain why cosmological SPH simulations predict more gas cooling than is expected based on simple analytic models (Pearce et al., 2001; Benson et al., 2001; Yoshida et al., 2001).

Note from Figure 5 that the low-resolution results for the geometric symmetrisation of the pressure terms are particularly discrepant compared to the high-resolution solutions. Even the simulation with 1472 particles cools all the gas in the initial collapse in this case, and for 4224 and 11536 particles there are still strong ‘shoulders’ in the accretion rate at early times around $`t2`$. Such a shoulder is also present for the other energy-based methods in the 4224 runs, but it is absent for the conservative entropy formulation.

To determine the cause of the difference in the low-resolution behaviour between the entropy and energy formulations of SPH, it is instructive to examine the radial structure of the cooling flows in the simulations. In Figure 6, we show profiles of the density, the temperature, the radial velocity, and the radial mass flux for one of our runs at time $`t=4.0`$, when the cooling flow has been operative for some time and the accretion is almost stationary. It is seen that there is an inflow region in the radial range $`20300h^1\mathrm{kpc}`$ where the gas is moving inwards at nearly constant radial velocity, and the density profile is steeply rising with slope $`2`$. It is important to note that this inflow is essentially adiabatic; the gas is not yet radiating efficiently, hence the temperature is strongly rising as the gas is compressed on its inbound path. Close to the very centre, the gas begins to cool efficiently, leading to an accelerated inflow, and eventually to a rapid drop in the temperature to around $`10^4K`$. It is clear that an underestimate of the compressional heating of the gas in the extended inflow region will artificially boost the overall accretion rate. In fact, we believe it is this effect which is responsible for the overcooling in the thermal energy runs at low resolution, as we now demonstrate.

In the left panel of Figure 7, we show the change $`\mathrm{\Delta }A(s)`$ in the entropy of individual SPH particles between two subsequent simulation outputs of the 1472 particle runs, for both the standard energy and entropy methods. Recall that $`A(s)`$ should change only as a result of radiative losses, or owing to generation of entropy by the artificial viscosity in shocks. We plot $`\mathrm{\Delta }A(s)`$ as a function of $`\rho ^{4/3}u^{1/2}`$, which is the relevant dependence of the entropy loss rate due to cooling in this temperature range. When we compare particles from the thermal energy and the entropy runs in this way, we see that $`A(s)`$ for the particles in the energy run decreases faster than expected based on the cooling rate, where this expectation is given by the locus of particles from the entropy run. This may be expressed differently by saying that the entropy of these particles is violated – it decreases faster than can be accounted for by the external sinks of entropy. It should be noted that in spite of this, the total energy is well conserved for these particles, i.e. the total energy declines at exactly the rate at which energy is radiated.

What we observe in Figure 7 is a discreteness effect arising from insufficient resolution in the cooling flow, where the adiabatic heating of the inflow is not adequately resolved. If the particle number is increased, the resulting false loss of entropy is avoided, as can be seen in the right panel of Figure 7, where we show the equivalent plot for the 11536 particle simulations. We note that the geometric symmetrisation of the thermal energy equation is particularly prone to this effect. At low resolution, an infalling particle that interacts with particles that are already much colder will not be heated sufficiently despite its strong compression, because the geometric mean heavily suppresses adiabatic work done between particle pairs at different temperatures.

In fact, we have found that this relative suppression of pair-wise forces in the geometric symmetrisation technique can be strong enough to render density discontinuities unstable. In order to highlight this effect, we have considered a test problem where two ‘phases’ of gas of equal pressure are brought into contact along a planar intersection, with the density of one of the phases being substantially higher than that of the other. Since we take the gas to be initially at rest in this problem, the resulting density discontinuity should remain stable for all times.

However, as a result of its inherent smoothing, it is clear that SPH will not be able to precisely maintain the sharpness of the initial density discontinuity. Instead, we expect that the gas at the intersection will undergo some small transient motion before it settles into a relaxed state, such that the phase boundary ends up being slightly washed out. We investigate this relaxation process with a three-dimensional realisation of a density discontinuity of strength $`\rho _2/\rho _1=125`$. For definiteness, we consider a periodic box of size $`4\times 1\times 1`$ where one half is filled with a $`2\times 7^3`$ particle grid at a fiducial density $`\rho _1=1`$, and the other half is filled with a $`2\times 35^3`$ grid. All particles are taken to have equal mass, such that the density in the latter half is $`\rho _2=125`$. Temperatures are assigned such that an equal pressure of $`P=1`$ in both phases (using $`\gamma =5/3`$) is obtained.

In Figure 8, we show a comparison of the time evolution of the resulting particle system obtained for three different SPH formulations, using $`N_{\mathrm{sph}}=32`$ smoothing neighbours in each case. In the top three panels, we show the evolution when the new conservative entropy formulation is used. The system shows the expected behaviour. Low-density particles adjacent to the high-density phase ‘see’ some of the dense particles, and are hence estimated slightly high in density. Likewise, the density estimates of particles close to the edge of the high-density phase become biased low, such that the phase transition is sligthly washed out. However, the phase boundary remains reasonably sharp, and its stability is maintained for arbitrarily long times.

In comparison, the behaviour resulting for the standard energy formulation of SPH is quite similar, but the sharpness of the initial density discontinuity is not preserved quite as well. In particular, a larger number of particles from the original low-density phase is affected by higher density estimates attained in the new equilibrium. If the system was subject to density-dependent cooling, this may cause an overestimate of the cooling rates of particles close to the boundary.

This may however be seen as a small effect compared to the relative performance of the geometric symmetrisation technique, which gives a solution that is drastically worse than the ones obtained for the other two methods. In fact, as the lower row of panels in Figure 8 shows, the phase boundary behaves in an unstable way in this case, resulting in an unphysical evolution of the system. Particles from the low-density phase are literally pressed into the high density phase. This happens because the geometric symmetrisation renders the repulsive forces between particle pairs of very different temperature largely ineffective. As a result, each of the $`7\times 7`$ columns of low-density particles orthogonal to the boundary penetrates as a finger into the high-density phase, causing a partial overlap of the fluids, and the gradual decay of the initial density discontinuity. Clearly, if the system was subject to cooling, the hot low-density phase would suffer from a substantial increase of its cooling rate because it would tend to overlap with the dense cold phase, thereby becoming subject to overestimates of density.

## 6 Cosmological results

Based on the results presented above, we expect a significant difference in the fraction of gas that can cool in halos of a given size in full cosmological simulations when SPH is formulated either in terms of the thermal energy equation or the entropy equation. We examine this further using moderate-sized test simulations at two different mass resolutions, containing a total of $`2\times 50^3`$ and $`2\times 100^3`$ dark matter and SPH particles, respectively. We consider a canonical $`\mathrm{\Lambda }`$CDM cosmology with parameters $`\mathrm{\Omega }_0=0.3`$, $`\mathrm{\Omega }_\mathrm{\Lambda }=0.7`$, $`\mathrm{\Omega }_b=0.04`$, $`h=0.67`$, and $`\sigma _8=0.9`$ within a periodic box of length $`11.3h^1\mathrm{Mpc}`$ per side. We include radiative heating and cooling processes for a primordial mix of helium and hydrogen in the way described by Katz et al. (1996). An external UV background field similar to that advocated by Haardt & Madau (1996) photoionizes the gas and reionization takes place at redshift $`z6`$ (for details, see Davé et al., 1999). For the purposes of the present study, we have not renormalized the ionizing background to match e.g. the redshift evolution of the mean opacity of the Lyman-$`\alpha `$ forest, as is commonly done (see e.g. Hernquist et al., 1996). While boxes of this size are too small to be representative of the universe at low redshifts, they suffice to investigate systematic differences arising from different formulations of SPH.

In Figure 9, we show the fraction of baryons residing in different phases as a function of redshift. These phases are defined as in Davé et al. (2001): Gas with temperatures above $`10^7\mathrm{K}`$ is referred to as ‘hot’, gas with temperatures between $`10^5\mathrm{K}`$ and $`10^7\mathrm{K}`$ is ‘warm’, and at lower temperatures gas is either labelled ‘diffuse’ if it is at overdensities $`\rho <1000\overline{\rho }`$, or ‘condensed’ (cold) for $`\rho >1000\overline{\rho }`$. The left panel in Figure 9 shows the result when the thermal energy equation is integrated with geometric mean symmetrisation, while the right panel is for our new conservative entropy method. In each panel, we show the results obtained for the $`2\times 50^3`$ resolution as thick lines, while thin lines give the results for the $`2\times 100^3`$ simulations.

When the two SPH methods are compared, it is obvious that there is a very large difference in the amount of gas that becomes cold and dense. The entropy method predicts almost a factor of 2 less cold gas for these simulations, which is reflected in a corresponding difference in the amounts of gas in the warm phase. The fraction of gas in the diffuse phase is nearly unaffected, however. There is virtually no hot gas in these simulations because the box size is too small to contain massive galaxy clusters. Note that the fraction of cold gas keeps growing for the geometric method by a large amount even at low redshift, while it stays nearly constant for the entropy method for redshifts below $`z2`$. This is likely related to the instability of density discontinuities in the geometric symmetrisation technique, as discussed at the end of Section 5. At low redshift, this allows large halos to keep cooling efficiently because their gas easily ‘overlaps’ with the phase of cold dense gas that has accumulated at the centres of these halos.

As far as convergence is concerned, the total amount of cold gas becomes slightly larger in both SPH variants when better mass resolution is employed. For example, the fraction of cold gas goes up from 20.0% to 23.7% in the entropy method when moving from the $`2\times 50^3`$ to the $`2\times 100^3`$ resolution. Note however that the higher resolution runs are able to probe halos of smaller mass than were previously resolved. The gas that cools in these newly seen objects can possibly account to a large extent for the difference in the total amount of cold gas. A better assessment of convergence is therefore obtained if we compare the content of cold gas in individual halos as a function of their mass. To this end, we find halos in the simulations as follows: We first construct a catalogue of group candidates using a friends-of-friends algorithm with linking length 0.2, applied to the dark matter only. For each dark matter group that contains at least 64 particles, we define a group centre as the position of the particle with the minimum gravitational potential. We then determine a spherical ‘virial radius’ around this centre such that the enclosed ‘virial mass’ of gas and dark matter has an overdensity of 180 with respect to the background. For each of the resulting halos, we also compute its fractional content of gas, $`(M_{\mathrm{gas}}/M_{\mathrm{vir}})\times (\mathrm{\Omega }_0/\mathrm{\Omega }_b)`$, normalized to the universal baryon fraction. Similarly, we define a fraction of cold gas, where cold is here taken to be gas colder than $`5\times 10^4\mathrm{K}`$, and we define a ‘hot’ gas content simply as the difference of the total and cold gas fractions.

In Figure 10, we show these gas fractions as a function of virial mass, comparing results for the $`2\times 50^3`$ and $`2\times 100^3`$ simulations when the geometric symmetrisation is used, or the new entropy formulation is employed. In the two panels on top, we have used logarithmic bins in mass to compute mean gas fractions for halos in each mass bin. This suppresses the scatter, which we separately illustrate for the entropy runs in the lower three panels by plotting indiviudal halos as symbols.

Overall, the agreement between the low- and high-resolution runs is gratifyingly good for the entropy method, as seen in the top right panel of Figure 10. Even at the resolution limit of $`6.15\times 10^{10}h^1\mathrm{M}_{}`$ of our $`2\times 50^3`$ halo catalogue, the fraction of cold gas is well reproduced, apart perhaps from a small bias towards a lower amount of cold gas, which is however expected based on our earlier results for the cooling in isolated gas spheres. Interestingly, at the low mass end, the high-resolution catalogue shows a drop of the total gas fraction below the universal value, while a similar effect is not observed for the low-resolution run. This depletion of the total gas content is presumably caused by the UV background that we included (Gnedin, 2000). At redshift $`z=0`$, this effect appears to be strong enough to start reducing the gas content in objects with virial temperatures below $`T8\times 10^4\mathrm{K}`$.

When the cold gas fractions for the geometric symmetrisation technique are considered, it becomes again clear that this method cools gas too efficiently. Here the low resolution run actually tends to produce more cold gas at a given mass scale than the high resolution run, again in agreement with the trends seen in the collapse simulations of isolated gas spheres discussed earlier. Note that the difference in the cold and total gas fractions between the two different resolutions is also markedly larger than seen for the entropy runs. Interestingly, the total baryon fraction is also pushed beyond the universal value in this method, indicating that the strong cooling also affects the accretion flow into the halos.

The difference in the phase space distributions of the gas between the different SPH techniques may also be shown directly in a density-temperature diagram. In Figure 11, we show where the SPH particles lie in the overdensity-temperature plane. Compared to the simulation that integrates the entropy, the plume of shock heated gas is clearly less populated in the run that evolved the energy equation. This gas has cooled and moved to the locus of dense gas at $`T10^4\mathrm{K}`$. The amount of diffuse, low-density gas is the same in the two simulations. However, the entropy run clearly exhibits much less scatter in the density-temperature relation of the photoionized IGM. Hui & Gnedin (1997) have studied the expected density-temperature relation of the photoionized intergalactic medium. At low overdensities ($`\delta <5`$), the gas is expected to follow a power-law ‘equation of state’

$$T=T_0(1+\delta _\mathrm{b})^\alpha ,$$
(22)
where $`\delta _b`$ is the baryonic overdensity. The amplitude $`T_0`$ and the slope $`\alpha `$ at a given epoch depend on the details of the reionization history, and on cosmology. Because gas with density around or below the mean is not expected to have been shocked, very little scatter is expected around this relation, as long as the ionising radiation field is homogeneous.

We examine the scatter in the equation of state of the IGM in Figure 12. In the top panels we compare the temperature-density relation for the two runs at redshift $`z=4`$, and in the bottom panels at redshift $`z=0`$. While at redshift $`z=4`$ the scatter in the thermal energy equation run is only marginally larger than that in the entropy run, it has grown substantially by redshift $`z=0`$. This happens because the adiabatic cooling of the particles due to the expansion of the universe is perturbed by noise in the SPH estimates. The particles see only a finite number of neighbours that move away from them and generate the adiabatic cooling. While they cool at nearly the correct rate on average, some of them will do so more rapidly or more slowly than average and, as a result, particles tend to diffuse away from the adiabat they are expected to follow. In addition, for a small number of smoothing neighbours and a nearly homogeneous particle distribution, there can be a small bias in the mean of the estimates of the local velocity divergences, resulting in a mean temperature that does not decline exactly in proportion to $`a^2`$, where $`a`$ is the scale factor. In fact, this happens for the thermal energy run, which arrives at an IGM equation of state with a slightly higher value of $`T_0`$. This is a result of a small systematic underestimate of the cooling rate due to the adiabatic expansion of the universe. Ultimately, this is again a manifestation of the violation of entropy conservation in this approach. Note that our code computes the rate of change of thermal energy for each SPH particle in physical coordinates, thus the situation can be compared to a homogeneously expanding polytrope where the smoothing lengths increase with the expansion. As Hernquist (1993) has shown, the entropy will be violated under these circumstances. On the other hand, if the entropy is integrated, the adiabatic cooling is automatically exact. Particles will only move to higher temperatures if they are shocked. The entropy computation thus proceeds essentially as in the ‘semi-analytic’ method outlined by Hui & Gnedin (1997), and should represent the more accurate result.

## 7 Discussion

We have derived a new formulation of SPH that, when appropriate, manifestly conserves both entropy and energy if the smoothing lengths are adjusted to a constant local mass resolution. Unlike previous attempts to include $`h`$-terms, the equations of motion retain a remarkably simple form. In our approach, $`h`$-terms do not have to be evaluated explicitly, mitigating additional sources of noise. This improves on previous schemes where these terms depend on the single particle with the maximum distance among the set of neighbours.

We have investigated the robustness of several formulations of SPH when applied to the problem of cooling in halos under conditions of poor resolution, and when dealing with strong explosions triggered by local energy injection, as it occurs in certain feedback scenarios. As such, these test problems are highly relevant for current attempts to model galaxy formation in cosmological SPH simulations.

We find that the standard formulations of SPH in terms of the thermal energy can lead to substantial overcooling in poorly resolved halos when cooling flows are not adequately resolved. This can lead to insufficient compressional heating in the accretion flow, or expressed differently, to a violation of entropy conservation in the adiabatic region of the flow. In the cooling problem, this leads to an accelerated accretion, an effect that becomes more substantial for halos resolved with few particles, and for temperatures in a range where the cooling function declines with temperature. The effect is particularly severe for a geometric symmetrisation of the pressure terms, which can also cause unstable behaviour of strong density discontinuities.

Point-like energy injection can lead to unphysical negative temperatures in the standard formulation of SPH. When the temperatures of these particles are prevented from becoming negative, a reasonable explosion still takes place once the initial energy has spread over a smoothing volume, but a substantial violation of total energy conservation results. When the pressure terms are symmetrised with a geometric mean, the integration of the thermal energy equation is more stable, but then unphysical behaviour may result for very strong gradients in energy density. A better alternative is to use the asymmetric form of the thermal energy equation, which avoids such artifacts.

On the other hand, when the entropy equation is integrated, point explosions are also treated very well. Even when the initial explosion energy is deposited in a smooth way, the blast waves exhibit less scatter than for an integration of the thermal energy equation. However, the total energy is not conserved well in this case, and shows fluctuations of order several percent, although there are no secular trends of energy violation. This is not surprising. Hernquist (1993) has shown that for standard formulations of SPH, simultaneous conservation of energy and entropy is not manifest. Depending on whether SPH is solved in terms of energy or entropy, one of the conservation laws is only strictly fulfilled in the continuum limit. However, in this paper we have derived a new formulation of SPH that solves this problem and conserves both energy and entropy. In addition, this new method is nearly as easy to implement numerically as the standard approach.

We believe that the results presented in this study clearly favour the entropy formulation of SPH for cosmological simulations of galaxy formation. It provides the best technique among the ones we have tested for modeling the process of point-like energy injection, which is relevant for certain feedback algorithms. It also avoids artificial overcooling in poorly resolved halos and reduces the scatter in the density-temperature relation of the gas in the low density Ly-$`\alpha `$ forest.

## Acknowledgements

We thank Simon White and Simone Marri for instructive discussions and critical comments that were helpful for the work on this paper. We are indebted to David Weinberg for useful suggestions, and to the anonymous referee for providing a highly useful report that helped to improve the paper. This work was supported in part by NSF grants ACI96-19019, AST-9803137, and PHY 9507695.


format: text

text:
# Signals for Non-Commutative QED in 𝑒𝛾 and 𝛾𝛾 Collisions

## Abstract

We study the effects of non-commutative QED (NCQED) in fermion pair production, $`\gamma +\gamma f+\overline{f}`$ and Compton scattering, $`e+\gamma e+\gamma `$. Non-commutative geometries appear naturally in the context of string/M-theory and gives rise to 3- and 4-point photon vertices and to momentum dependent phase factors in QED vertices which will have observable effects in high energy collisions. We consider $`e^+e^{}`$ colliders with energies appropriate to the TeV Linear Collider proposals and the multi-TeV CLIC project operating in $`\gamma \gamma `$ and $`e\gamma `$ modes. Non-commutative scales roughly equal to the center of mass energy of the $`e^+e^{}`$ collider can be probed, with the exact value depending on the model parameters and experimental factors. However, we found that the Compton process is sensitive to $`\mathrm{\Lambda }_{NC}`$ values roughly twice as large as those accessible to the pair production process.

Although string/M-theory is still developing, and the details of its connection to the Standard Model are still unclear, numerous ideas from string/M-theory have affected the phenomenology of particle physics. The latest of these ideas is non-commutative quantum field theory (NCQFT) Douglas ; ncqft . NCQFT arises through the quantization of strings by describing low energy excitations of D-branes in background EM fields. NCQFT generalizes our notion of space-time, replacing the usual, commuting, space-time coordinates with non-commuting space-time operators. Testable differences exist between QFT with commuting space-time coordinates and NCQFT.

At this time, the details of a general NCQFT model to compare to the Standard Model are just emerging ncsm . However, NCQED does exist and can be studied. NCQED modifies QED, with the addition of a non-Lorentz invariant, momentum dependent phase factor to the normal $`ee\gamma `$ vertex, along with the addition of cubic ($`\gamma \gamma \gamma `$) and quartic ($`\gamma \gamma \gamma \gamma `$) coupling, also, with non-Lorentz invariant momentum dependent phase factors. The Feynman rules for NCQED are given in frncqed ; hpr . Although the momentum dependent phase factors and higher dimensional operators in the Lagrangian arise naturally in NCQFT, the modifications, although similar, will in general, take on a different form than those for NCQED. We will see that the modifications of NCQFT to QED can be probed in $`\gamma \gamma f\overline{f}`$ and $`e\gamma e\gamma `$ collisions. For full details of‘¡ our analysis, please see Ref. us .

The essential idea of NCQFT is that in the non-commuting space time the conventional coordinates are represented by operators which no longer commute:

$$[\widehat{X}_\mu ,\widehat{X}_\nu ]=i\theta _{\mu \nu }\frac{i}{\mathrm{\Lambda }_{NC}^2}C_{\mu \nu }$$
(1)
Here we adopt the Hewett-Petriello-Rizzo parametrization hpr  where the overall scale, $`\mathrm{\Lambda }_{NC}`$, characterizes the threshold where non-commutative (NC) effects become relevant and $`C_{\mu \nu }`$ is a real antisymmetric matrix whose dimensionless elements are presumably of order unity. One might expect the scale $`\mathrm{\Lambda }_{NC}`$ to be of order the Planck scale. However, given the possibility of large extra dimensions add ; rs  where gravity becomes strong at scales of order a TeV, it is possible that NC effects could set in at a TeV. We therefore consider the possibility that $`\mathrm{\Lambda }_{NC}`$ may lie not too far above the TeV scale.

The $`C`$ matrix can be parameterized, following the notation of JoA , as

$$C_{\mu \nu }=\left(\begin{array}{cccc}0& C_{01}& C_{02}& C_{03}\\ C_{01}& 0& C_{12}& C_{13}\\ C_{02}& C_{12}& 0& C_{23}\\ C_{03}& C_{13}& C_{23}& 0\end{array}\right)$$
(2)
where $`_i|C_{0i}|^2=1`$. Thus, the $`C_{0i}`$ are related to space-time NC and are defined by the direction of the background E-field. Likewise, the $`C_{ij}`$ are related to the space-space non-commutativeness and are defined by the direction of the background B-field.

NCQED is beginning to attract theoretical and phenomenological interest hpr ; ncpheno ; mathews ; bghh . Hewett, Petriello and Rizzo hpr  have performed a series of phenomenological studies of NCQED at high energy, linear, $`e^+e^{}`$ colliders. They analyzed diphoton production ($`e^++e^{}\gamma +\gamma `$), Bhabha scattering ($`e^++e^{}e^++e^{}`$) and Moller scattering ($`e^{}+e^{}e^{}+e^{}`$). There are striking differences between QED and NCQED for all three processes; most interesting is significant structure in the $`\varphi `$ angular distribution.

Mathews mathews  and Baek, Ghosh, He and Hwang bghh  have also studied NCQED at high energy $`e^+e^{}`$ linear colliders. In the former case Mathews studied high energy Compton scattering while Baek et al., studied fermion pair production in $`\gamma +\gamma e^++e^{}`$. Independently of the aforementioned studies we studied Compton scattering and lepton pair production. Our study uses angular distributions to enhance the sensitivity of measurements to $`\mathrm{\Lambda }_{NC}`$, uses more realistic acceptance cuts, etc.

We consider linear $`e^+e^{}`$ colliders operating at $`\sqrt{s}=0.5`$ and 0.8 TeV appropriate to the TESLA proposal, tesla  $`\sqrt{s}=0.5`$, 1.0 and 1.5 TeV as advocated by the NLC proponents nlc , and $`\sqrt{s}=3.0`$, 5.0 and 8.0 TeV being considered in CLIC studies clic . In order to estimate event rates, we assume an integrated luminosity of $`L=500`$ fb-1 for all cases. We impose acceptance cuts on the final state particles of $`10^o\theta 170^o`$ and $`p__T>10GeV`$. Furthermore, all exclusion limits given below are for unpolarized electron and photon beams; the helicity structure of the NCQED cross section is identical to that in the SM, i.e. the fermion-photon couplings are vector-like, so polarization will not lead to an improvement in the exclusion limits.

In the pair production case, where only space-time NC enters, only one parameter, $`\alpha `$, remains in addition to $`\mathrm{\Lambda }_{NC}`$. We report exclusion limits for $`\alpha =0`$, $`\pi /4`$ and $`\pi /2`$. In the Compton scattering case, both space-space and space-time NC enter, leaving the two parameters, $`\alpha `$ and $`\gamma `$, in addition to $`\mathrm{\Lambda }_{NC}`$. We examine the two values $`\gamma =0`$ and $`\pi /2`$, and for each value of $`\gamma `$ give exclusion limits for $`\alpha =0`$, $`\pi /4`$ and $`\pi /2`$. We remind the reader that $`\alpha `$ relates to the direction of E, whereas $`\gamma `$ determines the orientation of B.

For the pair production process, the differential cross section for this process is given by:

$$\frac{d\sigma (\gamma \gamma f\overline{f})}{d\mathrm{cos}\theta d\varphi }=\frac{\alpha ^2}{2s}\left\{\frac{\widehat{u}}{\widehat{t}}+\frac{\widehat{t}}{\widehat{u}}4\frac{\widehat{t}^2+\widehat{u}^2}{\widehat{s}^2}\mathrm{sin}^2\left(\frac{k_1\theta k_2}{2}\right)\right\}.$$
(3)
The first two terms in the expression are the standard QED contributions, while the last term is due to the Feynman diagram with the cubic $`\gamma \gamma \gamma `$ coupling. The phase factor, $`\mathrm{sin}^2\left(\frac{k_1\theta k_2}{2}\right)`$ only appears in this new term. $`p_1`$ and $`p_2`$ are the momentum of the electron and positron, respectively, while $`k_1`$ and $`k_2`$ are the momenta of the incoming photons. $`\widehat{s}`$, $`\widehat{t}`$ and $`\widehat{u}`$ are the usual Mandelstam variables $`\widehat{s}=(k_1+k_2)^2`$, $`\widehat{t}=(k_1p_1)^2`$ and $`\widehat{u}=(k_1p_2)^2`$. The bilinear product in eqn. 6 simplifies to

$$\frac{1}{2}k_1\theta k_2=\frac{\widehat{s}}{4\mathrm{\Lambda }_{NC}^2}C_{03}.$$
(4)
The expression for the cross section is not Lorentz invariant due to the presence of the phase factor. Note that only space-time non-commutativity contributes and there is no $`\varphi `$ dependence in this case. As $`C_{03}=\mathrm{cos}\alpha `$, NCQED reproduces QED for pair production when $`\alpha =\pi /2`$, and also as $`\mathrm{\Lambda }_{NC}\mathrm{}`$.

The exclusion limits based on lepton pair production in $`\gamma \gamma `$ collisions and assuming an integrated luminosity of $`L=500fb^1`$ are summarized in Table I of Ref. us  for $`\alpha =0`$ and $`\pi /4`$. The values range from 220 GeV ($`\sqrt{s}=500GeV`$) to 2.7 TeV ($`\sqrt{s}=5.0TeV`$). These are based on the angular distribution which, as already noted, gives the highest limits. These limits could be improved by including three lepton generations in the final state and assuming some value for the lepton detection efficiency.

For the Compton scattering process, we find:

$$\frac{d\sigma (e^{}\gamma e^{}\gamma )}{d\mathrm{cos}\theta d\varphi }=\frac{\alpha ^2}{2s}\left\{\frac{\widehat{u}}{\widehat{s}}\frac{\widehat{s}}{\widehat{u}}+4\frac{\widehat{s}^2+\widehat{u}^2}{\widehat{t}^2}\mathrm{sin}^2\left(\frac{k_1\theta k_2}{2}\right)\right\}.$$
(5)
The first two terms in the expression are the standard, QED contribution, while the last term is due to the Feynman diagram with the cubic $`\gamma \gamma \gamma `$ coupling. As before, the phase factor only appears in this new term.

Here, $`p_1`$ and $`k_1`$ are the momenta of the initial state electron and photon, respectively, while $`p_2`$ and $`k_2`$ are the momenta of the final state electron and photon, respectively. $`\widehat{s}`$, $`\widehat{t}`$ and $`\widehat{u}`$ are the usual Mandelstam variables. In this case the phase factor simplifies to

$$\frac{1}{2}k_1\theta k_2=\frac{xk\sqrt{s}}{4\mathrm{\Lambda }_{NC}^2}[(C_{01}C_{13})\mathrm{sin}\theta \mathrm{cos}\varphi +(C_{02}+C_{23})\mathrm{sin}\theta \mathrm{sin}\varphi +C_{03}(1+\mathrm{cos}\theta )].$$
(6)
Compton scattering is sensitive to both $`\gamma `$ and $`\alpha `$, so it is complimentary to pair production studied here.

Fig. 1a shows the cross section $`\sigma `$ vs. $`\mathrm{\Lambda }_{NC}`$ for QED and NCQED with $`\alpha =0`$, $`\pi /4`$ and $`\pi /2`$, for a $`\sqrt{s}=0.5TeV`$ $`e^+e^{}`$ collider operating in $`e\gamma `$ mode. The QED (solid) curve includes the central QED value and $`\pm 1\sigma `$ bands (assuming 500 fb-1 of integrated luminosity). Fig. 1b shows the angular distribution, $`d\sigma /d\varphi `$, for QED and NCQED with $`\alpha =\pi /2`$, and $`\sqrt{s}`$ = $`\mathrm{\Lambda }_{NC}`$ = $`500GeV`$. The error bars in Fig. 1b assume 500 fb-1 of integrated luminosity. Note that there is no $`\varphi `$ dependence for $`\alpha =0`$ since for this case both E and B are parallel to the beam direction. In contrast, when $`\alpha =\pi /2`$, E is perpendicular to the beam direction which is reflected in the strong oscillatory behavior in the $`\varphi `$ distribution. The exclusion limits obtainable from Compton scattering are summarized in Table II of Ref. us  for $`L=500fb^1`$, with and without a 2% systematic error. The exclusion limits range from 545 GeV ($`\sqrt{s}=500GeV`$) to 7.4 TeV ($`\sqrt{s}=5.0TeV`$).

The pair production process is only sensitive to space-time NC and is therefore insensitive to $`\gamma `$. As $`\alpha `$ increases towards $`\pi /2`$ the deviations from SM decrease towards zero, with $`\alpha =\pi /2`$ being identical to the SM. On the other hand, the Compton scattering process is sensitive to both space-space and space-time NC as parametrized by $`\gamma `$ and $`\alpha `$. On the whole, we found that the Compton scattering process is superior to lepton pair production in probing NCQED. Despite significantly smaller statistics, the large modification of angular distributions (see Fig. 1b) leads to higher exclusion limits, well in excess of the center of mass energy for all colliders considered.


text:
# Centrality Dependence of Baryon and Meson Momentum Distributions in 𝑝𝐴 Collisions11footnote 1Talk presented in the 31th International Symposium on Multiparticle Dynamics, Datong, China

## Abstract

The proton and neutron inclusive distributions in the projectile fragmentation region of $`pA`$ collisions are studied in the valon model. Momentum degradation and flavor changes due to the nuclear medium are described at the valon level using two parameters. Particle production is treated by means of the recombination subprocess. Pion inclusive distributions can be calculated without any adjustable parameters.
PACS Number(s): 25.75.Dw, 13.85.Ni

The study of proton-nucleus $`(pA)`$ collisions is important because they are tractable intermediaries between $`pp`$ and $`AA`$ collisions, when intense interest exists in discovering the extent to which the dense medium created in an $`AA`$ collision differ from that of linear superpositions of $`pp`$ collisions. One of the properties of $`pA`$ collisions is that the momentum of the leading baryon in the projectile fragmentation region is degraded, a phenomenon commonly referred to, somewhat inappropriately, as baryon stopping. Such a transference of baryon number from the fragmentation to the central region contributes to the increase of matter density at mid-rapidity, thereby raising the likelihood of the formation of quark-gluon plasma. Thus it is important to understand the process of baryon momentum degradation and its dependence on nuclear size or centrality.

Since it is not feasible to perform first-principle parameter-free calculations of the momentum degradation at this point, experimental guidance is of crucial importance. Recently, several experiments have produced useful data on the subject, in particular, E910 and E941 at the AGS, and NA49 at the SPS . It is the $`x_F`$ dependences of the distributions of $`p\overline{p}`$ and $`n\overline{n}`$ that we shall focus on; moreover, their dependences on centrality will guide us in our determination of the nuclear effect on baryon momentum degradation.

Since it is questionable that the concept of color strings can be relevant in heavy-ion collisions where the abundance of color charges in the overlap region renders unlikely the development of constricted color flux tubes , our approach in this paper will be on the various levels of the constituents of the nucleon that are consistent with the parton model. More specifically, we shall use the valon model - to keep track of the momenta of the constituents and the recombination model  to describe the hadronization of the partons. The valons play a role in the collision problem as the constituent quarks do in the bound-state problem. Thus a nucleon has three valons which carry all the momentum of the nucleon, while each valon has one valence quark and its own sea quarks and gluons. Although soft processes are non-perturbative, the valon model (including recombination) nevertheless provides a systematic way of calculating all subprocesses that contribute to a particular inclusive process. Some of the subprocesses can be identified with certain diagrams in other approaches, e.g., baryon junction and diquark breaking terms . The valon distributions in the proton will be determined by fitting the parton distributions at low $`Q^2`$. The effect of the nuclear medium on the valon distribution of the proton projectile will involve two parameters, one characterizing the momentum degradation and the other flavor flipping. The color indices are all averaged over, since multiple gluon interactions, as the projectile traverses the target nucleus, are numerous and uncomputable. Their effects are, however, quantified in terms of the number, $`\nu `$, of target nucleons that participate in the $`pA`$ collision. For that reason it is important that the experimental data must have centrality selection expressible in terms of the average $`\overline{\nu }`$.

In the valon model a proton is considered to consist of three valons ($`UUD`$), which have the same flavors as the valence quarks ($`uud`$) that they individually contain. Thus a valon may be regarded as a parton cluster whose structure can be probed at high $`Q^2`$, but the structure of a nucleon itself in a low-$`p_T`$ scattering problem is described in terms of the valons. As in the parton model we work in a high-momentum frame so that it is sensible to use the momentum fractions of the constituents. Reserving $`x`$ for the momentum fraction of a quark, we use $`y`$ to denote the momentum fraction of a valon. In this paper $`y`$ never denotes rapidity. Let the exclusive valon distribution function be

$`G_{UUD}(y_1,y_2,y_3)=g(y_1y_2)^\alpha y_3^\beta \delta (y_1+y_2+y_31),`$ (1)
where $`y_1`$ and $`y_2`$ refer to the $`U`$ valons and $`y_3`$ the $`D`$ valon. The normalization factor $`g`$ is determined by requiring that the probability of finding these three valons in a proton be one, i. e.

$`{\displaystyle _0^1}𝑑y_1{\displaystyle _0^{1y_1}}𝑑y_2{\displaystyle _0^{1y_1y_2}}𝑑y_3G_{UUD}(y_1,y_2,y_3)=1.`$ (2)
Denote $`\frac{dz}{z}K_{NS}(z)`$ the distribution of the valence quark in a valon, and $`\frac{dz}{z}L(z)`$ that of the sea quarks, then the distribution for the favored quark (like $`u`$ quark in a $`U`$ valon) in a valon is $`K(z)K_{NS}(z)+L(z)`$. Thus the parton distribution in the proton can be written as

$`xu(x)`$ $`=`$ $`{\displaystyle 𝑑y_1𝑑y_2𝑑y_3G_{UUD}(y_1,y_2,y_3)\left[2K(x/y_1)+L(x/y_3)\right]},`$ (3)
$`xd(x)`$ $`=`$ $`{\displaystyle 𝑑y_1𝑑y_2𝑑y_3G_{UUD}(y_1,y_2,y_3)\left[2L(x/y_1)+K(x/y_3)\right]}`$ (4)
¿From previous studies, we use

$`K_{NS}(z)=z^a(1z)^b/B(a,b+1),`$ (5)
$`L(z)=\mathrm{}_o(1z)^5,`$ (6)
$`xu(x)`$ and $`xd(x)`$ have been given by CTEQ. By fitting $`xu(x)`$ and $`xd(x)`$ we can get $`\alpha ,\beta ,a,b,\mathrm{}_o`$

$`\alpha =0.70,\beta =0.25,a=0.79,b=0.26,\mathrm{}_0=0.083.`$
If we know the probability, $`F(x_1,x_2,x_3)`$, of finding a $`u`$ quark at $`x_1`$, another $`u`$ at $`x_2`$, and a $`d`$ at $`x_3`$, and the recombination function, $`R_p(x_1,x_2,x_3,x)`$, which is the probability for three quarks to form a proton at $`x`$, the momentum distribution of the produced proton in the projectile fragmentation region can be written as

$`{\displaystyle \frac{x}{\sigma _{in}}}{\displaystyle \frac{d\sigma ^p}{dx}}H_p(x)={\displaystyle \frac{1}{N}}{\displaystyle \frac{dx_1}{x_1}\frac{dx_2}{x_2}\frac{dx_3}{x_3}F(x_1,x_2,x_3)R_p(x_1,x_2,x_3,x)},`$ (7)
where $`N`$ is a normalization constant. From former studies

$`R_p(x_1,x_2,x_3,x)={\displaystyle \frac{x_1x_2x_3}{x^3}}G_{UUD}({\displaystyle \frac{x_1}{x}},{\displaystyle \frac{x_2}{x}},{\displaystyle \frac{x_3}{x}})`$ (8)
To calculate $`F(x_1,x_2,x_3)`$, we need to consider the physics involved in the passing of the incident proton through the nucleus. During the process, valons in the proton loss their energy and change their flavor contents. We describe the energy loss of valons by $`D(z_i,\nu _i)`$ through the following equation

$`G^{}(y_1^{},y_2^{},y_3^{})={\displaystyle \frac{dy_1dy_2dy_3}{y_1^{}y_2^{}y_3^{}}G(y_1,y_2,y_3)D(\frac{y_1^{}}{y_1},\nu _1)D(\frac{y_2^{}}{y_2},\nu _2)D(\frac{y_3^{}}{y_3}\nu _3)}`$ (9)
with $`G^{}(y_1^{},y_2^{},y_3^{})`$ depends on $`\nu _i`$ implicitly. The moments of $`D(z_i,\nu _i)`$ can be expressed, with a parameter $`\kappa `$, as

$`\stackrel{~}{D}(n,\nu _i)=\mathrm{exp}\left[\stackrel{~}{Q}(n)\nu _i\right]=\mathrm{exp}\left\{\kappa \nu _i\left[\psi (n)+\gamma _E\right]\right\}.`$
Assume that the flavor changes at each of the $`\nu _i`$ collisions are incoherent and the probability for a valon from $`U`$ to $`D`$ or from $`D`$ to $`U`$ in one valon-nucleon collision is $`q`$, then the probability for a valon to change its flavor after $`\nu _i`$ collisions can be proved to be $`q_{\nu _i}=\left\{1(12q)^{\nu _i}\right\}/2`$. Effectively, one can write valons as

$`U\stackrel{\nu _i}{}p_{\nu _i}U+q_{\nu _i}D,D\stackrel{\nu _i}{}p_{\nu _i}D+q_{\nu _i}U,`$ (10)
with $`p_{\nu _i}=1q_{\nu _i}`$, then $`F(x_1,x_2,x_3)`$ can be calculated. Therefore we can get the momentum distribution of produced proton in the projectile fragmentation region and compare the results with those from experiments. In doing so, one needs to calculate $`F(x_1,x_2,x_3)`$ for the produced protons and anti-protons and make average over all possible combinations of $`\nu _i`$, and choose parameters to fit the experimental data on $`p\overline{p}`$ production. We can also calculate the distribution for produced neutrons. By simultaneously fitting the experimental data on $`p\overline{p}`$ and $`n\overline{n}`$ in the proton fragmentation region, we get $`\kappa =0.62,q=0.37`$. This corresponds to an effective degradation length 17fm. From these two parameters, one can predict the momentum distribution for the pions in the same region. In above formulas, we used collision numbers at two different levels: $`\nu `$ at the nucleon level, and $`\nu _i`$ at the valon level. We need also to take into account the relation between the number of nucleon-nucleon collisions, $`\nu `$, and that of valon-nucleon collisions, $`\nu _i`$. Because of the space limit, the discussion is not contained in this talk, but can be found in Ref. . The fitted results and our prediction can also be found there.

In summary, the general shapes of momentum distribution of hadrons produced in the proton fragmentation region can be reproduced within the valon model with suitably chosen parameters and recombination function.


text:
# 1 Introduction

## 1 Introduction

Photon splitting is a process where the initial photon turns into a virtual electron-positron pair which scatters in the electric field of an atom and then transforms into two photons sharing the initial photon energy $`\omega _1`$. This is an example of a self-action of an electromagnetic field, which results also in such effects as coherent photon scattering (Delbrück scattering) and photon-photon scattering. The latter phenomenon was never observed experimentally.

Delbrück scattering was investigated in detail both theoretically and experimentally . It turned out that, for heavy atoms and high photon energy, its cross section calculated exactly in the parameter $`Z\alpha `$ ($`Z|e|`$ is the nucleus charge, $`\alpha =e^2/4\pi =1/137`$ is the fine-structure constant) drastically differs from that obtained in the lowest order in this parameter (Born approximation).

Recently, an essential progress in understanding of photon splitting phenomenon was achieved. In papers  various differential cross sections of high-energy photon splitting have been calculated exactly in the parameter $`Z\alpha `$. Similar to the case of Delbrück scattering, the exact cross section turns out to be noticeably smaller than that obtained in the Born approximation. So, the detailed experimental investigation of photon splitting provides a new sensitive test of QED when the effect of higher-order terms of the perturbation theory with respect to the external field is very important.

The observation of photon splitting is a difficult problem due to a smallness of its cross section as compared to those of other processes initiated by photons in a target. The following background processes are significant: double Compton effect off the atomic electrons ($`\gamma e\gamma \gamma e`$), and the emission of two hard photons from $`e^+e^{}`$ pair produced by the initial photon. The relative importance of these processes depends on the photon energy. For the energy $`\omega _1m`$, where the search of photon splitting was undertaken in two experiments , only double Compton scattering determines the background conditions. In these experiments, the photons from an intense radioactive source ($`Zn^{65}`$ with $`\omega _11.1`$ MeV in , and $`Co^{60}`$ with $`\omega _11.17,\mathrm{\hspace{0.17em}1.33}`$ MeV in ) were used. The combination of the coincidence and energy-summing detection technique was applied. The number of events considered as candidates for photon splitting exceeded the theoretical expectations by the factor of 300 in , and by the factor of 6 in .

At high photon energy $`\omega m`$ the emission of hard photons from $`e^+e^{}`$ pair becomes most important as a background process. In 1973 the experiment devoted to the study of Delbrück scattering of photons in energy region $`1÷7`$ GeV was performed . The bremsstrahlung non-tagged photon beam was used. Some events were assigned by authors of  to the photon splitting process. As shown in , the theoretical value for the number of photon splitting events under the conditions of the experiment was two orders of magnitude smaller than the experimental result. It was also argued that the events observed could be explained by the production of $`e^+e^{}`$ pair and one hard photon.

The first successful observation of photon splitting was performed in 1995-96 using the tagged photon beam of the energy $`120÷450`$ MeV at the VEPP-4M $`e^+e^{}`$ collider  in the Budker Institute of Nuclear Physics (Novosibirsk). Another goal of this experiment was a study of Delbrück scattering . The total statistics collected was $`1.610^9`$ incoming photons with BGO (bismuth germanate) target and $`410^8`$ without target for background measurements. The preliminary results were published in , . Here we present the new data analysis for this experiment.

## 2 Scheme of experiment

The experimental setup is shown in Fig. 1. Some ideas of this setup were suggested in . The main features of the experimental approach are:

* The use of high-quality tagged photon beam produced by backward Compton scattering of laser light off high-energy electron beam. Thereby, the energy of the initial photon is accurately determined.
* Strong suppression of the background processes by means of the detection of charged particles produced in the target and in other elements of the photon-beam line.
* The detection of both final photons to discriminate the photon splitting events from those with one final photon produced in Compton or Delbrück scattering.
* The requirement of the balance between the sum of the energies of the final photons and the energy of the tagged initial photon. This provides the additional suppression of the events with charged particles missed by the detection system.

At high energy of the initial photon $`\omega _1m`$, the photon splitting cross section is peaked at small angles between momenta of all photons ($`m/\omega _1`$). Therefore, a good collimation of the primary photon beam is required. The ROKK-1M facility  is used as the intense source of the tagged $`\gamma `$-quanta. The electron energy loss in the process of Compton scattering of laser light is measured by the tagging system (TS)  of the KEDR detector . The TS consists of the focussing spectrometer formed by accelerator quadrupole lenses, bending magnets, and 4 hodoscopes of the drift tubes. High-energy photons move in a narrow cone around the electron beam direction. The angular spread is of the order $`1/\gamma `$, where $`\gamma =E_{beam}/m`$ is the relativistic factor of the electron beam. The photon energy spectrum has a sharp edge at

$$\omega _{th}=\frac{4\gamma ^2\omega _{laser}}{1+4\gamma \omega _{laser}/m},$$
(1)
that allows one to perform the absolute calibration of the tagging system in a wide energy range. In experiment the laser photon energy was $`\omega _{laser}=1.165`$ eV, the electron beam energy $`E_{beam}=5.25`$ GeV, and $`\omega _{th}=450`$ MeV. The photon energy resolution provided by the tagging system depends on the photon energy and on the position of scattered electron in TS hodoscope: it was 0.8 % at $`\omega _1=450`$ MeV (at the center of the hodoscope) and $`5`$ % at $`\omega _1=120`$ MeV (at the edge of the hodoscope). The collimation of the photon beam is provided by two collimators spaced at 13.5 m. The last collimator, intended to strip off the beam halo produced on the first one, is made of four BGO (bismuth germanate) crystals as shown in the separate view in Fig. 1. After passing through the collimation system, the photon beam hits $`1X_0`$ thick BGO crystal target. In order to separate the photons passed the target without interaction from those scattered in the target, certain angular region around the photon beam direction ($`\theta 2.4`$ mrad) was enclosed by the dump. It is made of 13 $`X_0`$ thick BGO crystal installed in front of the photon detector. The only photons to be detected are those scattered outside the dump shadow. All active elements used in the beam line (collimators, target, dump, scintillating veto counter) set a veto signal in the trigger and their signals are used in the analysis for background suppression. The information from the target and beam dump is also used for measurement of the incoming photon flux. The liquid-krypton ionization calorimeter is used for the detection of the final photons. Its three-layer double-sided electrode structure enables one to get both (X and Y) coordinates for detected photons. The energy resolution of the calorimeter is $`2\%/\sqrt{\omega (\text{GeV})}`$. The liquid-krypton calorimeter is described in details in .

In the experiment, the detected final photons had the polar angles in the region $`2.4\text{mrad}\theta 20\text{mrad}`$. The corresponding cross section is called ”visible”. Fig. 2 shows the calculated energy dependence of the total ($`a`$) and visible ($`b`$) cross sections for various processes initiated by photons in BGO target: $`e^+e^{}`$ pair production, Compton scattering on atomic electrons, Delbrück scattering, and photon splitting. The calculation of the photon splitting cross section was performed using the results obtained in .

## 3 Results

In the event selection procedure the following constraints were applied:

* The absence of the signal caused by charge particles in all active elements of the photon-beam line.
* The balance of the tagged initial photon energy and the energy measured in the calorimeter within $`3\sigma `$ of its energy resolution.
* The existence of two separate tracks at least for one (X or Y) coordinate in the calorimeter strip structure.

The fulfillment of the latter requirement strongly suppresses the contribution of the processes with one photon in the final state which could imitate two photon events in the calorimeter.

The typical event which meets selection criteria is shown in Fig. 3. In this example two tracks are seen in both X and Y directions. The conversion of the first photon occurs in the Layer 1 while the second photon converts in the Layer 2.

The experimental results are presented in Tab. 1 and in Figs. 4, 5 together with the results of Monte-Carlo simulation based on the exact photon splitting cross section. The energy spectrum of the initial photons measured in the tagging system is shown in Fig.4(a). Tab. 1 and Fig. 5 present the data summed up over the initial photon spectrum. The errors shown in the Table are statistical ones. The systematic error is determined by the accuracy of the measurement of the number of initial photons and by the uncertainty in the reconstruction efficiency of photon splitting events. The estimation of these systematic errors gives 2 % and 5 %, respectively.

As seen from the Table 1 and Fig. 5a, the main part of the photon splitting events has, in agreement with the theory, a complanarity angle $`\phi `$ (the azimuth angle between final photon momenta) close to $`180^{}`$. The choice of the interval $`\phi >150^{}`$ allows us to improve the signal-to-background ratio (see, e.g., the last two rows of the Table 1). Just this $`\phi `$-interval was used to plot the distributions over polar angles in Fig. 5 and the dependence of the number of reconstructed photon splitting events on the initial photon energy $`E_{TS}`$ in Fig. 4(b). Note that for most of the events in this $`\phi `$-interval, the variable $`\stackrel{~}{x}=\theta _{min}/(\theta _{min}+\theta _{max})`$ is approximately equal to the ratio $`\text{min}(\omega _2,\omega _3)/\omega _1`$ since the main contribution to the cross section comes from the region $`|𝐤_2+𝐤_3|k_2,k_3`$, i.e. $`\phi 180^{}`$ and $`\omega _2\theta _2\omega _3\theta _3`$.

The results presented in the Table 1 and in Figs. 4(b), 5 show a good agreement between the theory and the experiment. More precisely, the total number of reconstructed events in the experiment (see Table 1) differs from the result of MC simulation by 1.5 standard deviations.

In order to demonstrate the role of the Coulomb corrections under the experimental conditions, we show in Fig. 6 the visible cross sections calculated exactly in $`Z\alpha `$ and in the Born approximation as a function of the initial photon energy. For all energies considered, the Born result exceeds the exact one by 20 % approximately. Therefore, the use of the Born cross section for MC simulation would lead to the disagreement between the theory and the experiment of 3.5 standard deviations. In other words, the experimental results are significantly closer to predictions of the exact theory than to those obtained in the Born approximation.

## 4 Conclusion

The results obtained confirm the existence of photon splitting phenomenon. They also make possible the quantitative comparison with the theoretical predictions. Moreover, the attained experimental accuracy allows one to distinguish between the theoretical predictions obtained with or without accounting for the Coulomb corrections. It turns out that the Coulomb corrections essentially improve the agreement between the theory and the experiment. We conclude that the experiment and the theory are consistent within the achieved experimental accuracy.

## Acknowledgments

We are grateful to the staff of the VEPP-4M accelerator complex for reliable work during the data taking. We thank A.N. Skrinsky and V.A. Sidorov for support of this experiment. Partial support by the RFBR (Grant 01-02-16926) is also gratefully acknowledged.