text:
# Specific heat of Ce0.8La0.2Al3 in magnetic fields: a test of the anisotropic Kondo picture

\[

## Abstract

The specific heat $`C`$ of Ce0.8La0.2Al3 has been measured as a function of temperature $`T`$ in magnetic fields up to 14 T. A large peak in $`C`$ at 2.3 K has recently been ascribed to an anisotropic Kondo effect in this compound. A 14-T field depresses the temperature of the peak by only 0.2 K, but strongly reduces its height. The corresponding peak in $`C/T`$ shifts from 2.1 K at zero field to 1.7 K at 14 T. The extrapolated specific heat coefficient $`\gamma =lim_{T0}C/T`$ increases with field over the range studied. We show that these trends are inconsistent with the anisotropic Kondo model.

\]

CeAl3 occupies a particularly important position in the history of heavy fermions. The first report of its unusual low-temperature specific heat a quarter of a century ago led to enormous interest in this and similar systems based on 4$`f`$ or 5$`f`$ elements. For years, CeAl3 was considered a canonical heavy fermion system, and it greatly influenced theoretical work in the field. Indeed, CeAl3 was viewed as a realization of a “standard model” based on the Kondo effect and Fermi-liquid theory. The hallmark properties of this compound are its specific heat and its electrical resistivity. The low-temperature specific heat $`C(T)`$ is greatly enhanced over that for a conventional metal, with a linear coefficient $`\gamma =lim_{T0}C/T1250`$ mJ/K2Ce mol. The resistivity above 10 K is described to very high very accuracy by a theory for a Kondo impurity in crystalline electric fields. Below 300 mK the resistivity has a Fermi-liquid form $`\rho =\rho _0+AT^2`$, with a strongly enhanced $`A`$ coefficient. The ratio $`A/\gamma ^2`$ is about $`10^5\mathrm{\Omega }`$ cm K2mol2J-2, a value close to that theoretically predicted for nonmagnetic Kondo lattices, and experimentally observed in many other heavy fermion compounds.

The description of CeAl3 in terms of a nonmagnetic Kondo lattice and a heavy Fermi liquid has been challenged by microscopic measurements such as muon spin resonance ($`\mu `$SR) and nuclear magnetic resonance (NMR). According to these measurements, either short-range magnetic order or strong antiferromagnetic correlations exist below 2 K. In addition, the specific heat itself has an unexplained feature: a maximum in $`C/T`$ near 0.4 K. A similar maximum is found in another heavy fermion system, CeCu2Si2, around the same temperature. In both compounds, this feature was initially attributed to coherence in the Kondo lattice. However, extensive studies of CeCu2Si2 gave rise to an alternative explanation based on weak magnetic ordering of heavy quasiparticles. A previous alloying study of CeAl3 has similarly pointed to a magnetic origin for the 0.4 K anomaly. When La is partially substituted for Ce in Ce1-xLaxAl3, this weak feature, observable in $`C/T`$ but not in $`C`$ for $`x=0`$, gradually evolves for $`x0.05`$ into a large peak in both $`C`$ and $`C/T`$. The highest-La-content alloy investigated in Ref. , Ce0.8La0.2Al3, has a pronounced maximum in $`C`$ near 2.3 K and a corresponding peak in the susceptibility at 2.5 K, reminiscent of an antiferromagnetic transition. The smooth and monotonic increase with $`x`$ of the temperature position and magnitude of the anomaly suggests that this feature has a common origin in pure and La-doped CeAl3. The apparent enhancement of the magnetic character of the anomaly upon La doping is consistent with Doniach’s Kondo necklace model, since doping increases the lattice constants, and therefore decreases the hybridization between $`f`$ and ligand states. However, an interpretation based on this model is somewhat undermined by the fact that anomalies in $`C`$, similarly pronounced to those produced by La substitution, can be induced by replacing Al atoms with either larger or smaller atoms.

Recent neutron scattering and $`\mu `$SR studies by Goremychkin et al.  on Ce0.8La0.2Al3 revealed the absence of magnetic Bragg peaks, and estimated the upper limit of any possible ordered moment to be $`0.05\mu _\text{B}`$. The response function deduced from time-of-flight measurements changes from a quasi-elastic form to an inelastic form around 3 K, the temperature range where features develop in the specific heat and the magnetic susceptibility. This result was attributed to weakly dissipative dynamics consistent with the anisotropic Kondo model (AKM). $`\mu `$SR spectra showed Lorentzian damping, with a temperature-dependent damping rate that diverges also around 3 K. The divergence was attributed to the development of static magnetic correlations, indicating the possibility of magnetic order of small moments, as seen in other heavy fermion systems.

In order to investigate further the applicability of the AKM to Ce0.8La0.2Al3, and to search for any contribution to the specific heat from small-moment ordering, we have studied the effect of magnetic fields up to 14 T on the linear coefficient $`\gamma `$ and on the temperatures $`T_M`$ and $`T_m`$ of the maxima in $`C`$ and $`C/T`$, respectively. We compare our data with numerical results for the specific heat of the AKM. The experimental data show qualitatively different trends from the model calculations, thereby casting considerable doubt on the validity of this theoretical description of Ce0.8La0.2Al3.

Our measurements used a polycrystal of Ce0.8La0.2Al3 from a previous alloying study. The sample was prepared in an arc furnace under an argon atmosphere, and annealed at 830$`^{}`$C for 3 weeks. Magnetic susceptibility and x-ray diffraction measurements revealed no sign of the secondary phases CeAl2 and Ce3Al11.

Figure 1 shows the specific heat of this alloy in fields of 0, 5, 10, and 14 T. The phonon contribution has been subtracted using the specific heat of LaAl3, and the remainder has been normalized to a mole of Ce. The same data are plotted as $`C/T`$ vs $`T`$ in Fig. 2. The main effect of the field is a strong reduction in the magnitude of the anomalies in $`C`$ and $`C/T`$. Also striking is the very weak field dependence of the temperature position of the anomalies. A pronounced peak in $`C`$ located at $`T_M2.3`$ K for $`H=0`$ is replaced by a shoulder near 2.1 K for $`H=14`$ T. The peak in $`C/T`$ also shifts slowly with field, $`T_m`$ decreasing from 2.1 K at 0 T to 1.7 K at 14 T (see Fig. 3). Note that the difference between $`T_M`$ and $`T_m`$ grows with applied field. A difference of the same order has been observed in zero field for Ce1-xLaxAl3 alloys with $`x<0.2`$, where $`T_mT_m`$ grows as $`x`$ becomes smaller. In this respect, an increase in the magnetic field has a similar effect to a decrease in $`x`$.

Another important result is an increase with field of $`C/T`$ values at low temperatures (below 1 K), signaling a partial restoration of the heavy fermion state present in pure CeAl3. It may be that the large nuclear moments of Al contribute to enhance $`C/T`$ at the lowest temperatures and the largest fields. Indeed, the 14-tesla $`C/T`$ data display a low-temperature tail which might be due to a nuclear hyperfine contribution $`\mathrm{\Delta }C/T1/T^3`$. None of the curves at lower fields show a similar upturn. Therefore, the linear specific heat coefficient $`\gamma `$ was extracted from a linear fit to $`C/T`$ vs $`T^2`$ below 1 K, except for the 14-T data, where $`\gamma `$ was determined from the slope of $`CT^2`$ vs $`T^3`$ below 1 K. As may be seen in Fig. 3, $`\gamma `$ seems to saturate in the range $`H10`$ T. (The error bars for $`\gamma `$ combine experimental and regression uncertainties.)

It is worth noting that $`C/T`$ for Ce0.8La0.2Al3 at 14 T and $`C/T`$ for CeAl3 in zero field coincide above 4 K to within the accuracy of the measurement. This is demonstrated in Fig. 4, which also includes the corresponding curve for Ce0.8La0.2Al3 at $`H=0`$. Since $`C/T`$ for the pure compound is only weakly field dependent above 4 K (for fields $`10`$ T), we can claim that the high-field ($`H14`$ T) specific heats for these two alloys converge in this temperature regime.

We now attempt to analyze our magnetic-field data in terms of the anisotropic Kondo model (AKM) for a single magnetic impurity. The model assumes an exchange interaction $`J_zS_zs_z+J_{}(S_xs_x+S_ys_y)`$ between the impurity spin $`𝐒`$ and the net conduction-electron spin $`𝐬`$ at the impurity site. Goremychkin et al. have proposed the AKM as a description for the thermodynamic properties of both Ce0.8La0.2Al3 and CeAl3. A strong dependence on field orientation in the magnetic susceptibility of CeAl3 single crystals is suggestive of anisotropic behavior corresponding to $`J_zJ_{}>0`$, with the magnetic $`z`$ direction being the crystallographic $`c`$ axis.

The AKM is known to be equivalent in the limit of low energies to a number of other models. In recent years a mapping onto the spin-boson model with Ohmic dissipation has been exploited  to deduce physical properties of dissipative two-level systems from numerical calculations performed for the AKM. These studies have shown that under certain conditions, the impurity contribution to the zero-field heat capacity of the AKM exhibits peaks both in $`C`$ and $`C/T`$, qualitatively similar to the data plotted in Figs. 1 and 2. Such peaks are found only for $`\varrho _0J_z1`$, where $`\varrho _0`$ is the conduction-band density of states at the Fermi energy. For a given value of $`\varrho _0J_z`$, the specific heat for all $`J_{}J_z`$ can be collapsed onto a universal scaling curve $`C/\gamma T`$ vs $`T/T_m`$. Moreover, the temperature of the peak in $`C/T`$ is given by $`T_m=\alpha R/\gamma `$ , where $`R`$ is the gas constant and $`\alpha `$ is a function of $`\varrho _0J_z`$ (only).

From the observed peak position $`T_m=2.1`$ K and a value $`\gamma =520`$ J/K2Ce mol extracted as described above, we deduce $`\alpha =\gamma T_m/R=0.13`$ for Ce0.8La0.2Al3 in zero field, in agreement with the estimate of Ref. . We have used $`\alpha `$ and $`T_m`$ as inputs for a numerical renormalization-group calculation of the specific heat of the AKM. Figure 5 shows the predicted effect of applying a uniform magnetic field along the magnetic $`z`$ axis, under the assumption that the impurity and the conduction electrons have $`g`$ factors $`g_i=g_e=2`$. (Changing the $`g`$ factors multiplies the field scale by an overall factor, but does not otherwise affect the results.)

The numerical data exhibit three main trends with increasing field: (1) The anomaly in $`C/T`$ becomes broader and lower. (2) The peak shifts markedly to higher temperatures. (3) $`C/T`$ decreases at all temperatures below the zero-field value of $`T_m`$; the fractional change in $`\gamma `$ is greater than that in the peak height, so that $`\alpha =\gamma T_m/R`$ decreases monotonically with increasing magnetic field, as shown in the legend of Fig. 5.

These numerical results are directly applicable only to single-crystal Ce0.8La0.2Al3 with a magnetic field along the $`c`$ axis. For comparison with our polycrystalline data, one must average over all possible field orientations. The Ising-like crystal-field ground state of Ce3+ in CeAl3  implies that $`g_i=0`$ for the basal-plane components of the magnetic field and, hence, that the specific heat of a polycrystal in field $`H`$ is an equally weighted average of the single-crystal results for all fields between zero and $`H`$. This averaging process preserves trends (1)–(3) above.

Trend (1) accords well with our measurements, but (2) and (3) both run counter to experiment. In Ce0.8La0.2Al3, $`T_m`$ does not rise with increasing field, but instead is weakly depressed, while $`C/T`$ undergoes a small increase at temperatures much below $`T_m`$. In particular, $`\gamma `$ rises sufficiently fast that $`\alpha `$ remains essentially constant up to a 14 T field (see Table I), in contrast to the prediction of the AKM.

The preceding comparisons seem to indicate significant shortcomings in the AKM as a description of Ce0.8La0.2Al3 in magnetic fields. One reason for the inadequacy of the AKM may be the neglect of magnetic correlations around the temperature of the maximum, as identified in the $`\mu `$SR studies of Ref. . It was noted above that the specific heat anomaly is reminiscent of an antiferromagnetic phase transition. The entropy under the peak in $`C/T`$ is a large fraction ($``$50$`\%`$) of $`R\mathrm{ln}2`$, and the linear specific heat coefficient $`\gamma =520`$ mJ/K2Ce mol is less than half that of pure CeAl3 ($`\gamma =1250`$ mJ/K2Ce mol). However, our field data suggest that any magnetic transition associated with the anomaly is rather unusual. We find that $`T_M`$ and $`T_m`$ are depressed in an applied field at a much lower rate than is the Néel temperature in Ce-based heavy fermion systems that order antiferromagnetically. In CeCu5.2Ag0.8, for example, $`T_N`$ is reduced from 0.7 K to 0 K in a field of about 2.5 T. In CePb3, which exhibits unconventional small-moment ordering at 1.1 K, a field of order 10 T depresses $`T_N`$ to zero.

In summary, we have measured the heat capacity of Ce0.8La0.2Al3 as a function of temperature in magnetic fields up to 14 T. The field strongly diminishes the peaks found around 2 K in both $`C`$ and $`C/T`$, but only weakly depresses the peak temperatures. The linear specific heat coefficient increases with field in the direction of the value for pure CeAl3, implying partial restoration of the heavy fermion state suppressed by La doping. We have analyzed our data in terms of the anisotropic Kondo model. The model predicts a shift of the peak in $`C/T`$ to higher temperatures with increasing field, accompanied by a significant reduction in $`C/T`$ at low temperatures. These two trends are at odds with experiment. Our results do not rule out an alternative theoretical picture based on small-moment magnetism. However, the field-insensitivity of the temperature of the heat-capacity peak remains to be understood within this scenario.

This work was supported in part by Department of Energy grant DE-FG02-99ER45748.


format: text

text:
# Memory and superposition in a spin glass

## I Introduction

The non-equilibrium nature of spin glasses and other glassy systems at low temperature implies that a comprehensive study of the dynamics usually requires sophisticated dc or ac relaxation experiments. Dc relaxation is typically used to investigate the aging properties in spin glasses. In such experiments, the sample is first cooled from above the transition temperature $`T_g`$ to a constant temperature. After a waiting time $`t_w`$, a small dc field is applied (or cut off), and the magnetization is recorded versus time. The relaxation occurs in a characteristic way, and shows a clear dependence on $`t_w`$. Relaxation of the low-frequency ac susceptibility is a prime tool to investigate memory effects in glassy magnets. When cooling from above $`T_g`$, a halt at constant temperature $`T_s`$$`<`$$`T_g`$ is made during $`t_s`$, allowing the system to relax towards its equilibrium state at $`T_s`$; both components of the ac-susceptibility then decay in magnitude. This equilibrated state becomes frozen in on further lowering the temperature, and is retrieved on re-heating. The weak low frequency ac-field employed in this kind of experiments does not affect the non-equilibrium processes intrinsic to the sample, but only works as a non-perturbing probe of the system.

In the present paper, an alternative method to investigate non-equilibrium dynamics is used; “simple” remanent magnetization measurements on a Ag(Mn) spin glass sample. The magnetization curves are recorded on heating after specific cooling protocols. The comparison of the different curves yields information on dynamic properties such as aging and memory phenomena.

## II Memory and superposition

The sample used in the experiments is a dilute magnetic alloy: the archetypical three dimensional spin glass Ag(11 at% Mn). This material is chosen as an established reference system to certify the validity of the dc method presented here.
The sample was prepared by melting pure Ag and Mn together at $`T=1000^{}`$C in an evacuated atmosphere. After annealing the sample at $`850^{}`$C for 72 h it was quenched to room temperature. The experiments were performed in a non-commercial low-field SQUID magnetometer . The dc magnetic field is generated by a small superconductive solenoid coil always working in persistent mode during measurements. The background field was less than 1 mOe.

As seen in Fig. 1, the zero field cooled (ZFC) magnetization exhibits a characteristic cusp at T=33K, which also roughly defines the spin glass temperature $`T_g`$ of our sample. Fig. 2 presents the temperature dependence of (a) the thermo-remanent (TRM) and (b) the isothermal remanent (IRM) magnetization. Both were measured on heating in zero field from a low temperature. The different curves have been recorded after cooling the sample using the protocol described below:

1. The sample is cooled in constant field $`H_o`$=0.1 Oe (TRM) or zero field (IRM) from a reference temperature above $`T_g`$ to a stop temperature $`T_s`$ ($`T_s`$ $`<`$ $`T_g`$) , where the sample is kept a stop time $`t_{s1}`$ without changing the field.

2. A field change is made: In the TRM case the field is either cut to zero (zero-field stop, ZFS) or kept constant (field stop, FS); for the IRM, $`H_o`$ is applied. The new magnetic field value is kept during a time $`t_{s2}`$.

3. The field is then shifted back to its initial value and the sample is immediately cooled to a lowest temperature $`T_{min}`$, where

4. the field is cut to zero (TRM) (or kept at zero (IRM)) and the remanent magnetization $`M_{\mathrm{TRM}}(T)`$ (denoted $`\mathrm{TRM}_{\mathrm{ZFS}}`$ and $`\mathrm{TRM}_{\mathrm{FS}}`$ for the case of ZFS and FS, respectively) or $`M_{\mathrm{IRM}}(T)`$ (denoted IRM) is measured on heating the sample at a constant heating rate.

The different curves in Fig. 2 are measured for two different values of $`t_{s1}`$: 0 s and 10000 s, respectively, and in all cases $`t_{s2}=10000`$ s; $`T_s`$=27K; $`T_{min}`$=20 K. The reference curves, without any field stops, are added. Figs. 2 (a) and (b) show that the 0 s $`\mathrm{TRM}_{\mathrm{ZFS}}`$ curve lies significantly below the 10000 s curve whereas the 0 s IRM curve lies well above the 10000 s one. These curves reflect the paramount influence that aging has on the magnetic relaxation in spin glasses. For comparison, Fig. 3 presents parallel dc-relaxation experiments (ZFC, FC and TRM), performed at the temperature $`T_s`$, using the same small magnetic field, for waiting times of $`t_w`$=0 s and $`t_w`$=10000 s; the results for $`t_w`$=1000 s are added to show the continuity of the waiting time dependence. The observation time, $`t`$, is defined as the time elapsed after the field application (cut off).
Both $`\mathrm{TRM}_{\mathrm{FS}}`$ curves in Fig. 2 (a) lie significantly above the reference curve, conveying that a considerable reinforcement of the spin structure occurs during the stop time $`t_{s1}+t_{s2}`$. The effect even overcomes the downward relaxation of the FC magnetization occurring during the stop at 27K (cf. Fig. 3). One notices that the difference between stopping 10000 and 10000+10000=20000 s at $`T_s`$ is significant but comparably small. Provided the experiments are made at low enough field, where there is a linear response to a field change, one can relate the results from different relaxation experiments through the principle of superposition. This implies e.g. that:

$$M_{\mathrm{ZFC}}(t_w,t)=M_{\mathrm{FC}}(0,t_w+t)M_{\mathrm{TRM}}(t_w,t).$$

Fig. 2 (c) shows difference plots, $`\mathrm{\Delta }\mathrm{TRM}_{\mathrm{ZFS}}`$ and $`\mathrm{\Delta }`$IRM, of the $`\mathrm{TRM}_{\mathrm{ZFS}}`$ and IRM curves shown in Figs. 2 (a) and (b) giving a direct measure of the frozen in excess magnetization due to aging. Also plotted in the figure is a difference plot of the $`\mathrm{TRM}_{\mathrm{FS}}`$ curves: $`\mathrm{\Delta }\mathrm{TRM}_{\mathrm{FS}}`$. The inset shows a plot of $`\mathrm{\Delta }\mathrm{TRM}_{\mathrm{ZFS}}`$-$`\mathrm{\Delta }`$IRM and $`\mathrm{\Delta }\mathrm{TRM}_{\mathrm{FS}}`$, the two curves look very similar and suggest an extended validity of the superposition principle to apply also to the temperature dependence of the frozen in excess magnetization of a spin glass. (IRM reflects $`M_{\mathrm{ZFC}}`$, $`\mathrm{TRM}_{\mathrm{ZFS}}`$ reflects $`M_{\mathrm{TRM}}`$, and $`M_{\mathrm{FC}}`$ corresponds to $`\mathrm{\Delta }\mathrm{TRM}_{\mathrm{FS}}`$.) The applicability of the “ordinary” superposition principal for the relaxation of our spin glass is illustrated in Fig. 3; the insert shows the sum of the ZFC and TRM relaxations obtained for $`t_w`$=0 s and the FC relaxation for the same waiting time. The curves are plotted without markers to be able to visually distinguish that there are two of them; as expected, they are closely equal to each other.

We have also studied the effect of simple field stops in the ZFC and FC magnetizations. Fig. 4 shows the ZFC, FC and TRM reference curves already presented in Fig. 1, as well as the same curves recorded on re-heating after a stop of duration $`t_s`$=104 s in constant field at $`T_s`$=27K. In the TRM and FC case, the field is thus kept to its $`H`$=0.1 Oe value, while for the ZFC, it remains zero. The FC curve is only weakly affected by the stop, whereas the ZFC and TRM curves are considerably affected. The curves reflect the reinforcement of the spin structure that occurs during the stop at constant temperature and that this reinforcement sustains when the temperature is recovered on re-heating. The insert of Fig. 4 shows the difference between the FS and reference curves for ZFC, TRM and FC, the excess magnetization gained due to the FS in TRM + FC corresponds to the magnetization “lost” in the ZFC. The extended validity of the principal of superposition is further supported by this apparent agreement.

Looking at the bump and dip in the TRM and ZFC magnetization curves of Fig. 4 and the derived excess magnetization in the inset, one observes that the influence of the stop at constant temperature during cooling is limited to a restricted temperature range around $`T_s`$. The width of this region may be assigned to the existence of an overlap between the spin configuration attained at $`T_s`$ and the corresponding state at a very neighboring temperature ($`T+\mathrm{\Delta }T`$). The two concepts that explain the width are then a chaotic nature of the spin glass equilibrium configuration and an overlap on short length scales between the equilibrium configurations at $`T`$ and $`T+\mathrm{\Delta }T`$. These concepts have been elaborated earlier . Further information on these matters can be obtained also by extending the current dc-method to include two stops when cooling, first at $`T_s`$ and then at $`T_s^{}`$$`<`$$`T_s`$. The results of this cooling protocol emphasize the chaotic nature of the equilibrium states at the two temperatures: the equilibration at $`T_s`$ has no influence on the non-equilibrium processes occurring only a few kelvins below, outside the region of overlap.
The actual spin configuration of a spin glass is controlled by its thermal history. In this paper we have shown that it is possible to extract substantial information on complex non-equilibrium phenomena linked to the evolving spin configuration by simply recording the temperature dependence of the remanent magnetization after employing certain cooling protocols. Introducing an excess magnetization in the spin glass by intermittently applying or removing a magnetic field during a stop at constant temperature it is possible to study memory phenomena and the influence of aging on the relaxation of the magnetisation. The magnitude of the excess magnetization is governed by the duration of the field cycle and the wait time at constant temperature before the magnetic field is changed. The logarithmic nature of the relaxation and the aging processes can be explored by systematically altering the duration of the field cycle and/or the wait time. An extended validity of the principal of superposition to apply to frozen in relaxation processes that are thermally activated on heating is disclosed by combining results from TRM, ZFC and FC experiments.

A note on the experimental procedure: the relaxation of spin glasses extends to eternal times and has a logarithmic character and the relaxation processes are thermally activated. These properties require that the cooling or heating rates in the experiments are controlled and kept the same in all different measurements to achieve comparable experimental data. The heating (and cooling) rate determines, due to the overlap between equilibrium states at different temperatures, an effective age of the system and defines an effective observation time (of order 10 s in our experiments) in the measurement of the remanent magnetization.

## III Conclusion

The present results from the here introduced dc method have been obtained on an archetypical 3D spin-glass. Corresponding remanent magnetization studies can advantageously be performed even using standard SQUID magnetometers. For example, the slow dynamics of the less conventional La1-xSrxCoO3 glassy system has been successfully investigated in this way, exposing characteristic memory and superposition features.

###### Acknowledgements.

Financial support from the Swedish Natural Science Research Council (NFR) is acknowledged. Marie Vennström is acknowledged for helping us fabricate the Ag(Mn) sample. D. N. H. Nam is thankful to SAREC/Sida and ISP.


text:
# A simulated 𝜏CDM cosmology cluster catalogue: the NFW profile and the temperature-mass scaling relations

## 1 Introduction

Clusters of galaxies are used to constrain cosmological parameters as they are the largest gravitationally bound systems and preserve imprints of the evolution of the universe. One of the most reliable predictors for the mass-function of collapsed objects is that of Press & Schechter (1974). Although based on naïve assumptions, it has been widely used because of its simple form and, more importantly, because the mass function is in excellent agreement with many N-body simulations (e.g. Lacey & Cole 1993; but see Gross et al. 1998; Governato et al. 1999; Jenkins et al. 2000).

In practice, cluster masses are not easily measured. For this reason, other properties such as X-ray luminosity or temperature are used to estimate the mass. The former of these is notoriously unreliable as it depends strongly upon emission from the core of the cluster; hence we will concentrate on the latter in this paper. A simple application of the virial theorem suggests that the mean temperature of an object that has collapsed to the virial radius, $`r_{\mathrm{vir}}`$, follows the scaling law

$$T_{\mathrm{vir}}\frac{M_{\mathrm{vir}}}{r_{\mathrm{vir}}}M_{\mathrm{vir}}^{2/3}(1+z_{\mathrm{form}}),$$
(1)
where $`M_{\mathrm{vir}}`$ is the mass within $`r_{\mathrm{vir}}`$ and $`z_{\mathrm{form}}`$ is the formation redshift. In practice, clusters are assumed to have formed at the redshift that they are observed, and the virial radius is defined as the radius of a sphere enclosing an overdensity of 180 (in critical-density cosmologies) relative to the mean density at that redshift.

Equation 1 is often assumed to be a rigorous theoretical prediction, but this is not the case. It assumes that the constant of proportionality in Equation 1 is the same for low- and high-mass clusters, but there is no reason why this need be true. For a power-law density fluctuation spectrum of dark-matter particles, the theory of self-similarity tells us that the population of clusters at one redshift is a scaled version of the population at another redshift, but that is not the same as saying that coeval low- and high-mass clusters are similar in form. Even if the clusters do have similar morphologies, we may still get different proportionality constants if we measure their properties within the virial radius, as this will correspond to different multiples of the half-mass radius in each case.

Another deficiency of Equation 1 is that it predicts only the virial temperature—there is empirical evidence that the observed (emission-weighted) X-ray temperature of the intracluster medium, $`T_X`$, is greater than the virial temperature (e.g. Edge & Stewart 1991a,b; Bahcall & Lubin 1994). Were these to be in constant ratio, then the scaling relations would be preserved, however this is unlikely to be the case. Physical processes such as shock-heating and radiative cooling act only on the gaseous component of clusters and not on the dark matter. Furthermore, these two processes have different dependencies upon density and so will break the expected self-similarity.

Hjorth, Oukbir & van Kampen (1998), in a sample of 8 clusters whose masses are inferred from weak and strong gravitational lensing and temperatures are measured from ASCA, confirm that (very roughly) $`T_XM/r`$. Most observers, however, report the direct relation between $`T_X`$ and $`M`$ without an explicit dependence upon $`r`$ and so that is what we will do here. The observations are discussed further in later sections when we can compare them with our own results.

The scaling relation of Equation 1 has been tested using hydrodynamical simulations with varying degrees of success. We will discuss these further in Section 6.1. On the basis of this, various authors have sought to constrain cosmological parameters from the observed X-ray cluster temperature and luminosity functions, for example, Lilje (1992), Viana & Liddle (1996), Eke, Cole & Frenk (1996), Henry (1997), Markevitch (1998), Eke et al. (1998), and Viana & Liddle (1999).

Because the scaling laws of galaxy clusters are used so extensively as the basis of cosmological models, we are undertaking a series of N-body, hydrodynamical simulations to investigate the X-ray properties of clusters of galaxies in detail. It is intended eventually to produce catalogues of clusters covering at least two decades in mass, for several different variants of the CDM cosmology, and with a variety of physical processes acting on the intracluster medium. In this current paper we report on the first of these catalogues for a non-radiative simulation in the $`\tau `$CDM cosmology. The catalogue is generated by combining three simulations (1283 particles each of gas and dark-matter) in which we resolve over 400 clusters at moderate resolution. Thus we have the correct boundary conditions for large-scale structure and a reasonable dynamic range with which to test the scaling relations.

We look at the circular velocity profiles of the clusters to see whether they all have the same form, such as the commonly-used NFW model. We find both that individual clusters show a wide range of profiles and that the parameters describing an average cluster are a function of mass. Thus, there is no a priori reason to suppose that the cluster population should follow the scaling law described in Equation 1. Nevertheless, it turns out, somewhat fortuitously, that the deviations of the mean cluster population from the scaling law are small, although the scatter from individual clusters can be quite large. When measured within a fixed physical radius, however, the dependence of temperature upon mass is steeper than that given by Equation 1.

Because this current simulation does not include radiative cooling, it cannot be an accurate reproduction of the real Universe. Nevertheless, we have deliberately kept the model simple in order to investigate how accurately, or otherwise, this simple model follows the self-similar scaling laws between temperature and mass111Because we are complete in mass, we quote $`T`$-$`M`$ rather than $`M`$-$`T`$ relations. Future papers will then look at the additional effects of cooling and non-gravitational heating.

We briefly describe the simulations and cluster identification in Section 2, and also discuss our choice of softening length. In Section 3, we look at the density profiles of the clusters and compare them with the NFW model. We derive predictions for the scaling relations in the NFW model in Section 4, then compare these with those of our simulated clusters in Section 5. We compare our results with previous work in Section 6 and summarise our conclusions in Section 7.

## 2 Numerical method

### 2.1 The simulations

We have carried out three simulations with 1283 particles each of gas and dark matter. The cosmological parameters were as follows: density parameter, $`\mathrm{\Omega }=1`$; cosmological constant, $`\mathrm{\Lambda }=0`$; power spectrum shape parameter, $`\mathrm{\Gamma }=0.21`$; and a linearly-extrapolated root-mean-square dispersion of the density fluctuations on a scale 8 $`h^1`$ Mpc, $`\sigma _8=0.60`$. The Hubble parameter is irrelevant as we did not allow the gas to cool; the baryon fraction was set to a low value, $`\mathrm{\Omega }_b=0.06`$ so that the gas makes only a minor contribution to the gravitational potential. The three simulations had different box-sizes, corresponding to different mass-resolutions, as listed in Table 1. We also carried out two further simulations of the middle-sized box to test the effect of changing the gravitational softening: these are shown in italics in the Table.

We use a parallel version of the Hydra N-body/SPH code (Couchman, Thomas & Pearce, 1995; Pearce & Couchman, 1997). The simulations were executed on the Cray T3E at the Edinburgh Parallel Computing Centre as part of the Virgo Consortium’s programme of investigations into large-scale structure.

### 2.2 Cluster identification

Initially we identify clusters in our simulation by searching for groups of dark matter particles within an isodensity contour of 180, as described in Thomas et al. (1998). We work with a preliminary catalogue of all objects with more than 30 particles, then retain only those which have a total mass within the virial radius exceeding $`M_{\mathrm{lim}}`$, corresponding to 500 particles of each species. The use of a small mass for the preliminary cluster selection ensures that our catalogue is complete. We have checked that using a different isodensity threshold, a different selection algorithm, or using gas particles instead of dark-matter particles to define the cluster, leads to an almost identical cluster catalogue—the only difference being the merger or otherwise of a small number of binary clusters.

We define the centre of the cluster to be the position of the densest gas particle. This will usually correspond to the peak of the X-ray emission and has the advantage that it is independent of the cluster selection method.

### 2.3 Substructure statistic

Much of the modelling that we will do on the clusters supposes that they are smooth and spherically-symmetric. In practice most clusters show some degree of substructure. Following Crone et al. (1996) and Thomas et al. (1998), we measure substructure by comparing the positions of the density maximum, $`𝒓_\mathrm{d}`$, and the centroid, $`𝒓_\mathrm{c}`$ of the cluster, where the latter is averaged over all particles within an isodensity contour of 180 times the background density. More specifically

$$S=\frac{|𝒓_\mathrm{d}𝒓_\mathrm{c}|}{r_{180}},$$
(2)
where $`r_{180}`$ is the radius of a sphere enclosing a mean density equal to 180 times the background density.

In some of the figures that follow, we indicate the clusters with the more prominent substructure, $`S>0.2`$, by plotting them with open symbols: 14 per cent of the clusters fall into this category.

### 2.4 Choice of softening

In any N-body simulation, it is necessary to introduce a gravitational softening in order that 2-body interactions do not become important. Thomas & Couchman (1992) estimated the minimum ratio of the 2-body relaxation time to the age of the Universe (which occurs in the core of the cluster, close to the softening radius) to be

$$\frac{t_{r,\mathrm{min}}}{t_0}1.3N_s^{\frac{1}{2}}\left(\frac{\mathrm{soft}}{10h^1\mathrm{kpc}}\frac{10h^1\mathrm{Mpc}}{\mathrm{box}}\right)^{\frac{3}{2}}$$
(3)
where we have put in parameters appropriate to the current simulations. Here $`N_s`$ is the number of particles within $`\sqrt{e}`$ times the softening length. The values given in Table 1 are for the smallest clusters extracted from each run; the relaxation time scales as $`T^{1/2}`$ for larger clusters because these have more particles within the softening length.

In a similar calculation to that of Thomas & Couchman, Steinmetz & White (1997) estimated the mean 2-body heating timescale for the cluster as a whole. They found that

$$t_r\sqrt{\frac{3}{\pi }}\frac{3N}{32\mathrm{ln}\mathrm{\Lambda }}t_c,$$
(4)
where $`N/2`$ is the number of particles within the half-mass radius, $`R_h`$, of the halo, $`\mathrm{ln}\mathrm{\Lambda }3`$ is the Coulomb logarithm, and $`t_c=2\pi R_h/v_c`$ is the orbital period at $`R_h`$, where $`v_c`$ is the circular velocity. The dependence upon the softening is much smaller than in Equation 3 which was measuring numerical relaxation in the core of the cluster. Putting in numbers appropriate to our halos, we find the heating timescales for the lowest-mass clusters in each run are approximately 10 times the age of the Universe (the heating rate scales roughly in inverse proportion to mass for higher-mass clusters). Thus we might expect the gas to be heated by up to 10 percent relative to the dark matter (the results of the tests, reported below, suggest that the heating is slightly smaller than this).

Because theoretical estimates of the numerical heating are so uncertain, we tested the sensitivity of our results to changes in the softening. Three simulations of the 112.9 $`h^1`$Mpc box were carried out which differed only in having softenings of 10, 50 and 100 $`h^1`$kpc, as shown in Table 1. For each cluster detected in the simulations, we measured the mean specific energy of particles within the virial radius (the radius of a sphere that encloses 180 times the mean density). Figure 1 shows $`\beta _{180}`$, the ratio of the mean specific energy of all particles to that of the gas particles.

As can be seen from the Figure, the effect of using too small a softening can be severe. Putting soft=$`10h^1`$kpc into Equation 3 gives $`t_{r,\mathrm{min}}/t_0`$0.12–0.24 for the lowest-to-highest mass clusters shown in the Figure. The Steinmetz & White criterion makes only a very slight distinction between the three runs with different softening (which enters Equation 4 only via the Coulomb logarithm) and so should be regarded as a necessary but not sufficient condition to prevent artificial heating.

Figure 2 shows the profile of $`\beta `$ as a function of radius (i.e. $`\beta (r)`$ is the ratio of the mean specific energy of all particles to the mean specific energy of gas particles within spherical shells about the cluster centre). Because individual profiles show large variations in $`\beta `$, we have averaged the profiles of the 70-80 smallest clusters, $`M<2\times 10^{14}h^1`$$`\mathrm{M}_{}`$, that show minimal substructure ($`S<0.2`$: see Section 2.3)—these are the clusters with the lowest values of $`t_{r,\mathrm{min}}/t_0`$. There is no evidence for core heating in the runs with softenings of 50 and 100 $`h^1`$kpc, whereas for the 10 $`h^1`$kpc softeing run $`\beta `$ is suppressed not only in the core but throughout the cluster. This suggests that much of the heating is going on in sub-halos before the cluster forms. As the properties of these sub-halos are similar to those of the final cluster core, this would explain why the Thomas & Couchman formula (Equation 3) is more appropriate than that of Steinmetz & White (Equation 4).

It is clear from Figures 1 and 2, that there is little difference between the results for softenings of 50 and 100 $`h^1`$kpc. They have similar values of $`\beta _{180}`$ and the radial variation in $`\beta `$ is minimal. We conclude that a softening of 50 $`h^1`$kpc is adequate for our purposes.

The softenings for each of the other two boxes were chosen to give a similar central value of $`t_{r,\mathrm{min}}/t_0`$ for the lowest-mass clusters as in the 50 $`h^1`$kpc softening run. The values of $`\beta _{180}`$ for each of the three production simulations are shown in Figure 3. Clusters with significant substructure, $`S>0.2`$, are shown as open symbols; smooth clusters with solid symbols. In this figure and all those that follow, clusters from the simulations of side 50.0, 112.9 and 153.0 $`h^1`$Mpc are shown as triangular, circular and square symbols, respectively.

Navarro & White (1993) and Pearce, Thomas & Couchman (1994) have shown that gas can gain energy from dark-matter during gravitational collapse even in the absence of numerical heating. The mechanism is quite straight-forward: both the gas and dark matter are stirred during the collapse, however the gas is able to thermalise its kinetic energy, then pick up more at the expense of the dark matter. There is no reason to suppose that this process should have occurred equally in low- and high-mass clusters. There are hints of a slight decline in $`\beta _{180}`$ at low masses in Figure 3, but to a good approximation, $`\beta _{180}`$ is independent of mass and equal to 0.94$`\pm 0.03`$. Thus the fractional heating is comparable to, but slightly smaller than, the estimate of 10 per cent that we made from the Steinmetz & White formula.

Because the gas has a greater specific energy than the dark matter, it is more extended and hence the baryon fraction within the virial radius is smaller than the global average. We measure it to be approximately $`(0.85\pm 0.02)\mathrm{\Omega }_b`$ in our simulated clusters, independent of mass.

## 3 Cluster profiles

In clusters, the dark matter is dynamically dominant. It is therefore important to have a good model of the dark-matter mass distribution. Navarro, Frenk & White (1995, 1996, 1997) showed that the density profiles of simulated clusters in a wide variety of cosmological models are well-described by the formula

$$\rho =\frac{\rho _0}{x(1+x)^2},$$
(5)
where $`r=ax`$ is the radius and $`\rho _0`$ and $`a`$ are free parameters, and went so far as to describe this as a ‘universal density profile’. This formula must break down at large radii as it predicts infinite mass, but it appears to hold out to the virial radius, defined (in critical-density models) as the radius, $`r_{180}`$, enclosing a mean overdensity of 180.

In a previous paper Thomas et al. (1998) found that the NFW formula was indeed a good approximation to the mass distribution of an average cluster, but that there was a wide dispersion in the rate at which the density was declining at the virial radius. Hence we introduce a more general profile

$$\rho =\frac{\rho _0}{x(1+x)^s},$$
(6)
where $`s`$ is a constant. Because the name is now so well-established in the literature, we call this the ‘generalised NFW’ model, although we note that the case $`s=3`$ was introduced first, by Hernquist (1990).

The resolution of the simulations presented in this paper is not sufficient to determine the density profile in the cores of the clusters. This has been looked at in depth by Moore et al. (1998) and by Jing & Suto (2000). There is now evidence to suggest that the inner density cusp rises more steeply than $`\rho r^1`$, but that is of little consequence for the overall dynamics of the cluster and so we will stick with Equation 6 here.

### 3.1 Best-fit circular speed profiles

The density profiles of individual clusters are often very noisy which makes them difficult to match to any given theoretical profile. It is better to use a cumulative profile such as the circular speed, $`v_\mathrm{c}`$, which is anyway a dynamically more relevant quantity. For the above distribution,

$$v_\mathrm{c}^2\frac{(1+x)^{s1}(s1)x1}{x(1+x)^{s1}}$$
(7)
when $`s2`$, and

$$v_\mathrm{c}^2\frac{(1+x)\mathrm{ln}(1+x)x}{x(1+x)}$$
(8)
when $`s=2`$.

An example of a circular-speed profile is given in Figure 4 together with the best-fit function of the form of Equation 7. We have fitted and plotted the curve between twice the softening and twice the virial radius. This particular cluster was chosen because the goodness-of-fit as measured by the mean-square deviation from the theoretical curve is the median value for all the clusters. As can be seen, there are kinks in the profile, showing evidence of substructure: this is typical of most of the clusters in our catalogue.

The best-fit values of $`s`$ as a function of cluster mass are shown in Figure 5. Note that $`s`$ has been limited to be less than or equal to 10.

The asymptotic density profile at large radii has a slope of $`(s+1)`$. However, this is not always representative of the slope at two virial radii, the outer radius to which we fit the rotation curve. Hence we plot with solid symbols in the Figure only those clusters for which the two differ by less than $`0.2s`$—this corresponds to a characteristic radius $`a<0.5r_{180}`$.

A quick glance at Figure 5 makes the suggestion that $`s=2`$ is a universal density profile seem surprising. However, it is hard to measure the density profiles in the outer parts of clusters with any degree of accuracy and the answer that one gets often depends upon the radial extent of the fit. At one virial radius, the slope of the density profile is far from $`(s+1)`$ and so the asymptotic slope is poorly constained. In addition, $`s`$ and $`a`$ are strongly correlated and it is often possible to get a reasonable fit by forcing $`s=2`$ and allowing $`a`$ to vary. Hence the statement that the profile within one virial radius is ‘consistent with and NFW profile’ is largely meaningless. It is for this reason that we choose to fit the profile within two virial radii instead.

We define clusters to be consistent with an NFW profile if the asymptotic slope of their density profile lies between $`2.8`$ and $`3.2`$ (i.e. $`1.8<s<2.2`$). Just under a quarter of the clusters meet this criterion.

The best-fit profiles of many clusters plotted with open symbols show a high value of $`s`$. However, this does not indicate steep density profiles at large radii because the best-fit core radii rise to compensate. Rather, it indicates that the functional form of the generalised NFW profile is a poor representation of the cluster. As an example consider the cluster shown in Figure 6.

This is a smooth cluster: visually it appears spherically-symmetric and its substructure parameter is small, $`S=0.038`$. In addition, the best-fitting ellipsoid (see Thomas et al. 1998) is amongst the most spherical of any cluster with axial ratios of 1.14:1.0:0.94. However, the density profile, even out to 0.8 virial radii, is poorly fit by an NFW profile (dashed line). Also, it shows a sharp change in slope at this radius that cannot be matched by any generalised NFW profile. The dotted line shows a generalised NFW model with $`s=10.`$ and $`a=3.1r_{180}`$. In order to reproduce the sharp decline in density at the virial radius, $`s`$ has to be very large, but this then leads to a density profile that is declining far too rapidly in the outer parts of the cluster (and would steepen even more at radii larger than those shown in the Figure). A better representation of the density profile in this case is given by the dotted line which corresponds to the function

$$\rho =\frac{\rho _0}{x(1+x^2)^{s/2}},$$
(9)
where $`x=r/a`$ as before, and $`s=3.2`$, $`a=0.40r_{180}`$. This value of $`s`$ is a much better estimate of the asymptotic slope of the density profile at large radii.

It can be seen from Figure 5 that the line $`s=3`$ roughly separates the solid from the open symbols in the upper half of the plot. Thus, where the outer slope of the density profile is well-defined, it generally lies between $`3`$ (an NFW profile) and $`4`$ (a Hernquist profile). The open symbols represent clusters, like that shown in Figure 6, that have a sharper break in their density profile than can be fit by a generalised NFW profile: these comprise about 30 per cent of the total cluster sample.

The open symbols that correspond to values of $`s`$ less than 2 are mostly clusters that show some degree of substructure, for which the density profile is not well-defined. These comprise another 16 percent of the cluster population.

Although the spread in $`s`$ is large, there is a weak trend for $`s`$ to increase with mass. To make this more evident, we define an average low-mass and an average high-mass cluster by selecting all relatively smooth clusters, $`S<0.2`$, in the mass ranges $`M<1.2\times 10^{13}h^1`$$`\mathrm{M}_{}`$and $`M>3.0\times 10^{14}h^1`$$`\mathrm{M}_{}`$. The resulting profiles are extremely well-fit by our theoretical model with slopes of $`s=2.0`$ and $`s=2.3`$, respectively.

It is clear from the above analysis that there is no ‘universal profile’ for dark matter halos. A substantial proportion of clusters show obvious substructure, and even those that don’t exhibit a wide variety of functional forms for the spherically-averaged density profiles of halos.

### 3.2 Best-fit NFW profiles and concentrations

The concept of a universal density profile is an attractive one. It makes modelling of observed clusters much simpler and it has the advantage that there is only one free parameter—the ratio of the characteristic radius in the NFW formula to the virial radius, $`x_{180}=r_{180}/a(s=2)`$.222NFW define $`x_{200}`$ to be the ‘concentration parameter’, presumably using 200 as an approximation for the virial overdensity; in this paper we will use the term to stand for $`x_{180}`$ instead—there is little difference between the two. Therefore, we wish to see how well one can approximate cluster properties by assuming that they all follow the NFW profile, in defiance of the results of Section 3.1.

The best-fit values of $`x_{180}`$ as a function of mass are shown in Figure 7. There is a general trend of decreasing concentration as one moves to higher masses, but once again the scatter is large. The solid line in Figure 7 corresponds to the function

$$x_{180}=3.7\left(\frac{M}{10^{15}h^1\mathrm{M}_{}}\right)^{0.1}.$$
(10)
The 86 per cent of the clusters with $`S<0.2`$ are spread equally above and below the line. We shall use this relation in the analysis that follows to see how well the simple NFW model predicts the measured scaling relations between temperature and mass.

We can compare our concentration parameters to those of Navarro, Frenk & White (1997), by measuring halo mass in terms of $`M_{}`$, the mass of a spherical region with root-mean-square linear density fluctuation of 1.69 today. For our simulations, $`M_{}3.3\times 10^{12}h^1\mathrm{M}_{}`$, and $`M_{180}/M_{}`$ ranges from 3 to 300. This gives concentrations in good agreement with their CDM models.

Interestingly, Figure 8 shows that the measured concentrations are highly-correlated with the substructure statistic. Thus low concentrations do not normally occur in spherically-symmetric clusters, but come from spherical averaging of clusters that have substructure. This suggests that the weak trend of decreasing concentration with increasing mass is driven primarily by the fact that substructure increases with mass, as shown in Figure 9.

We conclude this section by plotting in Figure 10 the temperature profiles for the mean clusters described above. The lower curve in each panel shows the gas temperature in units of keV$`/k`$, while the upper curve shows the total specific energy of the gas (thermal plus kinetic) in the same units. It should be noted that the gas is approximately isothermal within about one-third of the virial radius, but that its temperature drops rapidly in the outer parts of the cluster. This is a reflection of the fact that the specific energy of the dark matter declines in the outer regions of the cluster, because the ratio between the two is approximately constant out to two virial radii.

It is a matter of some debate in the literature as to whether the observed temperature profiles of clusters are isothermal or decline to large radii. Using ASCA data, Irwin & Bregman (2000) and White (2000) infer flat or slowly increasing temperature profiles within 0.3 $`r_{180}`$, whereas Markevitch et al. (1998) find that the temperature drops by a factor of two within 0.5 $`r_{180}`$. It is clear from our simulations that isothermal profiles beyond about 0.3 $`r_{180}`$ are inconsistent with our simple model and would probably require some form of non-gravitational heating.

## 4 Predicted scaling relations

In this section, we predict scaling relations for model clusters under the assumptions that they are spherically-symmetric, isolated and in hydrodynamical equilibrium. It is not that we are asserting that these assumptions are true, but rather we want to see how well such a simple model will perform and to test the sensitivity of the scaling relations to changes in the model parameters.

By scaling everything in terms of a characteristic radius, $`a`$, and density $`\rho _0`$, for each cluster, we can make all the variables dimensionless. We write $`r=ax`$ and $`\rho =\rho _0\stackrel{~}{\rho }`$, where

$$\stackrel{~}{\rho }=\frac{1}{x(1+x)^s}.$$
(11)
and the tilde indicates a dimensionless quantity. Similarly, for the mass within radius $`r`$,

$$m=4\pi \rho _0a^3\stackrel{~}{m}$$
(12)
where

$$\stackrel{~}{m}=_0^{r/a}\stackrel{~}{\rho }x^2\text{d}x.$$
(13)
When we observe clusters, we usually measure their properties at or within the virial radius, $`r_{\mathrm{vir}}=ax_{\mathrm{\Delta }_c}`$ at which their enclosed density equals some multiple, $`\mathrm{\Delta }_c`$, of the critical density of the universe. For the $`\tau `$CDM model that we consider in this paper, $`\mathrm{\Delta }_c=180`$, but we will set $`r_\delta =ax_\delta `$ and leave $`\delta `$ as a free parameter, so as to be completely general in our argument.

From the above definitions, we can write $`\rho _0`$ in terms of $`x_\delta `$,

$$\rho _0=60\delta _{180}\rho _c\frac{x_\delta ^3}{\stackrel{~}{m}_\delta },$$
(14)
where $`\rho _c`$ is the critical density, $`\delta _{180}=\delta /180`$ and $`\stackrel{~}{m}_\delta =\stackrel{~}{m}(x_\delta )`$.

### 4.1 The virial temperature-mass relation

We next calculate the velocity dispersion, $`\sigma `$ as a function of radius for the NFW profile. To do this, we assume that the velocity ellipsoid is isotropic. It then becomes convenient to write $`\sigma ^2=kT/\mu m_H`$, where $`k`$ is the Boltzmann constant, $`\mu m_H`$ is the mean particle mass for the gas, and $`T`$ is the ‘dynamical temperature’ of the halo. If the specific energy of the gas and dark-matter were the same, then $`T`$ would be equal to the gas temperature.

The spherically-symmetric Jeans’ Equation for a fluid with an isotropic velocity dispersion is

$$\frac{\text{d}(\stackrel{~}{\rho }\stackrel{~}{T})}{\text{d}x}=\frac{\stackrel{~}{m}\stackrel{~}{\rho }}{x^2},$$
(15)
where

$$T=\frac{G\mu m_H}{k}\mathrm{\hspace{0.17em}4}\pi \rho _0a^2\stackrel{~}{T}.$$
(16)
We are now in a position to calculate the mean dynamical temperature, $`T_\delta `$, within radius $`r_\delta `$. By definition,

$$\stackrel{~}{m}_\delta \stackrel{~}{T}_\delta =_0^{\stackrel{~}{m}_\delta }\stackrel{~}{T}\text{d}\stackrel{~}{m}=_0^{x_\delta }\stackrel{~}{\rho }\stackrel{~}{T}x^2\text{d}x.$$
(17)
Combining Equations 12 and 16 to eliminate $`a`$, and then substituting for $`\rho _0`$ from Equation 14, gives

$`kT_\delta `$ $`=`$ $`G\mu m_H(240\pi \rho _c\delta _{180})^{\frac{1}{3}}{\displaystyle \frac{x_\delta \stackrel{~}{T}_\delta }{\stackrel{~}{m}_\delta }}m_\delta ^{\frac{2}{3}}`$ (18)
$``$ $`16.0\delta _{180}^{\frac{1}{3}}{\displaystyle \frac{x_\delta \stackrel{~}{T}_\delta }{\stackrel{~}{m}_\delta }}\left({\displaystyle \frac{m_\delta }{10^{15}h^1\mathrm{M}_{}}}\right)^{\frac{2}{3}}\text{keV},`$
where we have taken $`\mu m_H10^{27}`$kg, corresponding to a fully-ionized cosmic mix of elements.

When $`\delta =\mathrm{\Delta }_c`$, then $`T_\delta `$ is the ‘virial temperature’. We will use the notation $`T_{180}`$ rather than $`T_{\mathrm{vir}}`$, to make explicit that we are referring to an overdensity of 180 (as in other cosmologies, the virial radius would correspond to a different overdensity).

If all clusters are self-similar (e.g. if they all follow the NFW model), and if we measure the total mass, $`m_\delta `$, and mean dynamical temperature, $`T_\delta `$, within $`r_\delta `$, then it is often assumed that they must obey the relation $`T_\delta m_\delta ^{2/3}`$. This is not strictly true, however, because of the presence of the $`x_\delta \stackrel{~}{T}_\delta /\stackrel{~}{m}_\delta `$ term which varies with concentration, $`x_\delta `$.

We show the variation of $`x_\delta \stackrel{~}{T}_\delta /\stackrel{~}{m}_\delta `$ with $`x_\delta `$ in Figure 11 for values of $`s=1.5`$, 2 and 3. We see that for values of $`x_{180}`$ in the range 1.5 to 10 that are typical when fitting the NFW profile to clusters, $`x_{180}\stackrel{~}{T}_{180}/\stackrel{~}{m}_{180}`$ lies in the narrow range 0.4–0.47. Fitting $`s=1.5`$ profiles instead would give smaller values of $`x_{180}`$, and fitting $`s=3`$ profiles larger values, so that these too give approximately the same answer. Thus the virial temperature-mass relation is robust to variations in the density profile.

### 4.2 The gas density profile and the X-ray temperature-mass relation

The virial temperature-mass relation predicted above is of little use observationally. Instead we need to use the emission-weighted, X-ray temperature. Previous studies such as Eke, Navarro & Frenk (1998) and Makino, Susaki & Suto (1998) have calculated the mean X-ray temperature and the luminosity of the intracluster medium assuming that it is isothermal and sitting in an NFW potential. Figure 10 shows that the gas is far from isothermal and so we instead make the assumption, as indicated by Figure 2, that the ratio of the specific energies of the gas and dark matter is everywhere constant and equal to $`\beta =0.94`$, i.e. $`T_{\mathrm{gas}}=T/\beta `$.

We will further assume that the dark matter is dynamically dominant. Then

$$\frac{1}{\stackrel{~}{\rho }_{\mathrm{gas}}}\frac{\text{d}(\stackrel{~}{\rho }_{\mathrm{gas}}\stackrel{~}{T}_{\mathrm{gas}})}{\text{d}x}=\frac{1}{\stackrel{~}{\rho }}\frac{\text{d}(\stackrel{~}{\rho }\stackrel{~}{T})}{\text{d}x},$$
(19)
which can be expanded to give

$$\frac{\text{d}\mathrm{ln}\stackrel{~}{\rho }_{\mathrm{gas}}}{\text{d}\mathrm{ln}x}=\frac{\text{d}\mathrm{ln}\stackrel{~}{\rho }}{\text{d}\mathrm{ln}x}(1\beta )\frac{\text{d}\mathrm{ln}\stackrel{~}{P}}{\text{d}\mathrm{ln}x},$$
(20)
where $`\stackrel{~}{P}=\stackrel{~}{\rho }\stackrel{~}{T}`$. This shows that the gas is able to support itself more effectively and hence has a shallower density gradient than the dark matter. For $`\beta =0.94`$, the asymptotic density gradient at large radii is $`2.7`$, compared to $`3`$ for the dark matter.

The mean X-ray temperature within radius $`r_\delta `$ is defined by

$$\stackrel{~}{L}_\delta \stackrel{~}{T}_{X,\delta }=_0^{\stackrel{~}{L}_\delta }\stackrel{~}{T}_{\mathrm{gas}}\text{d}\stackrel{~}{L},$$
(21)
where

$$\stackrel{~}{L}(x)=_0^x\stackrel{~}{\rho }_{\mathrm{gas}}^2\stackrel{~}{T}_{\mathrm{gas}}^{\frac{1}{2}}\text{d}x$$
(22)
and $`\stackrel{~}{L}_\delta =\stackrel{~}{L}(x_\delta )`$. Note that we have assumed that the emissivity of the gas scales as $`\stackrel{~}{T}_{\mathrm{gas}}^{1/2}`$, although this is only really true at high temperatures where bremsstrahlung dominates. However, the dominant contribution to the variation in the integrand in Equation 22 comes from $`\stackrel{~}{\rho }_{\mathrm{gas}}`$ and so this is a reasonable approximation.

The values of the normalisation factor in the X-ray temperature-mass relation, $`x_\delta \stackrel{~}{T}_{X,\delta }/\stackrel{~}{m}_\delta `$, are shown in Figure 12. They can be seen to rise more steeply at large concentrations than in Figure 11. This is because the X-ray temperature declines far less steeply with radius than the dynamical one used previously. Figure 12 suggests that the range of concentrations seen in our simulated clusters will add quite a lot of scatter to the mean X-ray temperatures within the virial radius, with the more concentrated clusters being hotter than the less-concentrated ones for a given mass by up to 40 per cent.

## 5 Measured scaling relations

### 5.1 The virial temperature-mass relation

The location of each of our clusters on the $`T_{180}`$-$`m_{180}`$ plane is shown in Figure 13. Also shown by the dashed line is the best-fit power-law

$$kT_{180}6.9\left(\frac{m_{180}}{10^{15}h^1\mathrm{M}_{}}\right)^{0.67\pm 0.02}\text{keV},$$
(23)
which corresponds to a value of $`x_\delta \stackrel{~}{T}_\delta /\stackrel{~}{m}_\delta =0.43`$.

The predicted relation from Equation 18, using the median concentrations, $`x_{180}`$, taken from Equation 10, is shown by the solid line in the Figure. It can be seen that this predicts temperatures that are slightly too low, although the spread in temperatures is consistent with the spread in $`x_{180}`$ seen in Figure 7.

It would be possible to reconcile the predicted and measured relations by using slightly larger values of $`x_{180}`$: Figure 8 shows that the most regular clusters are biased towards higher concentrations. An alternative explanation for the low virial temperatures is that the dark matter velocity ellipsoid is not isotropic, as radial orbits are less efficient at supporting the particles within a given potential than transverse ones. This is discussed further in Section 6.2.

### 5.2 X-ray versus virial temperatures

The ratios of the dynamical temperatures to the X-ray temperatures, averaged within the virial radii, are shown in Figure 14.

The prediction from the isotropic NFW model with concentrations given by Equation 10 is shown by the solid line. It can be seen that the theoretical curve mimics the data in that it predicts that low-mass clusters should have lower values of $`\beta _{X,180}`$ than high-mass ones. However the mean measured values of $`\beta _{X,180}`$ lie about 0.1 below the predictions. Once again, this could be reconciled by using higher values of $`x_{180}`$ or an anisotropic dark matter velocity ellipsoid.

Although virial temperatures are notoriously difficult to measure, there is empirical evidence that the observed, X-ray temperatures of the intracluster medium are greater than the virial temperatures (e.g. Edge & Stewart 1991a,b; Bahcall & Lubin 1994). The latter give a value of $`\beta _{X,180}=0.94`$ for the highest-mass clusters, which is not too dissimilar from the values that we find in Figure 14. However, the results of the next section suggest values that are lower than this.

### 5.3 X-ray temperature-mass relations within a fixed overdensity

The relationship between X-ray temperature and mass within the virial radius is shown in Figure 15. The dashed line shows the best-fit relation from Evrard, Metzler & Navarro (1996, hereafter EMN),

$$kT_{X,180}7.85\left(\frac{m_{180}}{10^{15}h^1\mathrm{M}_{}}\right)^{0.67\pm 0.02}\text{keV}.$$
(24)
They extracted 58 clusters in the temperature range 1–10 keV, from three different sets of cosmological simulations with a variety of cosmological parameters. We have used the information given in their paper to interpolate their results to an overdensity of 180—our results are in good agreement.

The temperature normalisation is slightly higher than for the dynamical temperature-mass relation and corresponds to $`x_\delta \stackrel{~}{T}_\delta /\stackrel{~}{m}_\delta =0.49`$. The prediction from Equation 18 using values of $`x_{180}`$ taken from Equation 10 is shown by the solid line on the Figure. Once again it gives temperatures that are slightly too low.

The most extensive observational investigation of the $`T_X`$-$`M_{180}`$ relation is by Horner, Mushotzky & Scharf (1999) who used many different ways to determine the mass within the virial radius. Their preferred measure, based on X-ray emissivity and temperature profiles, gives $`T_XM_{180}^{2/3}`$, while those obtained from the isothermal $`\beta `$-model and X-ray surface brightness deprojections are flatter, $`T_XM_{180}^{0.50.56}`$. They suggest that the slopes of these latter two measures are too low because they assume a dark matter density profile of $`\rho _{\mathrm{dark}}r^2`$ whilst observations and simulations suggest that it is steeper, $`\rho _{\mathrm{dark}}r^{2.4}`$.

The preferred relation from Horner, Mushotzky & Scharf (1999) is shown as a dotted line in Figure 15. From this it is clear that either the measured cluster masses are too low, or more likely dissipationless simulations predict X-ray temperatures that are smaller than the observed values. If we combine the results from Figures 14 and 15 then this argues that observed high-mass clusters have values of $`\beta _{X,180}2/3`$.

This high value of $`\beta _{X,180}`$ might be thought to provide evidence for heating of the intracluster medium, but in fact radiative cooling can have the same effect! Pearce et al. (2000) have shown that the removal of low-entropy gas by cooling in the core of the cluster can raise the X-ray temperature of the cluster by 20 to 40 percent. This cooling is relatively more important in low-mass clusters and so would have the effect of flattening the $`T_X`$-$`M`$ relation slightly.

Mohr, Mathiesen & Evrard (1999) find in a sample of 45 clusters that the ICM is more extended than the dark-matter and that $`T_XM_{\mathrm{gas}}^{1/2}`$. In our simulations, we also find that the gas is more extended than the dark matter, due to its higher specific energy. However, the baryon fraction within the virial radius is approximately 0.85, independent of mass, and so obtain a temperature-gas mass relation that parallels the one for the total mass, $`T_XM_{\mathrm{gas}}^{0.67}`$. Physical processes such as heating or radiative cooling are once again required (and act in the correct sense) to reconcile the observations and simulations.

In Figure 16, we show a similar plot to Figure 15, but for the X-ray temperature-mass relation within an overdensity contour of 1000.

Once again, the isotropic NFW model slightly underpredicts the X-ray temperatures (but agrees with the results of EMN). The best-fit power law is consistent with the self-similar prediction

$$kT_{X,1000}13.0\left(\frac{m_{1000}}{10^{15}h^1\mathrm{M}_{}}\right)^{0.67}\text{keV}.$$
(25)
The dotted line shows the observational results from Nevalainen, Markevitch & Forman (2000). The observed temperatures are again higher than the predictions from this non-radiative model but in a way that is now mass-dependent—this is consistent with our expectation that the effects of cooling and/or heating would be greater in lower-mass clusters.

### 5.4 The X-ray temperature-virial radius relation

From an observational point of view, it is much easier to measure the virial radius rather than the virial mass, and so it is surprising that more authors do not report the former. There is no extra information here, however, as the two are simply related and so it is possible to rewrite Equations 24 and 25 in terms of virial radius and X-ray temperature:

$$r_{180}1.9\left(\frac{kT_{X,180}}{10\text{keV}}\right)^{\frac{1}{2}}h^1\text{Mpc},$$
(26)
$$r_{1000}0.85\left(\frac{kT_{X,1000}}{10\text{keV}}\right)^{\frac{1}{2}}h^1\text{Mpc}.$$
(27)
These follow with the self-similar form in that they have slopes of $`\frac{1}{2}`$. $`r_{180}`$ agrees with the results of EMN; $`r_{1000}`$ has a slightly lower normalisation.

We note that the only observational measurement of $`r_{1000}`$ by Vikhlinin, Forman & Jones (1999) is at first sight in gross disagreement with our prediction. They have

$$r_{1000}0.415\left(\frac{kT_{X,1000}}{10\text{keV}}\right)^{0.57}h^1\text{Mpc}$$
(28)
which gives a value at 5 keV of 1.04 $`h^1`$Mpc—far in excess of the simulation results. The reason for this discrepancy is that Vikhlinin et al. have defined overdensity with respect to a mean baryon abundance of $`\mathrm{\Omega }_bh^2=0.010`$ which is much lower than either recent determination from primordial nucleosynthesis ($`\mathrm{\Omega }_bh^2=0.019`$, Tytler et al. 2000) or from X-ray measurements of the baryon fraction ($`\mathrm{\Omega }_bh^2=0.0594h^{1/2}`$, Ettori & Fabian 1999). We can turn the argument around and ask what baryon abundance would make our results compatible with Vikhlinin et al.. The radii they quote enclose a baryon overdensity of about 230 in our models. Hence we require

$$\mathrm{\Omega }_bh^20.010\frac{1000}{230}0.043,$$
(29)
This agrees with the value from Etorri & Fabian for $`h=0.54`$.

### 5.5 Temperature-mass relations within the Abell radius

If we measure cluster properties within a fixed radius, rather than a fixed overdensity, then the deviation from the self-similar scaling relations can be quite large. This is shown in Figure 17 where we plot, for within the Abell radius, emission-weighted temperature versus mass. The best-fit power law, shown as the dashed line, has a slope of 0.81.

The reason for the steeper slope is that, for low-mass clusters, the Abell radius is greater than the virial radius and so we are averaging properties over a larger volume than before. This has the effect of lowering the X-ray temperature slightly (because the X-ray temperature is heavily weighted by emission from the centre of the cluster this effect is small), but greatly increasing the mass. For high-mass clusters, however, the virial radii are similar to the Abell radii and so there is no change.

Because the X-ray temperatures are almost unaltered, the mass-temperature relation can be estimated by simply multiplying the right-hand-side of Equation 24 by the ratio of the mass within an Abell radius to that within the virial radius, $`M_{\mathrm{Abell}}/M_{180}`$. We have done this, using Equation 26 to estimate the virial radii and assuming that the mass profile follows the NFW model. The result, shown as the solid line in Figure 17, closely follows the data points.

Figure 18 is the same as for Figure 17, but for properties averaged within half an Abell radius. Once again the best-fit power law has a slope, 0.86, that is steeper than the self-similar relation, and the corrected NFW model provides an excellent fit to the data.

## 6 Discussion

### 6.1 Comparison with previous simulations

Previous authors have also looked at scaling relations in simulations of clusters without radiative cooling. All agree that the gas is more extended than the dark matter, but disagree on the details of the $`T_X`$-$`M`$ relation. The studies split into two types: those that simulate a small number of clusters in detail, and those that, like us, look at a large number of clusters at moderate resolution.

Navarro, Frenk & White (1995) simulate only 6 clusters, extracted from a low-resolution simulation, but choose these to span a factor of thirty in mass and adjust the particle mass so that each cluster contains several thousand particles within the virial radius. They find a $`T_X`$-$`M`$ relation in agreement with the virial relation and values of $`\beta `$ of about unity (where their $`\beta `$, as ours, includes a contribution from unthermalised motions in the gas) with no dependence upon mass.

EMN have taken the Navarro, Frenk & White results and compared them with two other sets of simulations of isolated clusters from Metzler (1995) and Mohr et al. (1995). These span a range of cosmologies and use two distinct codes, yet yield very similar results. It should be noted that EMN do not quote results for an overdensity of 180, but they give sufficient information in their paper to allow one to interpolate from their preferred overdensity of 500. Once we do so, we find a $`T_X`$-$`M`$ relation essentially identical to ours.

Eke, Navarro & Frenk (1998) perform a similar study for 10 clusters in a $`\mathrm{\Lambda }`$-CDM cosmology, ranging over a factor of only three in mass. Their results are consistent with the scaling relations but cover too small a dynamic range to provide a strong constraint.

By contrast, Bryan & Norman (1998) use a smaller box and simulate a large number of clusters at relatively poor resolution. They find that $`\beta `$ is an increasing function of cluster mass which suggests that numerical heating is having an effect. They also find a $`T_X`$-$`M`$ relation that is flatter than the self-similar one, but only for an open cosmology.

Yoshikawa, Jing & Suto (2000) have undertaken similar simulations to ourselves, but with a poorer spatial resolution and covering a smaller mass-range. They too find the expected slope of $`\frac{2}{3}`$ for the $`T_X`$-$`M`$ relation.

### 6.2 Beyond the isotropic NFW model

It can be seen from Figures 13 and 15 that the isotropic NFW model underpredicts the virial temperature slightly and the X-ray temperature quite a lot.

One possible explanation for the low virial temperatures is that the dark matter velocity ellipsoid is not isotropic, as radial orbits are less efficient at supporting the particles within a given potential than transverse ones.

We have checked that a ratio of the transverse to radial velocity dispersion of $`\sigma _t=0.9\sigma _r`$ is sufficient to bring the dashed and solid lines in Figures 13 and 15 into close agreement. We choose not to report on this in detail, however, because we believe that it overly-complicates the simple model that we are trying to test. It should be remembered that few of our clusters are free from some form of substructure and this too would invalidate the simple NFW model.

We prefer to conclude that the isotropic NFW model does pretty well, whilst bearing in mind that it underestimates the temperatures slightly.

## 7 Conclusions

We extract clusters, covering more than 2 decades in mass, from three simulations of the $`\tau `$CDM cosmology. This represents the largest, uniform catalogue of simulated clusters ever produced.

We fit the circular-speed profiles of the clusters with models corresponding to spherically-symmetric density profiles of the form

$$\rho =\frac{\rho _0}{x(1+x)^s},$$
(30)
where $`s=2`$ corresponds to the renowned NFW profile.

We find that the scatter in the best-fit values of $`s`$ is large, with under a quarter of the clusters being accurately fit by the NFW model. The others either have steeper outer slopes in their density profiles, show a sharp break in their density profiles that cannot be fit by the above form, or have significant substructure. However, the mean, smooth, low-mass cluster does have $`s2`$, whereas the equivalent high-mass cluster has a steeper profile with $`s2.8`$.

When we force $`s=2`$, then the cluster concentrations show a large scatter, but the median concentration declines slightly with mass. This is driven primarily by an anti-correlation between concentration and substructure, with substructure being more prevalent in high-mass clusters.

We investigate how well this median NFW model predicts the cluster temperature-mass scaling relations, by deriving theoretical relationships between both dynamical and emission-weighted, X-ray temperature and mass in the isotropic NFW model, as a function of cluster concentration. The virial temperature-mass relation, averaged within a spherical region enclosing an overdensity of 180, closely mimics the self-similar form and lies a few per cent above the NFW prediction.

The gas in our clusters is hotter than the dark matter but we do not attach much significance to this: there will be a small amount of true and numerical heating due to heating by clumps and individual particles of dark matter. More importantly, we have chosen in this paper to neglect the effects of radiative cooling and heating associated with metal enrichment of the intracluster medium.

When measured within spheres enclosing a fixed overdensity, the X-ray temperature versus virial mass relation has a slope of approximately 0.67, in agreement with the self-similar prediction and with previous work (but extending over a greater mass-range). The normalisation is lower than that of the observations. This is almost certainly because we have chosen to neglect cooling and/or heating.

The radius-temperature relation is in agreement with the measured value of Vikhlinin, Forman & Jones (1999) provided that the baryon fraction is large, as indicated by Ettori & Fabian (1999).

When averaged within an Abell radius (or half an Abell radius), the X-ray temperature versus mass relation is steeper, with a slope of 0.81 (or 0.86). This is because the Abell radius is a greater multiple of the virial radius in low-mass clusters as compared to high-mass ones and so increases their measured mass.

We have not commented on the X-ray luminosity in this paper, because our simulations are not able to fully resolve the X-ray emission in the cluster cores. This core gas has a cooling time that is anyway less than the age of the Universe and so cannot be correctly modelled using a non-radiative simulation. We note that the results of Pearce et al.(2000) suggest that radiative cooling will act so as to *raise* the temperature of the intracluster medium and so tend to bring our simulated clusters into agreement with the observations.

In future papers we will contrast the dynamics of the gas and dark matter in clusters, consider the effect of heating and cooling processes, and compare the results from different cosmological models.

## Acknowledgments

The simulations described in this paper were carried out on the Cray-T3E at the Edinburgh Parallel Computing Centre as part of the Virgo Consortium investigations of cosmological structure formation. Interaction between authors was aided by a NATO Collaborative Research Grant, CRG 970081. OM is supported by a DPST Scholarship from the Thai government; PAT is a PPARC Lecturer Fellow; LO is a Daphne-Jackson Fellow, funded by the Royal Society.