text:
# On the physical origin of dark matter density profiles

## 1 Introduction

The hierarchical clustering paradigm states that the growth of cold dark matter (CDM) haloes proceeds by accretion of smaller units from the surrounding environment, either by continuous infall or by discrete merging events.

Cosmological N-body simulations are a valuable tool to study the mass distribution of dark matter haloes and its evolution in the non-linear regime. In the early numerical work of Quinn et al. (1986) and Frenk et al. (1988), haloes showed an isothermal density profile ($`\rho r^2`$). Dubinski & Carlberg (1991) and Crone et al. (1994) had enough resolution in their simulations to detect the first evidence of departure from a pure power-law. Later on, Navarro et al. (1996, 1997, hereafter NFW) found that the density profile could be fitted by a simple analytical function

$$\rho (r)=\frac{\rho _\mathrm{s}}{(r/r_\mathrm{s})(1+r/r_\mathrm{s})^2}$$
(1)
in terms of a characteristic density $`\rho _\mathrm{s}`$ and a characteristic radius $`r_\mathrm{s}`$. This profile is steeper than isothermal at large radii and shallower near the centre. The logarithmic slope of the density profile, $`\alpha (r)\mathrm{d}\mathrm{log}(\rho )/\mathrm{d}\mathrm{log}(r)`$, tends to $`\alpha =3`$ for $`r\mathrm{}`$ and $`\alpha =1`$ for $`r0`$. It corresponds to the isothermal case, $`\alpha =2`$, at the characteristic radius only. Navarro et al. (1997) further showed that the two free parameters in equation (1) are not independent. Should this be true, the final mass distribution of objects of different scales could be described in terms of a one-parameter family of analytical profiles.

Similar results have been found in independent simulations with much higher mass and force resolution than the original NFW paper. However, there is still some controversy about the innermost value of $`\alpha `$ and its dependence on resolution. Moore et al. (1998, 1999), Ghigna et al. (1998, 2000) and Fukushige & Makino (1997, 2001) find steeper density profiles near the centre ($`\alpha 1.5`$), whereas other authors (Jing & Suto, 2000; Klypin et al., 2001) claim that the actual value of $`\alpha `$ may depend on halo mass, merger history, and substructure. Power et al. (2003) pointed out that the logarithmic slope becomes increasingly shallow inwards, with little sign of approaching an asymptotic value at the resolved radii. In that case, the precise value of $`\alpha `$ at a given cut-off scale would not be particularly meaningful. This result has been later confirmed by Fukushige et al. (2003) and Hayashi et al. (2003), and it is predicted by several analytical models (e.g. Taylor & Navarro, 2001; Hoeft et al., 2003).

Observed rotation curves of dwarf spiral and LSB galaxies (e.g. Flores & Primack, 1994; Moore, 1994; Burkert, 1995; Kravtsov et al., 1998; Borriello & Salucci, 2001; de Blok et al., 2001; de Blok & Bosma, 2002; Marchesini et al., 2002) seem to indicate that the shape of the density profile at small scales is significantly shallower than what is found in numerical simulations. This discrepancy has been often signalled as a genuine crisis of the CDM scenario, and several alternatives have been suggested, such as warm (Colín et al., 2000; Sommer-Larsen & Dolgov, 2001), repulsive (Goodman, 2000), fluid (Peebles, 2000), fuzzy (Hu et al., 2000), decaying (Cen, 2001), annihilating (Kaplinghat et al., 2000), or self-interacting (Spergel & Steinhardt, 2000; Yoshida et al., 2000; Davé et al., 2001) dark matter.

Unfortunately, it has proved remarkably hard to establish the inner slope of the dark matter distribution observationally (see e.g. Swaters et al., 2003). Some authors (van den Bosch & Swaters, 2001; Jimenez et al., 2003; Swaters et al., 2003) claim that a cuspy density profile with $`\alpha 1`$ is consistent with current observations, although a shallower slope is able to explain them as well. Yet, a value as steep as $`\alpha =1.5`$ can be confidently ruled out in most cases. According to Hayashi et al. (2003), only about 30 per cent of the rotational curves are actually inconsistent with the simulation data.

On cluster scales, X-ray analyses have led to wide ranging results, from $`\alpha =0.6`$ (Ettori et al., 2002) to $`\alpha =1.2`$ (Lewis et al., 2003) or even $`\alpha =1.9`$ (Arabadjis et al., 2002). Measurements based on gravitational lensing yield conflicting estimates as well, either in rough agreement with the results of numerical simulations (e.g. Dahle et al., 2003; Gavazzi et al., 2003), or finding much shallower slopes, $`\alpha =0.5`$ (e.g. Sand et al., 2002, 2003).

A conclusive theoretical prediction of the central mass distribution of CDM haloes is therefore an important check for any model of structure formation. The controversy regarding the ’universal’ density profile and its logarithmic slope at the centre has stimulated a great deal of analytical work. On one hand, we would like to find out not only the actual shape of the profile, but also the physical mechanisms behind the ’universality’ observed in numerical N-body simulations. On the other hand, it would be interesting to find and explain correlations with halo size, environment, power spectrum or even the nature of dark matter particles.

A number of plausible arguments about the radial structure of CDM haloes have been advanced during the last 30 years. The basic problem of the collisionless collapse of a spherical perturbation in an expanding background was first addressed in the two seminal papers by Gunn & Gott (1972) and Gunn (1977), where the cosmological expansion and the role of adiabatic invariance were first introduced in the context of the formation of individual objects. The next step was accomplished by Fillmore & Goldreich (1984) and Bertschinger (1985), who found analytical predictions for the density of collapsed objects seeded by scale-free primordial perturbations in a flat universe. Hoffman & Shaham (1985) generalised these solutions to realistic initial conditions in flat as well as open Friedmann models. Modifications of the self-similar collapse model to include more realistic dynamics of the growth process have been proposed (e.g. Avila-Reese et al., 1998; Henriksen & Widrow, 1999; Łokas, 2000; Kull, 1999; Subramanian et al., 2000). Several authors (e.g. Syer & White, 1998; Salvador-Sole et al., 1998; Nusser & Sheth, 1999; Manrique et al., 2003) argue that the central density profile is linked to the merging history of dark matter substructure, and baryons have been invoked both to shallow (e.g. El-Zant et al., 2001; El-Zant et al., 2003) and to steepen (Blumenthal et al., 1986) the dark matter profile.

In this paper, we present an analytical model for the assembly of CDM haloes based on the spherical collapse paradigm. Following the spirit of Hoffman & Shaham (1985), we will assume that objects do not form around *local* maxima of the primordial density field, but of the *smoothed* density field. We argue that all information about the initial conditions below the smoothing scale is lost during the merging process. We will show that setting the smoothing scale for a halo of a given mass is equivalent to specifying its formation time. Initial conditions are computed following the Gaussian random peaks statistics described by Bardeen et al. (1986, hereafter BBKS), and angular momentum is included in a phenomenological way. Our theoretical predictions are then compared with the results of high-resolution numerical simulations, showing that this model is able to accurately reproduce the mass distribution of individual objects.

The paper is structured as follows. Details of our implementation of spherical collapse are given in Section 2. Numerical experiments are described in Section 3, where both the mass and velocity distributions are investigated. We discuss our results in Section 4, and Section 5 summarises our main conclusions.

## 2 Spherical collapse

The assembly of dark matter haloes is a highly non-linear process, and strong simplifying assumptions must be made in order to tackle the problem analytically. Traditionally, there are two complementary paradigms: the spherical infall model (Gunn & Gott, 1972) and the Press-Schecter formalism (Press & Schechter, 1974), in which mergers play a dominant role (see e.g. Nusser & Sheth, 1999; Manrique et al., 2003).

### 2.1 The model

The most simple way of addressing the problem of structure formation is to assume spherical symmetry. As shown by Tolman (1934) and Bondi (1947), a spherically symmetric solution of the Einstein equations can be easily interpreted in terms of Newtonian dynamics. The equation of motion for a Lagrangian shell enclosing a mass $`M`$ can be derived from the conservation of energy

$$ϵ(r)\frac{E(r)}{m}\frac{\dot{r}^2}{2}G\frac{M(r)+\frac{4\pi }{3}\rho _\mathrm{\Lambda }r^3}{r}=ϵ_\mathrm{i}(r)$$
(2)
where $`r(t)`$ is the Lagrangian coordinate, $`M(r)=M_\mathrm{i}(r_\mathrm{i})`$ is the mass contained within the shell, and $`\rho _\mathrm{\Lambda }\frac{\mathrm{\Lambda }}{8\pi G}`$ is the vacuum energy density associated to a cosmological constant. Throughout this paper, the subscript ’i’ will be used to denote initial conditions.

If we assume the universe to be homogeneous, the mass is given by $`M_\mathrm{i}=\frac{4\pi }{3}\mathrm{\Omega }_\mathrm{m}^i\rho _c^ir_\mathrm{i}^3`$, where $`\rho _c^i\frac{3H_\mathrm{i}^2}{8\pi G}`$. Substituting in (2) and making a coordinate transformation

$$r(t)r_\mathrm{i}\alpha (t)$$
(3)
we obtain the Friedmann equation without radiation:

$$\frac{\dot{\alpha }^2}{H_\mathrm{i}^2}\mathrm{\Omega }_\mathrm{m}^i\alpha ^1\mathrm{\Omega }_\mathrm{\Lambda }^i\alpha ^2=1\mathrm{\Omega }_\mathrm{m}^i\mathrm{\Omega }_\mathrm{\Lambda }^i.$$
(4)
In a homogeneous universe, $`\alpha (t)`$ is independent of $`r_\mathrm{i}`$ and plays the role of a uniform cosmic expansion factor. Our model assumes that structures grow from spherically symmetric perturbations, defined as

$$\mathrm{\Delta }_\mathrm{i}(r_\mathrm{i})\frac{M(r_\mathrm{i})}{\frac{4\pi }{3}\mathrm{\Omega }_\mathrm{m}^i\rho _c^ir_\mathrm{i}^3}1$$
(5)
with $`0<\mathrm{\Delta }_\mathrm{i}1`$. These perturbations are introduced at an early epoch $`a_\mathrm{i}`$ by slightly displacing the spherical shells of matter. Keeping only terms to first order in $`\mathrm{\Delta }_\mathrm{i}(r_\mathrm{i})`$, the new positions and velocities are given by

$$r_\mathrm{i}^{}r_\mathrm{i}\left(1\frac{1}{3}\mathrm{\Delta }_\mathrm{i}(r_\mathrm{i})\right);v_\mathrm{i}^{}H_\mathrm{i}r_\mathrm{i}\left(1\frac{2}{3}\mathrm{\Delta }_\mathrm{i}(r_\mathrm{i})\right).$$
(6)
The mass enclosed by the perturbed shell is still $`M_\mathrm{i}`$, and the initial specific energy (using $`\mathrm{\Omega }_\mathrm{m}^i1`$ and $`\mathrm{\Omega }_\mathrm{\Lambda }^i0`$) is $`ϵ_\mathrm{i}\frac{5}{6}\mathrm{\Delta }_\mathrm{i}(H_\mathrm{i}r_\mathrm{i})^2`$. Energy conservation (2) leads to the equation of motion

$$\frac{5}{3}\mathrm{\Delta }_\mathrm{i}\frac{\dot{\alpha }^2}{H_\mathrm{i}^2}\mathrm{\Omega }_\mathrm{m}^i\alpha ^1\mathrm{\Omega }_\mathrm{\Lambda }^i\alpha ^2.$$
(7)
According to this equation, the evolution of a single spherical shell would be similar to that of a Friedmann universe. For high enough values of $`\mathrm{\Delta }_\mathrm{i}`$, the shell reaches a maximum radius $`r_\mathrm{m}`$ at a turn-around time $`t_\mathrm{m}`$ and then re-collapses. In an Einstein-de Sitter universe ($`\mathrm{\Omega }_\mathrm{m}=1`$, $`\mathrm{\Omega }_\mathrm{\Lambda }=0`$), equation (7) can be solved analytically, yielding

$$r_\mathrm{m}=\frac{3r_\mathrm{i}}{5\mathrm{\Delta }_\mathrm{i}};t_\mathrm{m}=\frac{\pi }{2H_\mathrm{i}(5\mathrm{\Delta }_\mathrm{i}/3)^{3/2}}.$$
(8)
This is also a valid approximation for shells reaching their turn-around radius $`r_\mathrm{m}`$ before the cosmological constant term starts to dominate the expansion. Since the shells need at least another $`t_\mathrm{m}`$ to virialise, expression (8) can be applied in a $`\mathrm{\Lambda }`$CDM universe to estimate the maximum expansion radius and time for the innermost part of a virialised halo. For the outer shells, the equation of motion (7) must be integrated numerically in order to find the trajectories up to the maximum radius.

In the absence of shell-crossing, shells would reach the origin at $`T=2t_\mathrm{m}`$. Since they are assumed to be composed of collisionless CDM particles, they would simply pass through the centre, describing an oscillatory motion with amplitude $`r_\mathrm{m}`$ and period $`T`$. However, equation (7) holds as long as the enclosed mass $`M_\mathrm{i}`$ remains constant. As a shell re-collapses, its particles will cross the orbits of inner shells, and energy will not be a constant of motion any more.

After turn-around, our model assumes that the particles belonging to a shell oscillate (or, more generally, orbit) within the gravitational potential of the dark matter halo. Since CDM particles are expected to spend most of the time in the outermost regions of their orbits (particularly when angular momentum is taken into account), we approximate the mass distribution of the halo by a simple power law

$$M(r)=M_\mathrm{i}\left(\frac{r}{r_\mathrm{m}}\right)^{\alpha _\mathrm{M}(r_\mathrm{m})}$$
(9)
where $`\alpha _\mathrm{M}(r)\mathrm{d}\mathrm{log}M(r)/\mathrm{d}\mathrm{log}r`$ is the local value of the logarithmic slope of the mass profile, evaluated at the maximum radius of the orbit.

At first sight, this might seem similar to the classical approach based on self-similarity (Bertschinger, 1985), but in that case the final mass distribution is indeed assumed to be a power law, whereas in our model this ansatz is only an approximation to compute the local potential. The final mass profile is obtained self-consistently, adding the contributions from all shells, each of them with an individual value of $`\alpha _\mathrm{M}(r_\mathrm{m})`$.

The probability of finding a particle inside radius $`r`$ is proportional to the fraction of time it spends within $`r`$:

$$P(r,r_\mathrm{m})=\frac{1}{t_\mathrm{m}}_0^r\frac{\mathrm{d}x}{v_r(x)}=\frac{1}{t_\mathrm{m}}_0^r\frac{\mathrm{d}x}{\sqrt{\mathrm{\Phi }(r_\mathrm{m})\mathrm{\Phi }(x)}}.$$
(10)
We evaluate numerically the value of $`\alpha _\mathrm{M}(r_\mathrm{m})`$ in order to compute the potential. Taking different prescriptions, such as an isothermal profile ($`\alpha _\mathrm{M}=1`$ for every shell) or a Keplerian potential ($`\alpha _\mathrm{M}=0`$) does not lead to qualitative variations in the probability $`P(r,r_\mathrm{m})`$ and the resulting mass distribution.

If phase mixing is considered to be efficient, particles initially on the same shell will be spread out from $`r=0`$ to $`r=r_\mathrm{m}`$ a short time after $`t_\mathrm{m}`$. For the sake of simplicity, we will consider that phase mixing is instantaneous, so immediately after turn-around the shell is transformed into a density distribution whose cumulative mass is proportional to $`P(r,r_\mathrm{m})`$.

This means that recently accreted particles contribute to the mass within $`x_\mathrm{m}<r_\mathrm{m}`$ (i.e. shell-crossing). For shells whose maximum radius was $`x_\mathrm{m}`$, the enclosed mass changes from $`M_\mathrm{i}(x_\mathrm{m})`$ to

$$M(x_\mathrm{m})=M_\mathrm{i}(x_\mathrm{m})+M_{\mathrm{add}}(x_\mathrm{m})$$
(11)
where $`M_{\mathrm{add}}(x_\mathrm{m})`$ accounts for particles belonging to outer shells. To compute $`M_{\mathrm{add}}(x_\mathrm{m})`$ (see Zaroubi & Hoffman, 1993), we must integrate the contribution of every shell whose maximum radius is larger than $`x_\mathrm{m}`$, up to the current turn-around radius $`R_{\mathrm{ta}}`$:

$$M_{\mathrm{add}}(x_\mathrm{m})=_{x_\mathrm{m}}^{R_{\mathrm{ta}}}\frac{\mathrm{d}M_\mathrm{i}(r)}{\mathrm{d}r}P(x_\mathrm{m},r)dr.$$
(12)
To compute the evolution of the shell after shell-crossing, we apply adiabatic invariance (Gunn, 1977). If the potential evolves on a timescale much longer than the orbital period of the inner particles, their dynamics admits an adiabatic invariant

$$J_r=\frac{1}{2\pi }v_r(r)dr$$
(13)
also known as the radial action. For a power-law potential, the radial action $`J_r`$ is proportional to $`\sqrt{x_\mathrm{m}M(x_\mathrm{m})}`$. When we increase the mass by an amount $`M_{\mathrm{add}}`$, the maximum radius of the inner shell must decrease in order to keep $`J_r`$ constant. The final radius, $`x_0`$, is given by the implicit equation $`x_0=F(x_0)x_\mathrm{m}`$, where

$$F(x_0)=\frac{M_\mathrm{i}}{M_\mathrm{i}+M_{\mathrm{add}}(x_0)}=\frac{M_\mathrm{i}}{M_\mathrm{i}+M_{\mathrm{add}}[F(x_0)x_\mathrm{m}]}$$
(14)
and whose solution must be obtained numerically for each shell.

To summarise, the numerical procedure to compute the final radius $`r_0(t_0)`$ of a Lagrangian shell of matter involves the following steps:

1. Set $`M_\mathrm{i}`$ (or, equivalently, $`r_\mathrm{i}`$). Start by the outer shell.
2. Integrate the equation of motion (7) up to $`t_\mathrm{m}`$.
3. If $`t_\mathrm{m}>t_0`$, the shell is still expanding: $`r_0=r(t_0)`$.
4. If $`t_\mathrm{m}<t_0`$, solve (14) to compute $`r_0=F(r_0)r_\mathrm{m}`$, and add the contribution of this shell to $`M_{\mathrm{add}}(r)`$.
5. Repeat for the next shell towards the centre.

### 2.2 Initial conditions

The model described above allows us to compute the mass distribution arising from a primordial fluctuation $`\mathrm{\Delta }_\mathrm{i}(r_\mathrm{i})`$, but does not say anything about the shape of this function or its physical origin. Never the less, it is important to note that the final density profile is entirely determined by this initial condition. Thus, in the spherical collapse paradigm, the case for a universal density profile can be reformulated in terms of universality in the primordial fluctuations that set $`\mathrm{\Delta }_\mathrm{i}(r_\mathrm{i})`$.

Hoffman & Shaham (1985) suggested that, according to the hierarchical scenario of structure formation, haloes should collapse around maxima of the smoothed density field. The statistics of peaks in a Gaussian random field has been extensively studied in the classical paper by BBKS. A well-known result is the expression for the radial density profile of a fluctuation centered on a primordial peak of arbitrary height:

$$\frac{\delta _\mathrm{f}(r)}{\sigma _0}=\nu \psi (r)\frac{\theta (\gamma ,\gamma \nu )}{\gamma (1\gamma ^2)}\left[\gamma ^2\psi (r)+\frac{R_{}^2}{3}^2\psi (r)\right]$$
(15)
where $`\psi (r)\xi (r)/\sigma _0`$ is the normalised two-point correlation function, $`\sigma _0\xi (0)^{1/2}`$ is the rms density fluctuation, and $`\nu \sigma _0`$ is the height of the peak. The quantities $`\gamma \sigma _1^2/(\sigma _2\sigma _0)`$ and $`R_{}\sqrt{3}\sigma _1/\sigma _2`$ are related to the moments of the power spectrum,

$$\sigma _j^2\frac{1}{2\pi ^2}_0^{\mathrm{}}P(k)k^{2(j+1)}dk,$$
(16)
and the function $`\theta (\gamma ,\gamma \nu )`$ parametrizes the second derivative of the density fluctuation near the peak. BBKS suggest the approximate fitting formula

$$\theta (\gamma ,\gamma \nu )\frac{3(1\gamma ^2)+(1.2160.9\gamma ^4)\mathrm{exp}[\frac{\gamma }{2}(\frac{\gamma \nu }{2})^2]}{\left[3(1\gamma ^2)+0.45+\frac{\gamma \nu }{2}\right]^{1/2}+\frac{\gamma \nu }{2}}$$
(17)
valid to 1 per cent accuracy in the range $`0.4<\gamma <0.7`$ and $`1<\gamma \nu <3`$, which is the scale relevant for both galaxies and galaxy clusters.

Expression (15) is often quoted in the literature (e.g. Hoffman, 1988; Łokas & Hoffman, 2000; Del Popolo et al., 2000; Hiotelis, 2002) as the initial condition $`\mathrm{\Delta }_\mathrm{i}(r_\mathrm{i})`$. However, we argue that $`\delta _\mathrm{f}(r)`$ denotes the smoothed density profile around a local maximum of the smoothed density field,

$$\delta _\mathrm{f}(\stackrel{}{r})=W_\mathrm{f}(\stackrel{}{r}\stackrel{}{x})\delta (\stackrel{}{x})\mathrm{d}^3\stackrel{}{x}$$
(18)
where the function $`W_\mathrm{f}(\stackrel{}{r}\stackrel{}{x})`$ is a smoothing kernel that depends on a certain filtering scale $`R_\mathrm{f}`$.

The smoothed profile $`\delta _\mathrm{f}(r)`$ given by (15) is in general not equal to the mean value of the actual overdensity, $`\delta (r)`$, that must be integrated to compute $`\mathrm{\Delta }_\mathrm{i}(r_\mathrm{i})`$:

$$\mathrm{\Delta }_\mathrm{i}(r_\mathrm{i})=4\pi _0^{r_\mathrm{i}}\delta (\stackrel{}{x})x^2dx.$$
(19)
Comparing this expression with (18), we see that $`\mathrm{\Delta }_\mathrm{i}(r_\mathrm{i})`$ is equivalent to $`\delta _\mathrm{f}(0)`$ as long as $`W_\mathrm{f}`$ is taken to be a spherical top hat of radius $`r_\mathrm{i}`$.

To sum up, we are interested in the physical density profile $`\mathrm{\Delta }_\mathrm{i}(r_\mathrm{i})`$ around a local maximum of the smoothed density field. To locate the maximum, we use a Gaussian smoothing kernel

$$W_\mathrm{f}(\stackrel{}{r}\stackrel{}{x})=(2\pi R_\mathrm{f}^2)^{3/2}\mathrm{exp}\left(\frac{|\stackrel{}{r}\stackrel{}{x}|^2}{2R_\mathrm{f}^2}\right)$$
(20)
in order to avoid the oscillations in Fourier space that arise from a top hat filter. We set the scale of the fluctuation, $`R_\mathrm{f}`$, and impose the condition that $`r=0`$ is a maximum of $`\delta _\mathrm{f}(\stackrel{}{r})`$.

Then, we compute $`\mathrm{\Delta }_\mathrm{i}(r_\mathrm{i})=\delta _{r_\mathrm{i}}(0)`$ by applying a top hat smoothing of radius $`r_\mathrm{i}`$. BBKS show that the probability distribution of $`\delta _{r_\mathrm{i}}(0)`$ is a Gaussian with mean

$$\delta _{r_\mathrm{i}}(0)=\nu _\mathrm{f}\frac{\sigma _{0h}^2}{\sigma _{0f}}\frac{\gamma _\mathrm{f}\theta (\gamma _\mathrm{f},\nu _\mathrm{f})}{1\gamma _\mathrm{f}^2}\frac{\sigma _{0h}^2}{\sigma _{0f}}\left(1\frac{\sigma _{1h}^2\sigma _{0f}^2}{\sigma _{0h}^2\sigma _{1f}^2}\right).$$
(21)
The moments

$$\sigma _{jx}^2\frac{1}{2\pi ^2}_0^{\mathrm{}}P_x(k)k^{2(j+1)}dk$$
(22)
are computed from

$$P_\mathrm{f}(k)P(k)\mathrm{exp}\left[(kR_\mathrm{f})^2\right]$$
(23)
and

$$P_h(k)P(k)\mathrm{exp}\left[\frac{(kR_\mathrm{f})^2}{2}\right]\frac{3[\mathrm{sin}(kr_\mathrm{i})kr_\mathrm{i}\mathrm{cos}(kr_\mathrm{i})]}{(kr_\mathrm{i})^3}$$
(24)
where $`P(k)`$ is the $`\mathrm{\Lambda }\mathrm{CDM}`$ power spectrum that we used to generate the initial conditions for our simulations, evaluated at time $`a_\mathrm{i}`$.

### 2.3 Angular momentum

Two decades after the seminal paper by Gunn & Gott (1972), White & Zaritsky (1992) pointed out (see also the pioneering work of Gurevich & Zybin, 1988) that angular momentum would prevent the orbits of CDM particles from reaching the origin.

For pure radial orbits, the mass in the centre is dominated by $`M_{\mathrm{add}}`$ (i.e. particles from the outer shells) when the profile at turn-around is shallower than isothermal (Fillmore & Goldreich, 1984). The final density distribution is $`\rho (r)r^2`$, independent on the initial logarithmic slope. More recently, Łokas (2000) and Łokas & Hoffman (2000) found a similar result considering non-self-similar primordial fluctuations based on the peak formalism.

Angular momentum arises in linear theory (White, 1984) from the misalignment between the asymmetry of the collapsing region (i.e. the inertia tensor) and the tidal forces it experiences during the linear regime. In addition, violent relaxation (Lynden-Bell, 1967) transforms the radial infall energy into random velocity dispersion, which has a tangential component.

The total amount of angular momentum acquired by the system is, however, an open question. The usual approach (Avila-Reese et al., 1998; Nusser, 2001; Hiotelis, 2002) consists in assigning a specific angular momentum at turn-around

$$j\sqrt{GMr_\mathrm{m}}.$$
(25)
With this prescription, the orbit eccentricity $`e`$ is the same for all particles in the halo (Nusser, 2001). The pericentric radius is given by

$$r_{\mathrm{min}}=\frac{1e}{1+e}r_{\mathrm{max}}$$
(26)
where $`r_{\mathrm{max}}`$ is computed from adiabatic invariance (14), following the procedure explained in Section 2.1. Angular momentum is taken into account by adding the usual term $`j^2/(2r^2)`$ to the gravitational potential $`\mathrm{\Phi }(r)`$ in equation (10). Although this changes the actual value of the radial action, $`J_r`$ is still proportional to $`\sqrt{r_\mathrm{m}M(r_\mathrm{m})}`$ and (14) can be used to compute the collapse factor $`F(r)`$. The assumption of spherical symmetry leads to angular momentum conservation, which implies constant $`e`$ during the contraction.

The effect of angular momentum is illustrated in Figure 1. We plot the mass distribution at $`z=0`$ corresponding to a $`3\sigma `$ peak in the primordial density field, smoothed on $`1h^1`$ Mpc scales. As noted by White & Zaritsky (1992) and Hiotelis (2002), angular momentum plays a key role during secondary infall, decreasing the amount of mass contributed by recently accreted shells in the innermost regions of the halo.

In agreement with Hiotelis (2002), we find that the predicted density profile becomes shallower in the centre as the amount of angular momentum is increased (lower $`e`$). Pure radial orbits give rise to a steep profile, somewhat similar to the form proposed by Moore et al. (1999), although the exact shape is slightly different. When the eccentricity is low, the halo develops a constant density core inconsistent with the results of numerical simulations. A value $`e=0.5`$ yields a mass distribution virtually indistinguishable from the NFW formula (dotted line in Figure 1) down to kpc scales. However, the logarithmic slope predicted by the model keeps increasing towards the centre, while NFW tends to an asymptotic value, $`\alpha =1`$.

## 3 Numerical Experiments

In order to test the analytical model, we carried out a series of high-resolution N-body simulations of cluster formation with the adaptive mesh code ART (Kravtsov et al., 1997). Slightly lower spatial resolution experiments have been run from the same initial conditions with the Tree-SPH gasdynamical code Gadget (Springel et al., 2001; Springel & Hernquist, 2002). The radial structure of both dark and baryonic components has been addressed in Ascasibar et al. (2003). For a detailed description of the numerical experiments, the reader is referred to Ascasibar (2003).

A sample of cluster-size haloes has been selected from an initial low-resolution ($`128^3`$ particles) pure N-body simulation of a 80 $`h^1`$ Mpc box in a $`\mathrm{\Lambda }`$CDM universe ($`\mathrm{\Omega }_\mathrm{m}=0.3`$; $`\mathrm{\Omega }_\mathrm{\Lambda }=0.7`$; $`h=0.7`$; $`\sigma _8=0.9`$). Higher resolution has been achieved by means of the multiple-mass technique (see Klypin et al., 2001, for details), using 3 levels of mass refinement. This is equivalent to an effective resolution of $`512^3`$ CDM particles ($`3.1\times 10^8h^1`$ M) in the highest refinement level. The minimum cell size allowed in the ART runs was $`1.2h^1`$ kpc. All simulations were started at $`z=50`$.

This procedure has been applied to 15 objects, ranging from $`3\times 10^{13}`$ to $`2\times 10^{14}h^1`$ M. In order to explore a broader mass range, a smaller box (25 $`h^1`$ Mpc) has been simulated, and six galaxy-size haloes have been added to the cluster sample. Mass resolution is $`1.2\times 10^6h^1`$ M for these objects, with a minimum cell size of 0.2 $`h^1`$ kpc.

All CDM haloes have been classified according to their dynamical state according to a substructure-based criterium. We use the Bound Density Maxima galaxy finding algorithm (see e.g. Colín et al., 1999; Klypin et al., 1999) and look for the most massive subhalo within the virial radius. We label as major merger any object in which the subhalo is heavier than 50 per cent the mass of the main object; if the mass is between 10 and 50 per cent, it is classified as a minor merger; otherwise, we assume the object to be a relaxed system in virial equilibrium.

### 3.1 Mass distribution

The spherically-averaged mass distribution of our numerical CDM haloes is shown in Figure 2. All profiles have been rescaled by their best-fitting characteristic density and radius. As can be seen in the Figure, the NFW formula provides a reasonable approximation to the mass distribution, although some deviations occur at the innermost part, as well as beyond the virial radius.

In Figure 3, the logarithmic slopes of both mass and density profiles are plotted according to our dynamical classification. We do not find any evidence of an asymptotic behavior up to the resolution limit, in agreement with recent numerical studies (Power et al., 2003; Fukushige et al., 2003; Hayashi et al., 2003). An extrapolation of our results suggests that the central slopes in relaxed haloes could be indeed less steep than the NFW fit, as predicted by analytical models based on the velocity dispersion profile (Taylor & Navarro, 2001; Hoeft et al., 2003). More resolution is clearly needed in order to establish a firm conclusion.

As pointed out by several authors (e.g. Jing & Suto, 2000; Klypin et al., 2001; Ascasibar et al., 2003), the mass distribution near the centre might depend on the dynamical state of the halo. At the resolution of the present simulations, relaxed haloes are well described by NFW model, but merging systems display steeper profiles. A conjecture that deserves further investigation is whether violent relaxation and subsequent equipartition of energy may drive the inner part of the system into a nearly isothermal state. That would provide a natural mechanism for the final profile to be independent on the details of the merging history.

Another difference between merging and relaxed systems is that mergers are typically more extended than relaxed haloes of the same mass. We plot in Figure 4 the mass-concentration relation for our sample, compared to the model of Bullock et al. (2001),

$$c(M,a)=K\frac{a}{a_c(M)}$$
(27)
where the collapse time $`a_c`$ is defined as the epoch at which the typical collapsing mass, $`M_{}(a_c)`$, equals a fixed fraction $`F`$ of the halo mass at epoch $`a`$. According to Bullock et al. (2001), we set $`K=3`$ and $`F=0.001`$ in order to account for the most massive haloes. The scatter around the relation is fitted by $`K=2`$ and $`K=4.7`$.

Our sample of dark matter haloes is consistent with the expected trend, although our cluster-size haloes seem to be slightly more concentrated than the theoretical model. On the other hand, major mergers display systematically lower concentrations, since their collapse time is very close to the present. However, the value of $`a_c`$ is not well defined in a merging system. In particular, during the first stages of the merger the profile corresponds to an old object, perturbed by an approaching satellite. It is only after virialisation when the profile relaxes to its final form, and $`a_c`$ increases accordingly.

### 3.2 Angular momentum

In addition to the radial mass distribution, the kinetic structure of numerical haloes can offer interesting insights into the formation of galaxies and galaxy clusters. In particular, we are interested in the specific angular momentum of dark matter particles in order to set the eccentricity in the analytical model.

We separate the velocity field into a random component (i.e. velocity dispersion) and ordered rotation (i.e. average $`j`$). The velocity dispersion of a collisionless CDM halo is related to the total mass distribution by virtue of the Jeans equation (see e.g. Binney & Tremaine, 1987). For a spherically symmetric system with isotropic velocity dispersion, neglecting infall,

$$\frac{1}{\rho }\frac{\mathrm{d}(\rho \sigma _r^2)}{\mathrm{d}r}=\frac{GM}{r}$$
(28)
where $`\rho `$ is the local density, $`\sigma _r`$ is the radial velocity dispersion, and $`M`$ is the mass enclosed within radius $`r`$. An approximate velocity dispersion profile could be derived by substituting a given mass distribution (e.g. NFW) and setting an arbitrary normalisation $`\sigma _r(0)`$. The contribution of random motions to the angular momentum of the CDM particles would be given by the tangential velocity dispersion. In the isotropic case, this amounts to $`j^2=2r^2\sigma _r^2(r)`$. Note that random motions do not contribute to the total angular momentum of the halo (i.e. $`j=0`$).

A more empirical approach to the velocity dispersion profile has been followed by Taylor & Navarro (2001). They realised that the coarse-grained phase-space density of a sample of numerical galaxy-size haloes followed a power law

$$\frac{\rho }{\sigma ^3}r^\beta $$
(29)
with $`\beta =1.875`$ over more than two and a half decades in radius. Rasia et al. (2003) obtained a similar result for cluster-size haloes, although their best-fitting slope is $`\beta =1.95`$. The normalisation, though, has not been given in any of the two studies.

We have investigated the phase-space structure of the CDM haloes in our sample, including both galaxies and galaxy clusters. The average profile $`\rho /\sigma ^3`$ is plotted in Figure 5. We find that all our haloes are well described by

$$\frac{\rho }{\sigma ^3}=10^{1.46\pm 0.04}\frac{\rho _c}{v_{\mathrm{vir}}^3}\left(\frac{r}{r_{\mathrm{vir}}}\right)^{1.90\pm 0.05}$$
(30)
where $`v_{\mathrm{vir}}^2GM_{\mathrm{vir}}/r_{\mathrm{vir}}`$. This result is in fair agreement with the slopes reported by Taylor & Navarro (2001) and Rasia et al. (2003). The scatter around the average profile is remarkably low, taking into account that our haloes span four orders of magnitude in mass, and that we considered all systems (even major mergers) in the analysis. Indeed, we do not find any evidence that the slope of the relation depends on mass or dynamical state. Therefore, we claim that expression (30) can be regarded as a ’universal’ phase-space density profile with only one free parameter, which is the virial radius of the halo.

We now show in Figure 6 that angular momentum is indeed dominated by random motions. We compare the contribution to the tangential velocity of CDM particles from $`<j>`$ and $`<j^2>`$ over spherical shells. The average bulk rotation velocity of our halos (solid symbols in Figure 6) is about $`0.1v_{\mathrm{vir}}`$, although we find a significant increase near the centre. Never the less, it is worth mentioning that the separation in bulk and random velocity is prone to numerical errors at the innermost regions. The specific angular momentum grows roughly linearly with radius (rigid body rotation), in agreement with previous numerical work (Barnes & Efstathiou, 1987; Bullock et al., 2001; van den Bosch et al., 2002; Chen & Jing, 2002). In any case, the tangential velocity of individual dark matter particles is always dominated by the velocity dispersion, for all haloes, at any radius.

As was discussed in Section 2.3, the orbit eccentricity is an important ingredient of our analytical model, since it measures how deep a CDM particle is able to penetrate within the gravitational potential. The ideal case $`e=0`$ corresponds to a configuration in which all particles would orbit at constant radius with no mixing between concentric shells, while $`e=1`$ (often assumed in spherical infall models) implies that every particle goes through the very centre of mass of the halo. For large values of the eccentricity, the central regions are mainly composed of particles coming from the outer shells, that happen to be near the pericentre of their orbits.

The eccentricity profile of our haloes is plotted in Figure 7. For each object, we compute the spherically-symmetric potential derived from its mass distribution. The pericentric and apocentric radii of every CDM particle are estimated from its position and velocity, and the eccentricity is computed as

$$e=\frac{r_{\mathrm{max}}r_{\mathrm{min}}}{r_{\mathrm{max}}+r_{\mathrm{min}}}.$$
(31)
The general trend is that particle orbits are slightly more radial as we move out to the current turn-around radius of the halo, $`R_{\mathrm{ta}}`$. Moreover, we find a systematic dependence on the dynamical state. As can be seen in Figure 8, major mergers are well described by constant eccentricity up to the virial radius <sup>1</sup><sup>1</sup>1For our dark matter haloes, $`r_{\mathrm{vir}}/R_{\mathrm{ta}}`$ is typically of the order of $`0.20.3`$., while relaxed systems are more consistent with a power-law profile. Minor mergers are somewhat in the middle. The average profile can be fitted by a power law, but the slope is shallower than for relaxed systems.

A least-square fit to the relaxed population yields

$$e(r_{\mathrm{max}})0.8\left(\frac{r_{\mathrm{max}}}{R_{\mathrm{ta}}}\right)^{0.1}$$
(32)
for $`r_{\mathrm{max}}<0.1R_{\mathrm{ta}}`$. We note that our approximation to compute the eccentricity breaks down at larger radii, since the mass distribution (and thus, the gravitational potential) are by no means static beyond the virial radius.

## 4 Model results

We now attempt to reproduce the density profiles of our simulated haloes with the model described in Section 2. In order to keep the number of free parameters to a minimum, we assumed a constant eccentricity, $`e0.5`$, which corresponds to the average over the whole radial range. This value implies that particles are able to sink into the dark matter potential up to one third of their maximum radius, or, equivalently, that all mass within radius $`r`$ comes from shells whose turn-around radius was $`3r`$. The effect of a power law profile will be considered in Section 4.2.

### 4.1 Constant eccentricity

Once the eccentricity is set to $`e=0.5`$, the model has two free parameters, which describe the primordial density peak: its height, $`\nu `$, and the smoothing scale, $`R_\mathrm{f}`$. These parameters set the mass and formation time of the halo. For a given mass, a higher peak on a smaller scale corresponds to an earlier formation time and a steeper density profile near the centre.

For each halo, we performed a $`\chi ^2`$ minimisation of the mass profile. Best-fitting values of $`R_\mathrm{f}`$ and $`\nu `$ are plotted in Figure 9, as well as the areas where $`\sqrt{\chi ^2/(dof)}0.12`$. Most cluster-size haloes can be described as very high-$`\nu `$ peaks, while galaxies seem to have collapsed around less extreme fluctuations. This means that clusters are expected to form *earlier* than galaxies, in the sense that the seeds of their dark matter haloes are already in place at a higher redshift.

An exception to this rule are major mergers. While some of these systems still have an early collapse time, some others have collapsed approximately at the present epoch. The first class corresponds to objects that have not relaxed yet, and their core still corresponds to the old halo. As relaxation completes and substructure is erased, the best-fitting parameters move along the lines of constant mass in the $`\nu R_\mathrm{f}`$ diagram. The smoothing scale rises sharply, and $`\nu `$ decreases accordingly.

As can be seen in Figure 9, there is never the less a certain degeneracy between $`R_\mathrm{f}`$ and $`\nu `$, which prevents a reliable determination of the formation time. The exact value of the best-fitting parameters may vary within the shaded area, depending on the details of the fit. For instance, we tend to get higher peaks on smaller scales as we give more weight to the inner parts of the profile.

This can be understood when we compare the numerical data with the results of our spherical collapse model. We chose to fit the radial range $`0.1<r/r_{\mathrm{vir}}<1`$ in order to test the quality of the extrapolations, both towards the centre and to large radii. Individual circular velocity profiles are shown in Figure 10, together with the best fits provided by our model and NFW formula.

Although the mass distribution is fairly well described in general terms, the central density is usually underestimated by both models. When we fit the innermost parts, we are biased towards more concentrated distributions. Indeed, if only data within $`0.1r_{\mathrm{vir}}`$ are considered, the best-fitting profiles are rather unrealistic.

At large radii, the spherical collapse model reproduces the upturn in the circular velocity, while NFW drops to zero. This is due to the fact that the density in the NFW model vanishes at infinity, while it tends asymptotically to the mean value in the spherical collapse model. Adding a constant density background to the NFW formula is enough to bring the profile in agreement with the numerical data.

We assess the accuracy of both NFW and spherical collapse models in Figure 11, where the quantity

$$\frac{\mathrm{\Delta }M}{M}\frac{M_{\mathrm{model}}M_{\mathrm{data}}}{M_{\mathrm{data}}}$$
(33)
is plotted as a function of radius.

The most important difference between both models is that the accuracy of NFW fit improves significantly for relaxed systems, particularly concerning the extrapolation towards $`r0`$. The fits based on our spherical collapse model are on average less accurate, and the scatter between individual dark matter haloes is a little bit larger. Yet, the uncertainty is always lower than 20 per cent for relaxed objects, and only slightly larger for merging systems, where the prediction of the spherical collapse model is indeed very similar to the NFW fit.

### 4.2 Variable eccentricity

Although constant eccentricity might provide a useful approximation for merging systems, it is evident from Figure 8 that it is not a good description of relaxed haloes, where the eccentricity profile $`e(r)`$ increases significantly from the centre to the turn-around radius. Moreover, the assumption of constant eccentricity leads to systematic differences between the mass distribution predicted by our model and the numerical data, as shown in Figure 11.

Therefore, we would like to investigate the consequences of using a more elaborate prescription for $`e(r)`$. For constant eccentricity, high values of $`e`$ lead to steeper density profiles in the central regions, whereas the outer parts of the halo remain largely unaffected (see Figure 1). We plot in Figure 12 the density profile resulting from our power-law fit (32) for $`\nu =3`$ and $`R_\mathrm{f}=1h^1`$ Mpc. The mass distribution is very similar to that obtained with $`e=0.5`$, but the shape is slightly different. There is a small increase in the density between 10 and 100 $`h^1`$ kpc, but the profile flattens near the centre due to the more circular orbits.

This flattening can be clearly seen on the bottom panel of Figure 12, where we plot the logarithmic slope of the density profile. It is interesting to note that the spherical collapse model predicts a finite density at $`r=0`$ when non-radial motions are included, albeit the size of the ’core’ is extremely small. As was discussed in Section 3.1, a similar trend is shown by the relaxed haloes in our sample.

At our resolution limit, the net effect of using expression (32) is just to increase the density at small radii. Then our model gives a somewhat poorer description of merging systems, but the quality of the fit improves considerably for relaxed haloes. As can be seen in Figure 13, the accuracy of our model is comparable to NFW formula when a realistic prescription is used for the eccentricity of particle orbits.

## 5 Discussion and conclusions

The radial density profile of dark matter haloes has been investigated within the framework of the spherical collapse theory. We have shown that the model described in Section 2 is able to reproduce the mass distribution of realistic CDM haloes. Although the final profile cannot be cast in a simple analytical form, it provides not only a mere phenomenological fit, but a physically-motivated description of the density distribution in terms of the primordial initial conditions.

However, it is not easy to understand how the assumption of spherical symmetry could be able to describe the hierarchical assembly of cosmological structures. Instead of continuous infall of spherical shells, the formation of CDM haloes observed in numerical experiments takes place in a discrete and anisotropic way. Most of the matter is accreted in clumps, along the preferred directions set by the filamentary large scale structure.

None the less, the very complicated coalescence process looks very regular in energy space (Zaroubi et al., 1996). Moreover, Moore et al. (1999) have shown that the final density profile is not very sensitive to the details of the merging history, by comparing the mass distribution of a galaxy cluster halo with a re-simulation in which the power spectrum was truncated at $`4`$ Mpc scales.

We suggest, as a plausible explanation for the success of the spherical infall model, that merging is implicitly taken into account through the smoothing scale $`R_\mathrm{f}`$. Cosmological structures do not form around maxima of primordial density field, but of the *smoothed* density field (Hoffman & Shaham, 1985). Some memory of the initial conditions will be lost during major mergers, particularly at the innermost regions of the resulting halo. Contrary to the common view (see e.g. BBKS), we argue that $`R_\mathrm{f}`$ has a precise physical interpretation; below the mass scale defined by $`R_\mathrm{f}`$, all information about the primordial substructure would have been erased by relaxation processes.

In the outer regions, matter is accreted in a more gentle way. Minor mergers do not significantly alter the dynamical structure of the halo. The mass of infalling clumps is much less than the average over a spherical shell at large radii. Thus, the density profile and the accretion rate are determined by the amount of matter available to the halo, which is ultimately set by the primordial initial conditions.

Never the less, there is still an additional degree of freedom, related to the amount and distribution of angular momentum within the dark matter halo. Angular momentum sets the shape of the density profile at the inner regions. For pure radial orbits, the core is dominated by particles from the outer shells. As the angular momentum increases, these particles remain closer to their maximum radius, resulting in a shallower density profile.

We found that angular momentum is dominated by the tangential component of the velocity dispersion. Random motions are well described by a ’universal’ phase-space density profile over several orders of magnitude in radius. This profile is a power law with slope $`\beta 1.9`$, in agreement with the results of Taylor & Navarro (2001) for galactic-size haloes and Rasia et al. (2003) for clusters. We found that this profile is valid for all objects in our sample (spanning a mass range from $`10^{11}`$ to $`10^{15}`$ M) as long as the virial radius is used as a scale parameter.

Although it would be desirable to establish a link between the initial conditions and the specific angular momentum, we resorted to a phenomenological approach. According to the results of our numerical simulations, the orbit eccentricity can be set to a constant value of $`e=0.5`$ as a first order approximation. However, relaxed haloes are better described by a power-law profile in terms of their turn-around radius. The details of the mass distribution depend only moderately on the prescription assumed for the angular momentum. Therefore, the spherical collapse model provides a valuable tool to predict the structure of virialised dark matter haloes, given the power spectrum of primordial fluctuations. Indeed, we claim that the physical origin of the mass distribution observed at the present day is related to the shape of the primordial density peaks. ’Universal’ profiles with two free parameters arise naturally from Gaussian random peak statistics, since the primordial fluctuations are fully specified by their height, $`\nu `$, and smoothing scale, $`R_\mathrm{f}`$.

## Acknowledgments

The authors would like to thank Andreas Faltenbacher and Yehuda Hoffman for useful discussions. This work has been partially supported by the MCyT (Spain), by the Acciones Integradas Hispano-Alemanas and by Deutscher Akademischer Austauschdienst DAAD (Germany). All simulations have been performed in the supercomputers of the Liebniz-Rechnenzentrum LRZ (Germany).


format: text

text:
# DESY 03-198SFB/CPP-03-56December 2003 Scaling test for Wilson twisted mass QCD

## 1 Introduction

Formulating QCD on a space-time lattice admits a substantial amount of freedom in discretising the continuous derivative with restrictions coming from obeying principles such as gauge invariance, locality and unitarity. Naturally, as long as universality holds, all these formulations should provide consistent results in the limit when the discretisation is removed. The standard Wilson formulation of lattice QCD  is a simple realisation of such a discretised version of QCD and has been used for a long time in lattice simulations. However, it has been realised that this formulation has a number of severe problems: it shows large discretisation effects  that are linear in the lattice spacing $`a`$, it violates chiral symmetry strongly and develops unphysical small eigenvalues of the corresponding lattice Wilson-Dirac operator, even at rather large values of the quark mass.

The problem of discretisation effects can be overcome following the Symanzik improvement program , introducing the well known clover term  to get full, on-shell O$`(a)`$ improved results, if the improvement is performed non-perturbatively . In this approach then also the chiral properties are improved, although chiral symmetry breaking and O$`(a^2)`$ lattice spacing effects are still left in the theory and have to be extrapolated away. Despite the fact that non-perturbatively improved Wilson fermions clearly diminish discretisation errors, they, unfortunately, show the same – if not worse – problem of the appearance of small eigenvalues of the lattice Dirac operator .

In order to solve the problem of small eigenvalues, it has been proposed in  to use the so-called twisted mass formulation of QCD. In this approach the mass term in the Dirac operator is chirally twisted , see also . When using such a chirally twisted lattice action in combination with the clover term and a non-perturbatively tuned value of the improvement coefficient, the theory is O$`(a)`$ improved and the corresponding lattice Dirac operator is safe against developing small eigenvalues.

Recently, it has been realised  that all the above properties of non-perturbatively improved, twisted mass QCD can also be obtained when the clover term is completely omitted. If a special value of the twisting angle is chosen one theoretically obtains O$`(a)`$ improved results without adding improvement terms. At the same time, the lattice Dirac operator is still protected against the appearance of small eigenvalues by construction. This curious observation receives a special importance for simulations with dynamical quarks: In  it was found that dynamical fermion simulations with non-perturbatively improved Wilson fermions show signs of first order phase transitions which render simulations very difficult and induce large cut-off effects.

Although by varying the form of the gauge actions  the first order phase transition seems to vanish, it is unclear, whether eventually such phenomena will reappear. The problem of the presence of the first order phase transition may be related to the fact that the clover term in the fermion action generates an adjoint gauge action. Therefore, if O$`(a)`$ improvement can be achieved without the clover term, as anticipated in , the problems connected to such phase transitions should be completely eliminated. As a consequence, the potential of twisted mass fermions in general may be very large. Since the small eigenvalues are regulated by the twisted mass parameter, simulations at much lower quark masses than used today could be performed with promising, but hard to estimate, advantages to explore the chiral limit of lattice QCD.

In this paper we provide a first test of the conjecture of O$`(a)`$ improvement of Wilson twisted mass QCD in quenched numerical simulations. For this purpose we performed a scaling test of the vector meson mass $`m_\mathrm{V}`$ and the pseudoscalar decay constant $`F_{\mathrm{PS}}`$ at a fixed value of the physical pseudoscalar mass $`m_{\mathrm{PS}}`$. We compare the results with those that have been obtained for standard Wilson fermions, see  and references therein, and non-perturbatively improved clover fermions .

## 2 Twisted mass QCD with Wilson quarks

In twisted mass QCD (tmQCD) as formulated in  the twisted mass action in the continuum reads as follows:

$$S_F[\psi ,\overline{\psi }]=d^4x\overline{\psi }(D_\mu \gamma _\mu +m_0+i\mu _q\gamma _5\tau ^3)\psi ,$$
(2-1)
where $`D_\mu `$ denotes the usual covariant derivative, $`m_0`$ is the standard bare quark mass, $`\tau ^3`$ is the third Pauli matrix acting in flavour space and $`\mu _q`$ is the twisted mass parameter, also referred to as the twisted mass.

An axial transformation,

$$\psi ^{}=\mathrm{exp}(i\omega \gamma _5\tau ^3/2)\psi ,\overline{\psi }^{}=\overline{\psi }\mathrm{exp}(i\omega \gamma _5\tau ^3/2),$$
(2-2)
with a real rotation angle $`\omega `$ leaves the form of the action invariant and merely changes the mass parameters into $`m_0^{}`$ and $`\mu _q^{}`$,

$`m_0^{}`$ $`=`$ $`m_0\mathrm{cos}(\omega )+\mu _q\mathrm{sin}(\omega )`$
$`\mu _q^{}`$ $`=`$ $`m_0\mathrm{sin}(\omega )+\mu _q\mathrm{cos}(\omega ).`$ (2-3)
The standard action ($`\mu _q^{}=0`$) is obtained by setting $`\mathrm{tan}\omega =\mu _q/m_0`$. Note that $`\tau ^3`$ is traceless and therefore the transformation (2-2) does not couple to the fermion determinant anomaly.

In order to have an O$`(a)`$ improved twisted mass lattice action, it appears to be natural to discretise the Dirac operator adding appropriate improvement terms to the standard Wilson-Dirac operator. Indeed, it has been shown in  that by adding the usual clover term, full on-shell O$`(a)`$ improvement can be obtained.

In  it has been realised later that using simply the standard lattice Dirac operator $`D_W`$ one can obtain O$`(a)`$ improved physical observables without adding improvement terms. More precisely, it is possible to obtain O$`(a)`$ improved lattice results by employing only the standard massless Wilson-Dirac operator,

$$D_W=\frac{1}{2}\{\gamma _\mu (_\mu ^{}+_\mu )ar_\mu ^{}_\mu \}$$
(2-4)
under the condition that one averages over physical observables that are obtained from simulations at positive and negative values of the Wilson parameter $`r`$. A less general, but similar suggestion was made by the authors of refs. . Instead of using positive and negative values of $`r`$, quark masses with different signs may be used: the bare quark mass $`m_0`$ can be written as

$$m_0=m_c(r)+m_q,\text{with}m_c(r)=m_c(r),$$
(2-5)
where $`m_c(r)`$ is the critical quark mass, and $`m_q`$ is the subtracted bare quark mass. Averaging physical observables obtained from simulations at positive and negative subtracted bare quark masses, again O$`(a)`$-improvement is obtained.

Let us shortly sketch the arguments leading to this surprising result. One first has to observe that with

$$_5:\{\begin{array}{cc}\psi \hfill & \psi ^{}=\gamma _5\psi \hfill \\ \overline{\psi }\hfill & \overline{\psi }^{}=\overline{\psi }\gamma _5\hfill \end{array}$$
(2-6)
the following combined transformation

$$_5^{\mathrm{𝐬𝐩}}_5\times [rr]\times [m_qm_q]$$
(2-7)
is a so called spurionic symmetry of the ordinary Wilson action. Another symmetry of the Wilson (and Wilson tmQCD) action is $`_5\times 𝒟_d`$. In the continuum, the transformation $`𝒟_d`$ has the effect of changing the sign of all the space-time coordinates and multiplies each local term $`_i`$ in the Lagrangian density by the factor $`(1)^{d[_i]}`$, where $`d[_i]`$ is the naive dimension of $`_i`$. The lattice version of this transformation is more involved and we refer to ref.  for details.

Taking now the parity properties of multiplicatively renormalisable operators under $`_5^{\mathrm{𝐬𝐩}}`$ and $`_5\times 𝒟_d`$ into account, one can show – using the Symanzik expansion – that one gets O$`(a)`$ improvement when averaging over two simulations with positive and negative Wilson coefficient $`r`$ (Wilson average (WA)). In addition, from the spurionic symmetry $`_5^{\mathrm{𝐬𝐩}}`$ of the Wilson (and Wilson tm) action one can obtain O$`(a)`$ improved physical observables when averaging, at a fixed value of $`r`$, over two simulations with positive and negative $`m_q`$, as defined in eq. (2-5) (mass average (MA)), taking into account this time the $`_5`$-parity of the operators. Studying the chiral properties of the scalar condensate with Wilson fermions, a similar suggestion was made by the authors in ref. . In the special case of choosing $`\omega =\pm \pi /2`$, such an averaging procedure is done automatically. A change of the sign of $`r`$ is equivalent to $`\omega \omega +\pi `$. Hence, for $`\omega =\pm \pi /2`$ all the quantities that are even under $`\omega \omega `$ are automatically improved without any averaging procedure.

It is the aim of this paper to check this conjecture in practical simulations by performing a scaling test for the vector meson mass $`m_\mathrm{V}`$ and the pseudoscalar decay constant $`F_{\mathrm{PS}}`$ at $`\omega =\pi /2`$. The main goal is to test whether the results for $`m_\mathrm{V}`$ and $`F_{\mathrm{PS}}`$ are consistent with the anticipated leading O$`(a^2)`$ behaviour and that the linear $`a`$ dependence is indeed cancelled. In addition, it is an interesting and important question, what the size of the remaining lattice spacing effects arising in O$`(a^2)`$ will be.

Let us list a few properties of the composite fields in the twisted mass formulation before going to the numerical results. Due to the transformation rule (2-2) one also has to transform the composite fields defined in the usual way,

$`S^0(x)`$ $`=`$ $`\overline{\psi }(x)\psi (x),P^\alpha (x)=\overline{\psi }(x)\gamma _5{\displaystyle \frac{\tau ^\alpha }{2}}\psi (x),`$
$`A_\mu ^\alpha (x)`$ $`=`$ $`\overline{\psi }(x)\gamma _\mu \gamma _5{\displaystyle \frac{\tau ^\alpha }{2}}\psi (x),V_\mu ^\alpha (x)=\overline{\psi }(x)\gamma _\mu {\displaystyle \frac{\tau ^\alpha }{2}}\psi (x).`$ (2-8)
As an example we give here the relations for the axial and vector currents in the “physical basis” (primed quantities) and “the twisted basis” (unprimed quantities),

$$A_\mu ^\alpha =\{\begin{array}{cc}\mathrm{cos}(\omega )A_\mu ^\alpha +ϵ^{3\alpha \beta }\mathrm{sin}(\omega )V_\mu ^\beta \hfill & \text{(}\alpha =1,2\text{)},\hfill \\ A_\mu ^3\hfill & \text{(}\alpha =3\text{)},\hfill \end{array}$$
(2-9)
$$V_\mu ^\alpha =\{\begin{array}{cc}\mathrm{cos}(\omega )V_\mu ^\alpha +ϵ^{3\alpha \beta }\mathrm{sin}(\omega )A_\mu ^\beta \hfill & \text{(}\alpha =1,2\text{)},\hfill \\ V_\mu ^3\hfill & \text{(}\alpha =3\text{)}.\hfill \end{array}$$
(2-10)
Note that in eq. (2-9) for $`\alpha =1,2`$ and $`\omega =\pi /2`$ the role of the axial and vector currents are just interchanged. Of particular interest is the PCVC relation, which takes the following form in the twisted basis:

$$_\mu ^{}V_\mu ^\alpha =2\mu _qϵ^{3\alpha \beta }P^\beta ,$$
(2-11)
where $`_\mu ^{}`$ is the usual backward derivative. Through a vector variation of the action one obtains the point-split vector current as defined in . This current is protected against renormalisation and using the point-split vector current, the PCVC relation is an exact lattice identity. This implies that $`Z_P=Z_\mu ^1`$, where $`Z_\mu `$ is the renormalisation constant for the twisted mass $`\mu _q`$. This will become important in the extraction of the pseudoscalar decay constant $`F_{\mathrm{PS}}`$ as described below.

## 3 Numerical Tests

In this section, we describe our numerical results for testing the scaling behaviour of Wilson tmQCD in the quenched approximation. We started our investigation by performing a comparative benchmark study of different solvers for obtaining the quark propagator. We found the CGS algorithm  to be superior to the BiCGstab and the CG algorithms. We therefore used the CGS algorithm throughout this work. Gauge field configurations were generated by standard heat-bath and over-relaxation techniques.

### 3.1 Mass average

In order to test the predictions of ref. , we started with the mass average procedure. To this end, we selected a value of $`\beta =5.85`$, set the Wilson parameter $`r=1`$ and performed simulations on $`12^3\times 24`$ lattices at positive and negative values of $`m_q=1/2(1/\kappa 1/\kappa _c)`$. While for $`m_q=+0.02725`$, the propagator computations went smoothly, for $`m_q=0.02725`$ the computation of the quark propagator was exceedingly expensive. The reason for this behaviour can be traced back to the spectrum of the Wilson tmQCD operator as can be seen in fig. 1.

Comparing fig. 1(a) with fig. 1(b), in the case of negative $`m_q`$ one has to deal with extremely small eigenvalues of the operator $`(D_W+m_0)^{}(D_W+m_0)`$. Clearly, these very small low-lying eigenvalues lead to a poor convergence of the solver employed and hence to very costly simulations. Projecting out these small modes does also not help in this respect since the computation of the eigenvalues is again costly.

We therefore proceeded to the “self-averaging” case of choosing $`\omega =\pi /2`$, for which it has been shown in  that one gets O$`(a)`$ improvement even without the need of any averaging for all quantities that are even under $`\omega \omega `$. In our practical implementation we have used the twisted basis. Hence, the choice $`\omega =\pi /2`$ corresponds to set $`m_q=0`$ and $`\mu _q0`$. In this situation, the corresponding Wilson-Dirac operator is protected against small eigenvalues and we do not expect difficulties with the simulations.

### 3.2 Scaling of the vector meson mass

In order to verify the prediction of ref.  we computed the vector meson mass $`m_\mathrm{V}`$ and the pseudoscalar decay constant $`F_{\mathrm{PS}}`$ for the following values of $`\beta `$: 5.85, 6.0, 6.1, 6.2. We used periodic boundary conditions. The corresponding lattice volumes were $`14^3\times 28`$, $`16^3\times 32`$, $`20^3\times 40`$ and $`24^3\times 48`$, respectively.

For our simulation in the twisted basis at $`m_0=m_c`$ we had to determine the critical hopping parameter $`\kappa _c`$ for each value of $`\beta `$. At all the $`\beta `$ values of our simulations we made our own determination of the value of $`\kappa _c`$ from the intercept in $`\kappa `$ at zero pion mass. The values of $`\kappa _c`$ are given in table 1. Note that these critical values of $`\kappa `$ have an intrinsic uncertainty of O$`(a)`$. This is, however, sufficient for obtaining fully O$`(a)`$ improved results in Wilson tmQCD .

In order to fix the physical situation in our scaling test, we kept $`r_0m_{\mathrm{PS}}`$ fixed for all values of $`\beta `$. For this purpose we determined the value of $`\mu _q`$ to fix $`r_0m_{\mathrm{PS}}=1.79`$. The corresponding values of $`\mu _q`$ for each value of $`\beta `$ and all our simulation parameters are given in table 1.

We computed the standard 2-point correlation functions at zero momentum for the pseudoscalar and axial operators (which in the twisted basis at $`\omega =\pi /2`$ gives the correct operator to extract the vector meson mass),

$$f_P^\alpha (t)=\underset{\stackrel{}{x}}{}P^\alpha (x)P^\alpha (0)$$
(3-1)
$$f_A^\alpha (t)=\frac{1}{3}\underset{i=1}{\overset{3}{}}\underset{\stackrel{}{x}}{}A_i^\alpha (x)A_i^\alpha (0)$$
(3-2)
with $`P^\alpha (x)`$ and $`A^\alpha (x)`$ given in eqs. (2-8) and $`x=(\stackrel{}{x},t)`$. In order to obtain a non-vanishing result, the flavour index has to be the same in these correlation functions and we will choose $`\alpha =1`$ in the following. The pion mass could be extracted easily from the exponential decay of the correlation function $`f_P^\alpha (t)`$. For the vector meson mass, we performed two mass fits for the ground state mass and the first excited state. We checked the stability of the fit by changing the value of $`t_{\mathrm{min}}`$ where the fit started. As a cross check, we also determined the effective ground state mass and found consistent results. All errors were computed by a jackknife analysis. The numerical results at our simulation points are collected in table 2.

In fig. 2 we show our results for the vector meson mass as a function of $`a^2`$ represented by the open circles. In addition, we also show results from non-perturbatively O$`(a)`$ improved Wilson fermions (filled circles) . Finally, we added results for standard Wilson fermion simulations, see  and references therein, as they were available in the literature. We remark that the published data were not always at exactly the same value of $`r_0m_{\mathrm{PS}}`$ that we used for our Wilson tmQCD simulations. In such cases we performed an interpolation to the desired value of $`r_0m_{\mathrm{PS}}`$. The error from this (small) interpolation is negligible for the results presented here.

We performed a simple extrapolation of the Wilson tmQCD and the O$`(a)`$ improved data of the form $`r_0m_\mathrm{V}=r_0m_\mathrm{V}^{\mathrm{cont}}+b(a/r_0)^2`$, with $`r_0m_\mathrm{V}^{\mathrm{cont}}`$ the continuum value of the vector meson mass and $`r_00.5\text{fm}`$. For the pure Wilson results we replaced the quadratic term with a term proportional to $`(a/r_0)`$ in the extrapolation. Let us remark that our data for $`m_\mathrm{V}`$ for Wilson tmQCD has about a factor of four less statistics than the data from O$`(a)`$-improved Wilson fermions, which is reflected in the larger error bars. Nevertheless, it is evident that the Wilson tmQCD results show a very similar scaling behaviour as the O$`(a)`$-improved Wilson fermions. This becomes even clearer when we compare with the unimproved pure Wilson data for $`m_\mathrm{V}`$ that we show as open squares in fig. 2. Here large lattice artefacts are seen and the scaling behaviour is much worse than with Wilson tmQCD or O$`(a)`$-improved Wilson fermions. We also note that the data for Wilson tmQCD are rather flat as a function of $`a^2`$ indicating that also higher order lattice spacing effects are suppressed. Clearly, it would be desirable to test these promising results in more precise simulations using a much higher statistics.

### 3.3 Scaling test for $`F_{\mathrm{PS}}`$

In order to extract the pion decay constant $`F_{\mathrm{PS}}`$ as another quantity to test scaling, we start with the standard definition of $`F_{\mathrm{PS}}`$ (again fixing the flavour index $`\alpha =1`$),

$$0|A_0^1|PS=m_{\mathrm{PS}}F_{\mathrm{PS}}.$$
(3-3)
In the twisted basis, at $`\omega =\pi /2`$, the axial current is related to the vector current by the transformation eq. (2-2), see eqs.(2-9,2-10), and so we can write

$$_\mu 0|V_\mu ^2|\mathrm{PS}=F_{\mathrm{PS}}m_{\mathrm{PS}}^2.$$
(3-4)
Using the vector Ward identity in eq. (2-11), we can finally relate the divergence of the vector current to the pseudoscalar density and obtain

$$F_{\mathrm{PS}}m_{\mathrm{PS}}^2=_\mu 0|V_\mu ^2|\mathrm{PS}=2\mu _q0|P^1|PS.$$
(3-5)
For asymptotic Euclidean times, the pseudoscalar correlation function $`f_P^1`$ assumes the form

$$f_P^1(t)=\frac{|0|P^1|\mathrm{PS}|^2}{2m_{\mathrm{PS}}}\left(e^{m_{\mathrm{PS}}t}+e^{m_{\mathrm{PS}}(Tt)}\right),atT.$$
(3-6)
Thus, by fitting the pseudoscalar correlation function for large time separations, we can obtain $`m_{\mathrm{PS}}`$ and the amplitude $`|0|P^1|\mathrm{PS}|^2/m_{\mathrm{PS}}`$ from which we then compute the desired matrix element $`|0|P^1|\mathrm{PS}|`$. Hence, we have all necessary ingredients to determine $`F_{\mathrm{PS}}`$ from eq. (3-5), without the need of any renormalisation factor, since $`Z_P=Z_\mu ^1`$, as we have mentioned in the previous section. The error estimate of the so computed value of $`F_{\mathrm{PS}}`$ is performed with a jackknife procedure. We note that the discussion of how to get $`F_{\mathrm{PS}}`$ with tmQCD resembles very closely the strategy one would follow in the continuum or with overlap fermions.

In fig. 3 we show our results for $`r_0F_{\mathrm{PS}}`$ at a fixed value of $`r_0m_{\mathrm{PS}}=1.79`$. We also add results from O$`(a)`$-improved Wilson fermions  in the figure. We performed an extrapolation of the Wilson tmQCD and the O$`(a)`$-improved data of the form $`r_0F_{\mathrm{PS}}=r_0F_{\mathrm{PS}}^{\mathrm{cont}}+b^{}(a/r_0)^2`$, with $`r_0F_{\mathrm{PS}}^{\mathrm{cont}}`$ the continuum value of the pion decay constant in units of $`r_00.5\text{fm}`$. While the O$`(a)`$-improved Wilson results do show cut-off effects, presumably due to the particular value of the improvement coefficient $`c_A`$ employed in these calculations, the results for Wilson tmQCD are essentially flat, indicating that also the higher order lattice spacing effects are substantially smaller.

## 4 Conclusion and Outlook

In this paper we tested the striking idea of using the standard Wilson-Dirac operator with a twisted mass and twist angle $`\omega =\pi /2`$ to obtain full O$`(a)`$ improvement for physical quantities that are even under $`\omega \omega `$. The results of our – quenched – study for the vector meson mass $`m_\mathrm{V}`$ and the pion decay constant $`F_{\mathrm{PS}}`$ are very encouraging. It seems that in this setup the lattice spacing effects are substantially reduced with respect to standard Wilson fermions and consistent with vanishing O$`(a)`$ discretisation errors. At the same time it seems that also lattice spacing effects that come in higher orders in $`a`$ are small, a result that could not be anticipated before.

Thus it seems that Wilson twisted mass QCD has the potential to solve many problems at once:

* It automatically reduces lattice spacing effects, consistent with O$`(a)`$ improvement.
* It protects against very small low-lying eigenvalues of the lattice Dirac operator.
* It can avoid unwanted phase transitions as appeared in dynamical simulations with clover improved Wilson fermions and Wilson gauge action.

In the light of this, dynamical simulations with twisted mass QCD appear to be very promising to approach the physical point at realistic values of quark masses. In some way, the twisted mass plays a similar role as an infrared cut-off as the quark mass in the staggered fermion or overlap approach does. The very good performance of staggered dynamical simulations gives hope that also in twisted mass QCD the simulations can be accelerated as compared to standard, improved or unimproved, Wilson fermions. Naturally, with Wilson twisted mass QCD the problem of unphysical taste degrees of freedom is completely avoided, which, to us, is a clear advantage of the tmQCD idea.

## Acknowledgements

We thank Michele Della Morte, Roberto Frezzotti and Stefan Sint for many valuable discussions. This work was supported in part by the European Union Improving Human Potential Programme under contracts HPRN-CT-2002-00311 (EURIDICE) and by the DFG through the SFB/TR9-03 (Aachen-Berlin-Karlsruhe). We are thankful to the John von Neumann-Institute for Computing for providing the necessary computer resources to perform this project.