text:
# Thermal conduction in cosmological SPH simulations

## 1 Introduction

Clusters of galaxies provide a unique laboratory to study structure formation and the material content of the Universe, because they are not only the largest virialized systems but are also believed to contain a fair mixture of cosmic matter. Among the many interesting aspects of cluster physics, their X-ray emission takes a particularly prominent role. It provides direct information on the thermodynamic state of the diffuse intracluster gas, which makes up for most of the baryons in clusters.

While the bulk properties of this gas are well understood in terms of the canonical $`\mathrm{\Lambda }`$CDM model for structure formation, there are a number of discrepancies between observations and the results of present hydrodynamical simulations. For example, a long standing problem is to understand in detail the scaling relations of observed clusters, which deviate significantly from simple self-similar predictions. In particular, poor clusters of galaxies seem to contain gas of higher entropy in their centres than expected (Ponman et al., 1999; Lloyd-Davies et al., 2000). This has been interpreted either to be evidence for an entropy injection due to non-gravitational processes (Loewenstein, 2000; Wu et al., 1999; Metzler & Evrard, 1994), or as a sign of the selective removal of low-entropy gas by gas cooling (Voit et al., 2002; Wu & Xue, 2002). Both processes combined could influence the thermodynamic properties of the ICM in a complex interplay (Tornatore et al., 2003; Borgani et al., 2003).

Another interesting problem occurs for the radial temperature profiles of clusters. Most observed clusters show a nearly isothermal temperature profile, often with a smooth decline in their central parts (Allen et al., 2001; Johnstone et al., 2002; Ettori et al., 2002). Nearly isothermal profiles are also obtained in adiabatic simulations of cluster formation (e.g. Frenk et al., 1999). However, clusters in simulations that include dissipation typically show temperature profiles that increase towards the centre (e.g. Lewis et al., 2000), quite different from what is observed.

Perhaps the biggest puzzle is that spectroscopic X-ray observations of the centres of clusters of galaxies have revealed little evidence for cooling of substantial amounts of gas out of the intracluster medium (e.g. David et al., 2001), even though this would be expected based on their bolometric X-ray luminosity alone (Fabian, 1994). The apparent absence of strong cooling flows in clusters hence indicates the presence of some heating source for the central intracluster medium. Among the proposed sources are AGN, buoyant radio bubbles (Churazov et al., 2001; Enßlin & Heinz, 2001), feedback processes from star formation (Bower et al., 2001; Menci & Cavaliere, 2000), or acoustic waves (Fujita et al., 2003).

Recently, Narayan & Medvedev (2001) have proposed that thermal conduction may play an important role for the cooling processes in clusters. The highly ionised hot plasma making up the ICM in rich clusters of galaxies should be efficient in transporting thermal energy, unless heat diffusion is inhibited by magnetic fields. If conduction is efficient, then cooling losses in the central part could be offset by a conductive heat flow from hotter outer parts of clusters, which forms the basis of the conduction idea.

Indeed, using simple hydrostatic cluster models where cooling and conductive heating are assumed to be locally in equilibrium, Zakamska & Narayan (2003, ZN henceforth) have shown that the central temperature profiles of a number of clusters can be well reproduced in models with conduction (see also Fabian et al., 2002; Voigt et al., 2002; Brüggen, 2003; Voigt & Fabian, 2003). The required conductivities are typically sub-Spitzer (Medvedev et al., 2003). This suggests that thermal conduction may play an important role for the thermodynamic properties of the ICM.

On the other hand, it has been frequently argued (Chandran & Cowley, 1998; Malyshkin & Kulsrud, 2001) that magnetic fields in clusters most likely suppress the effective conductivity to values well below the Spitzer value for an unmagnetised gas. Note that rotation measurements show that magnetic fields do exist in clusters (e.g. Vogt & Enßlin, 2003). However, little is known about the small-scale field configuration, so that there is room for models with chaotically tangled magnetic fields (Narayan & Medvedev, 2001), which may leave a substantial fraction of the Spitzer conductivity intact. The survival of sharp temperature gradients along cold fronts, as observed by Chandra in several clusters (Markevitch et al., 2000; Vikhlinin et al., 2001; Ettori & Fabian, 2000), may require an ordered magnetic field to suppress conduction. For a more complete review of the effects of magnetic fields on galactic clusters, see Carilli & Taylor (2002) and references therein.

It is clearly of substantial interest to understand in detail the effects conduction may have on the formation and structure of galaxy clusters. In particular, it is far from clear whether the temperature profile required for a local balance between cooling and conduction can naturally arise during hierarchical formation of clusters in the $`\mathrm{\Lambda }`$CDM cosmology. This question is best addressed with cosmological hydrodynamical simulations that fully account for the coupled non-linear dynamics of the gas subject to radiative cooling and thermal conduction.

In this paper, we hence develop a new numerical implementation of conduction and include it in a modern TreeSPH code for structure formation. Smoothed particle hydrodynamics (Lucy, 1977; Monaghan & Lattanzio, 1985; Monaghan, 1992) is a powerful numerical tool to investigate gas dynamics, which has found widespread application in astrophysics. In cosmological simulations, gas densities vary over many orders of magnitude, and can change rapidly as function of position and time. Classical mesh-based Euclidean approaches to hydrodynamics have difficulty adjusting to this high dynamic range unless sophisticated adaptive mesh refinements (AMR) methods are used. In contrast, the Lagrangian approach of SPH, thanks to its matter-tracing nature, guarantees good spatial resolution in high-density regions, while spending little computational time on low-density regions of space, where a coarser spatial resolution is usually sufficient.

Another useful aspect of SPH is that the equations that describe physical processes can often be translated into an SPH form rather intuitively. Below we discuss this in detail for the conduction equation, highlighting in particular how a number of practical problems with respect to stability can be overcome. The final formulation we propose is robust and explicitly conserves thermal energy, even when individual timesteps for each particle are used.

After examining idealised test problems to validate our implementation of conduction, we apply our code to realizations of the static cluster models of ZN, investigating in particular, to what degree conduction may balance cooling in these clusters, and for how long this approximate equilibrium can be maintained. We also discuss results for a first set of cosmological simulations of cluster formation. We here compare simulations that follow only adiabatic gas dynamics, or only cooling and star formation, with corresponding ones that also include conduction. As we will see, thermal conduction can lead to a substantial modification of the final thermodynamic properties of rich clusters.

The outline of this paper is as follows. In Section 2, we discuss the basics of the conduction equation, and our numerical approach for discretising it in SPH. In Section 3, we show first test runs of conducting slabs, which we also use to illustrate various issues of numerical stability. In Section 4, we consider isolated clusters, constructed with an initial equilibrium model, while in Section 5 we present a comparison of results obtained in cosmological simulations of cluster formation. Finally, we summarise and discuss our findings in Section 6.

## 2 Thermal conduction in SPH

### 2.1 The conduction equation

Heat conduction is a transport process for thermal energy, driven by temperature gradients in the conducting medium. Provided the mean free path of particles is small compared to the scale length of the temperature variation, the local heat flux can be described by

$$𝒋=\kappa \mathbf{}T,$$
(1)
where $`T(𝒓)`$ gives the temperature field, and $`\kappa `$ is the heat conduction coefficient, which may depend on local properties of the medium. For example, in the case of an astrophysical plasma, we encounter a strong dependence of $`\kappa `$ on the temperature itself.

The rate of temperature change induced by conduction can simply be obtained from conservation of energy, viz.

$$\rho \frac{\text{d}u}{\text{d}t}=\mathbf{}\mathbf{}𝒋,$$
(2)
where $`u`$ is the thermal energy per unit mass, and $`\rho `$ denotes gas density. Eliminating the heat flux, this can also be written directly in terms of $`u`$, giving the heat conduction equation in the form

$$\frac{\text{d}u}{\text{d}t}=\frac{1}{\rho }\mathbf{}\mathbf{}(\kappa \mathbf{}T)\text{.}$$
(3)
Note that the temperature is typically simply proportional to $`u`$, unless the mean particle weight changes in the relevant temperature regime, for example as a result of a phase transition from neutral gas to ionised plasma. In the thin astrophysical plasmas we are interested in, the strong temperature dependence of the Spitzer conductivity (see below) suppresses conduction in low-temperature gas heavily. For all practical purposes we can set $`u=kT/[(\gamma 1)\mu ]=c_vT`$, where $`\mu =0.588m_\mathrm{p}`$ is the mean molecular weight of a fully ionised gas with the primordial mix of helium and hydrogen, and $`c_v`$ is the heat capacity per unit mass.

### 2.2 Spitzer Conductivity

Spitzer (1962) derived the classical result for the heat conductivity due to electrons in an ionised plasma. It is given by

$$\kappa _{\mathrm{sp}}=1.31n_e\lambda _ek\left(\frac{kT_e}{m_e}\right)^{1/2},$$
(4)
where $`n_e`$ is the electron density, and $`\lambda _e`$ the electron mean free path. Interestingly, the product $`n_e\lambda _e`$ depends only on the electron temperature $`T_e`$,

$$\lambda _en_e=\frac{3^{3/2}(kT_e)^2}{4\sqrt{\pi }e^4\mathrm{ln}\mathrm{\Lambda }},$$
(5)
provided we neglect the very weak logarithmic dependence of the Coulomb logarithm $`\mathrm{\Lambda }`$ on electron density and temperature, which is a good approximation for clusters. We will set $`\mathrm{ln}\mathrm{\Lambda }=37.8`$, appropriate for the plasma in clusters of galaxies (Sarazin, 1988). The Spitzer conductivity then shows only a strong temperature dependence, $`\kappa _{\mathrm{sp}}T^{5/2}`$, and has the value

$$\kappa _{\mathrm{sp}}=8.2\times 10^{20}\left(\frac{k_\mathrm{B}T}{10\text{keV}}\right)^{5/2}\frac{\text{erg}}{\text{cm}\text{s}\text{keV}}\text{.}$$
(6)
Note that the presence of magnetic fields can in principle strongly alter the conductivity. Depending on the field configuration, it can be suppressed in certain directions, or even in all directions in cases of certain tangled fields. The field configuration in clusters is not well understood, and it is currently debated to what degree magnetic fields suppress conduction. We will assume that the modification of the conductivity can be expressed in terms of an effective conductivity, which we parameterise as a fraction of the Spitzer conductivity.

Even in the absence of magnetic fields, the Spitzer conductivity can not be expected to apply down to arbitrarily low plasma densities. Eventually, the scale length of the temperature gradient will become comparable or smaller than the electron mean free path, at which point the heat flux will saturate, with no further increase when the temperature gradient is increased (Cowie & McKee, 1977). This maximum heat flux $`j_{\mathrm{sat}}`$ is given by

$$j_{\mathrm{sat}}0.4n_\mathrm{e}k_\mathrm{B}T\left(\frac{2k_\mathrm{B}T}{\pi m_\mathrm{e}}\right)^{1/2}\text{.}$$
(7)
In order to have a smooth transition between the Spitzer regime and the saturated regime, we limit the conductive heat flux by defining an effective conductivity (Sarazin, 1988) in the form

$$\kappa _{\mathrm{eff}}=\frac{\kappa _{\mathrm{sp}}}{1+4.2\lambda _\mathrm{e}/l_T}.$$
(8)
Here $`l_TT/|\mathbf{}T|`$ is the characteristic length-scale of the temperature gradient.

### 2.3 SPH formulation of conduction

At first sight, equation (3) appears to be comparatively easy to solve numerically. After all, the time evolution generated by the diffusion equation smoothes out initial temperature variations, suggesting that it should be quite ‘forgiving’ to noise in the discretisation scheme, which should simply also be smoothed out.

In practice, however, there are two problems that make it surprisingly difficult to obtain stable and robust implementations of the conduction equation in cosmological codes. The first has to do with the presence of second derivatives in equation (3), which in standard SPH kernel-interpolants can be noisy and sensitive to particle disorder. The second is that an explicit time integration method can easily lead to an unstable integration if large local gradients arise due to noise. We will discuss our approaches to solve these two problems in turn.

A simple discretisation of the conduction equation in SPH can be obtained by first estimating the heat flux for each particle applying standard kernel interpolation methods, and then estimating the divergence in a second step. However, this method has been shown to be quite sensitive to particle disorder (Brookshaw, 1985), which can be traced to the effective double-differentiation of the SPH-kernel. In addition, this method has the practical disadvantage that an intermediate result, the heat flux vectors, need to be computed and stored in a separate SPH-loop.

It is hence advantageous to use a simpler SPH discretisation of the Laplace operator, which should ideally involve only first order derivatives of the smoothing kernel. Such a discretisation has been proposed before (Brookshaw, 1985; Monaghan, 1992), and we here give a brief derivation of it in three dimensions.

For a well-behaved field $`Y(𝒙)`$, we can consider a Taylor-series approximation for $`Y(𝒙_j)`$ in the proximity of $`Y(𝒙_i)`$, e.g.

$`Y(𝒙_j)Y(𝒙_i)`$ $`=`$ $`\mathbf{}Y|_{𝒙_i}(𝒙_j𝒙_i)+`$ (9)
$`{\displaystyle \frac{1}{2}}{\displaystyle \frac{^2Y}{x_sx_k}}|_{𝒙_i}(𝒙_j𝒙_i)_s(𝒙_j𝒙_i)_k+`$
$`𝒪(𝒙_j𝒙_i)^3\text{.}`$
Neglecting terms of third and higher orders, we multiply this with

$$\frac{(𝒙_j𝒙_i)\mathbf{}_iW(𝒙_j𝒙_i)}{|𝒙_j𝒙_i|^2}\text{,}$$
(10)
where $`W(𝒙)=W(|𝒙|)`$ is the SPH smoothing kernel. Note that we choose this kernel to be spherically symmetric and normalised to unity. The expression in equation (10) is well behaved for $`𝒙_j𝒙_i`$ under these conditions. Introducing the notations $`𝒙_{ij}=𝒙_j𝒙_i`$ and $`W_{ij}=W(𝒙_j𝒙_i)`$, we integrate over all $`𝒙_j`$ and note that

$$𝒙_{ij}\frac{𝒙_{ij}\mathbf{}_iW_{ij}}{|𝒙_{ij}|^2}d^3𝒙_j=0,$$
(11)
$$(𝒙_{ij})_s(𝒙_{ij})_k\frac{𝒙_{ij}\mathbf{}_iW_{ij}}{|𝒙_{ij}|^2}d^3𝒙_j=\delta _{sk}.$$
(12)
So the term linear in $`\mathbf{}Y`$ drops out, and the terms involving off-diagonal elements of the Hesse matrix of $`Y`$ vanish, so that the sum over the second order term simply reduces to $`\mathbf{}^2Y`$. We hence end up with

$$\mathbf{}^2Y|_{𝒙_i}=2\frac{Y(𝒙_j)Y(𝒙_i)}{|𝒙_{ij}|^2}𝒙_{ij}\mathbf{}_iW_{ij}\mathrm{d}^3𝒙_j\text{.}$$
(13)
This analytical approximation of the Laplacian can now be easily translated into an SPH kernel interpolant. To this end, we can replace the integral by a sum over all particles indexed by $`j`$, and substitute the volume element $`\text{d}^3𝒙_j`$ by its discrete SPH analogue $`m_j/\rho _j`$. The values of the field $`Y(𝒙)`$ at the particle coordinates can be either taken as the value of the intrinsic particle property that is evolved, $`Y(𝒙_i)=Y_i`$, or as a kernel interpolant of these values, $`Y(𝒙_i)=Y_i`$, where for example

$$Y_i=\underset{j}{}Y_j\frac{m_j}{\rho _j}W(𝒙_{ij}).$$
(14)
We then end up with a discrete SPH approximation of the Laplace operator in the form

$$\mathbf{}^2Y|_i=2\underset{j}{}\frac{m_j}{\rho _j}\frac{Y_jY_i}{|𝒙_{ij}|^2}𝒙_{ij}\mathbf{}_iW_{ij}.$$
(15)
We now consider how this can be applied to the thermal conduction problem, where the conductivity may also show a spatial variation. Using the identity

$$\mathbf{}\left(\kappa \mathbf{}T\right)=\frac{1}{2}\left[\mathbf{}^2(\kappa T)T\mathbf{}^2\kappa +\kappa \mathbf{}^2T\right],$$
(16)
we can use our result from equation (15) to write down a discretised form of equation (3):

$$\frac{\text{d}u_i}{\text{d}t}=\underset{j}{}\frac{m_j}{\rho _i\rho _j}\frac{(\kappa _j+\kappa _i)(T_jT_i)}{|𝒙_{ij}|^2}𝒙_{ij}\mathbf{}_iW_{ij}\text{.}$$
(17)
This form is antisymmetric in the particles $`i`$ and $`j`$, and the energy exchange is always balanced on a pairwise basis, i.e. conservation of thermal energy is manifest. Also, it is easy to see that the total entropy can only increase, and that heat always flows from higher to lower temperature.

The conductivities $`\kappa _i`$ and $`\kappa _j`$ in equation (17) are effectively arithmetically averaged. Cleary & Monaghan (1999) proposed to make the replacement

$$\frac{\kappa _i+\kappa _j}{2}\kappa _{ij}=\frac{2\kappa _i\kappa _j}{\kappa _i+\kappa _j}\text{.}$$
(18)
They showed that this ensures a continuous heat flux even in cases when the heat conductivity exhibits a discontinuity, as for example along the interface between different phases. It is clear then that this modification should also behave better when the conductivity changes extremely rapidly on small scales, as it can happen for example in ICM gas when cool particles get into direct contact with comparatively hot neighbours. Indeed, we found this symmetrisation to give numerically more robust behaviour, particularly in simulations that in addition to heat conduction also followed radiative cooling. Note that since we have $`\mathrm{min}(\kappa _i,\kappa _j)2\kappa _i\kappa _j/(\kappa _i+\kappa _j)2\mathrm{min}(\kappa _i,\kappa _j)`$, the Cleary & Monaghan average stays always close to the smaller of the two conductivities involved, to within a factor of two.

### 2.4 Numerical implementation details

We have implemented heat conduction in a new version of the massively parallel TreeSPH-code GADGET (Springel et al., 2001b), which is a general purpose code for cosmological structure formation. Unlike earlier public releases of the code, the present version, GADGET-2, uses the ‘entropy formulation’ of SPH proposed by Springel & Hernquist (2002) which conserves both energy and entropy (in regions without shocks) for fully adaptive smoothing lengths. In this formulation, an entropic function

$$A=(\gamma 1)\frac{u}{\rho ^{\gamma 1}}$$
(19)
is evolved as independent variable for each particle, instead of the thermal energy per unit mass. Note that $`A=A(s)`$ is only a function of the thermodynamic entropy $`s`$ per unit mass.

For consistency with this formalism, we need to express the heat conduction equation in terms of entropy. This is easily accomplished in an isochoric approximation, noting that the temperature can be expressed as $`T=\mu /k_BA\rho ^{\gamma 1}`$, where $`\mu `$ is the mean molecular weight. This results in

$$\frac{\text{d}A_i}{\text{d}t}=\frac{2\mu }{k_B}\frac{\gamma 1}{\rho _i^{\gamma 1}}\underset{j}{}\frac{m_j\kappa _{ij}}{\rho _i\rho _j}\left(\frac{A_j}{\rho _j^{1\gamma }}\frac{A_i}{\rho _i^{1\gamma }}\right)\frac{𝒙_{ij}\mathbf{}_iW_{ij}}{|𝒙_{ij}|^2}\text{.}$$
(20)
One problematic aspect of the heat conduction equation is that small-scale numerical noise in the temperature field can generate comparatively large heat flows, simply because this noise can involve small-scale gradients of sizable magnitude. Since we are using an explicit time integration scheme for the hydrodynamical evolution, this immediately raises the danger of instabilities in the integration. Unless extremely small timesteps or an implicit integration scheme are used, the energy exchange between two particles due to a small-scale temperature difference can become so large that the explicit time integration “overshoots”, thereby potentially reversing the sign of the temperature difference between the two particles in conductive contact. This is not only incorrect, but makes it possible for the temperature differences to grow quickly in an oscillatory fashion, causing an instable behaviour of the integration.

We have found that a good method to avoid this problem is to use a kernel interpolant for the temperature field (or entropy field) in the discretisation of the heat conduction equation (20), instead of the individual particle temperature values themselves. The interpolant represents a smoothed version of the noisy sampling of the temperature field provided by the particle values, so that on the scale of the SPH smoothing length, small-scale noise in the heat flux is strongly suppressed. Particles will still try to equilibrate their temperatures even within the smoothing radius of each particle, but this will happen at a damped rate. Heat-conduction due to temperature gradients on larger scales is unaffected however. In a later section of this work, we will explicitly demonstrate how this improves the stability of the time integration, particularly when individual and adaptive timesteps are used.

For definiteness, the interpolant we use is a smoothed version $`\overline{A}_i`$ of the entropy, defined by

$$\rho _i^\gamma \overline{A}_i=\underset{j}{}m_j\rho _j^{(\gamma 1)}A_jW_{ij}\text{.}$$
(21)
We then replace $`A_i`$ and/or $`A_j`$ on the right-hand-side of equation (20) with the interpolants $`\overline{A}_i`$ and $`\overline{A}_j`$. Note that the weighting by $`\rho ^{(\gamma 1)}`$ ensures that we obtain a value of $`\overline{A}`$ that corresponds to a smoothed temperature field, as required since conduction is driven by gradients in temperature and not entropy. However, since the density values need to be already known to evaluate this interpolant, this unfortunately requires an additional SPH loop, which causes quite a bit of computational overhead. This can in principle be avoided if the weighting with densities in equation (21) is dropped, which according to our tests appears to be sufficiently accurate in most situations.

Note that although we use the SPH interpolant for the entropy values in equation (20), we still compute the values for the particle conductivity $`\kappa `$ based on the intrinsic particle temperature and not on its smoothed counterpart.

While the above formulation manifestly conserves thermal energy, this property may get lost when individual and adaptive time steps are used, where for a given system step only a subset of the particles is evolved. There can then be particle pairs where only one particle is active in the current system step, while the other is not, so that a ‘one-sided’ heat conduction may occur that causes a (brief) violation of energy conservation during the step. While most of the resulting energy imbalance will only be a temporary fluctuation that will be compensated as soon as the ‘inactive’ particle in the pair is evolved, the very strong temperature dependence of the conductivity may produce sizable errors in this situation, particularly when coarse timestepping is used. We therefore decided to implement an explicitly conservative scheme for the heat exchange, even when adaptive and individual timesteps are used.

To this end, we define a pairwise exchange of heat energy as

$$E_{ij}\frac{2\mu }{k_B}\frac{m_im_j\kappa _{ij}}{\rho _i\rho _j}\left(\frac{A_j}{\rho _j^{1\gamma }}\frac{A_i}{\rho _i^{1\gamma }}\right)\frac{𝒙_{ij}\mathbf{}_iW_{ij}}{|𝒙_{ij}|^2}\text{.}$$
(22)
A simple translation into a finite difference scheme for the time evolution would then be

$$m_iu_i^{}=m_iu_i+\mathrm{\Delta }t_i\underset{j}{}E_{ij},$$
(23)
where $`u_i^{}`$ and $`u_i`$ can also be expressed in terms of the corresponding entropy values $`A_i^{}`$ and $`A_i`$. However, for individual and variable timesteps, only a subfraction of particles will be ‘active’ in the current system timestep. These particles have individual timesteps $`\mathrm{\Delta }t_i`$, while the ‘inactive’ particles can be formally assigned $`\mathrm{\Delta }t_i=0`$ for the step. Equation (23) then clearly does not guarantee detailed energy conservation in each system step.

We recover this property in the following way. We update the energy of particles according to

$$m_iu_i^{}=m_iu_i+\frac{1}{2}\underset{jk}{}\mathrm{\Delta }t_j(\delta _{ij}\delta _{ik})E_{jk}$$
(24)
in each system step. In practice, the double-sum on the right hand side can be simply computed using the usual SPH loop for active particles. Each interacting pair found in the neighbour search for active particle $`i`$ is simply used to change the thermal energy of particle $`i`$ by $`\mathrm{\Delta }t_iE_{ij}/2`$, and that of the neighbouring particle $`j`$ by $`\mathrm{\Delta }t_iE_{ij}/2`$ as well. Note that if particle $`j`$ is active, it carries out a neighbour search itself and will find $`i`$ itself, so that the total energy change of $`i`$ due to the presence of $`j`$ is given by $`(\mathrm{\Delta }t_i+\mathrm{\Delta }t_j)E_{ij}/2`$. Energy is conserved by construction in this scheme, independent of the values if individual timesteps of particles. If all particles have equal steps, equation (24) is identical to the form of (23). As before, in our SPH implementation we convert the final thermal energy change to a corresponding entropy change. These entropy changes are applied instantaneously for all particles at the end of one system step, so that even particles that are inactive at the current system step may have their entropy changed by conduction. The resulting low order of the time integration scheme is sufficient for the diffusion equation.

If conduction is so strong that the timescale of conductive redistribution of thermal energy is short compared to relevant dynamical timescales of the gas, the usual hydrodynamic timestep selected by the code based on the Courant-criterion may become too coarse to follow conduction accurately. We therefore introduce an additional timestep criterion in simulations that follow thermal conduction. To this end, we limit the maximum allowed timestep to some prescribed fraction of the conduction time scale, $`A/|\dot{A}_{\mathrm{cond}}|`$, namely

$$\mathrm{\Delta }t_{\mathrm{cond}}=\alpha \times \frac{A}{|\dot{A}_{\mathrm{cond}}|}\text{,}$$
(25)
where $`\alpha `$ is a dimensionless accuracy parameter. In our simulations, we typically employed a value of $`\alpha =0.25`$, which has provided good enough accuracy at a moderate cost of CPU time.

Finally, in order to account for a possible limitation of conduction by saturation, we compute the gradient of the (smoothed) temperature field in the usual SPH fashion, and then use it to compute $`l_T`$ and the saturation-limited conductivity based on equation (8).

## 3 Illustrative test problems

In this section, we verify our numerical implementation of thermal conduction with a number of simple test problems that have known analytic solutions. We will also investigate the robustness of our formulation with respect to particle disorder, and (initial) noise in the temperature field.

### 3.1 Conduction in a one-dimensional problem

We first consider two solid slabs with different initial temperatures, brought into contact with each other at time $`t=0`$. The slabs were realized as a lattice of SPH particles, with dimensions $`100\times 10\times 10`$, and an equidistant particle-spacing of $`1\text{cm}`$. To avoid perturbations of the 1D-symmetry, the simulation volume was taken to be periodic along the two short axes, while kept non-periodic along the long axis, allowing us to study surface effects, if present, on the sides of the slabs that are not in contact. Note that the test was carried out with our 3D code. All particle velocities were set to zero initially, and kept at this value to mimic a solid body. The thermal conductivity was set to a value corresponding to $`\alpha =\kappa /(c_v\rho )=1\text{cm}^2\text{s}^1`$ throughout both materials, independent of temperature.

In Figure 1, we compare the time evolution of our numerical results for this test with the analytical solution of the same problem, which is given by

$$u(z,t)=u_0+\frac{\mathrm{\Delta }u}{2}\mathrm{erf}\left(\frac{zz_m}{\sqrt{4\alpha t}}\right),$$
(26)
where $`z_m`$ gives the position of the initial difference of size $`\mathrm{\Delta }u`$ in thermal energy, and $`u_0`$ is the mean thermal energy. We see that the numerical solution tracks the analytic result very nicely.

We have also repeated the test for a particle configuration corresponding to a ‘glass’ (White, 1996), with equally good results. For a Poisson distribution of particles, we noted however a small reduction in the speed of conduction when a small number of neighbours of 32 is used. Apparently, here the large density fluctuations due to the Poisson process combined with the smallness of the number of neighbours leads to somewhat poor coupling between the particles. This effect goes away for larger numbers of neighbours, as expected. Note however that in practice, due to the pressure forces in a gas, the typical configuration of tracer particles is much more akin to a glass than to a Poisson distribution.

Next, we examine how robust our formulation is with respect to small-scale noise in the temperature field. To this end, we repeat the above test using the glass configuration, but we perturb the initial thermal energies randomly with fluctuations at an rms-level of $`\sigma =80\text{erg g}^1`$. Further, we increase the maximum timestep allowed for particles to $`\mathrm{\Delta }t=1.0\text{s}`$.

We discussed previously that we can either use the particle values of the temperatures (or entropy) in the right-hand-side of equation (22), or kernel interpolants thereof. Here we compare the following different choices with each other:

1. Basic formulation: Use particle values $`A_i`$ and $`A_j`$.
2. Smoothed formulation: Use $`\overline{A}_i`$ and $`\overline{A}_j`$.
3. Mixed formulation: Use $`A_i`$ and $`\overline{A}_j`$.

The mixed formulation (C) may at first seem problematic, because its pair-wise antisymmetry is not manifest. However, since we use equation (24) to exchange heat between particles, conservation of thermal energy is ensured also in this case. When all particles have equal timesteps, formulation (C) would correspond to an arithmetic average of (A) and (B).

In Figure 2, we compare the result obtained for these three formulations after a simulation time of 5 sec. The artificially perturbed initial conditions are shown in the top left panel. Interestingly, while all three different formulations are able to recover the analytic solution in the mean, they show qualitatively different behaviour with respect to the imprinted noise. When the intrinsic particle values for the temperatures are used (top right panel), very large pair-wise gradients on small scales occur that induce large heat exchanges. As a result, the particles oscillate around the mean, maintaining a certain rms-scatter which does not reduce with time, i.e. an efficient relaxation to the medium temperature does not occur. Note that the absolute size of the scatter is larger in the hot part of the slab. This is a result of the timestep criterion (25), which manages to hedge the rms-noise to something of order $`\alpha T`$. Reducing the timestep parameter $`\alpha `$ can thus improve the behaviour, while simultaneously increasing the computational cost significantly. For coarse timestepping (or high conductivities) the integration with this method can easily become unstable.

Formulation (B), which uses the smoothed kernel-interpolated temperature field for both particles in each pair, does significantly better in this respect (bottom left panel). However, the particle temperatures only very slowly approach the local mean value and hence the analytical solution in this case. This is because the smoothing here is quite efficient in eliminating the small-scale noise, meaning that a deviation of an individual particle’s temperature from the local mean is decaying only very slowly.

The mixed formulation (C), shown in the bottom right, obviously shows the best behaviour in this test. It suppresses noise quickly, matches the analytical solution very well, and allows the largest timesteps of all schemes we tested. In this formulation, particles which have a large deviation from the local average temperature try to equalise this difference quickly, while a particle that is already close to the mean is not ‘pulled away’ by neighbours that have large deviations. Apparently, this leads to better behaviour than for schemes (A) and (B), particularly when individual timesteps are used. We hence choose formulation (C) as our default method.

### 3.2 Conduction in a three-dimensional problem

As a simple test of conduction in an intrinsically three-dimensional problem, we consider the temporal evolution of a point-like thermal energy perturbation. The time evolution of an initial $`\delta `$-function is given by the three-dimensional Green’s function for the conduction problem,

$$G(x,y,z,t)=\frac{1}{(4\pi \alpha t)^{3/2}}\mathrm{exp}\left(\frac{x^2+y^2+z^2}{4\alpha t}\right)\text{,}$$
(27)
where $`c_v=u/T`$ is the heat capacity. We again consider a solid material, realized with a glass configuration of $`50^3`$ simulation particles with a mass of $`1\text{g}`$ each and a mean particle spacing of $`1\mathrm{cm}`$. We choose a conductivity of $`\kappa c_v=1\text{cm}^2\text{s}^1`$, as before. We give the material a specific energy per unit mass of $`u_0=1000\text{erg}\text{g}^1`$, and add a perturbation of $`10000.0\text{erg}\text{cm}^3\text{g}^1\times G(x,y,z,t_0=10\text{sec})`$. These initial conditions correspond to a $`\delta `$-function perturbation that has already evolved for a brief period; this should largely eliminate timing offsets that would arise in the later evolution if we let the intrinsic SPH smoothing wash out a true initial $`\delta `$-function.

Figure 3 shows the specific energy profile after evolving the conduction problem for an additional 10 or 20 seconds, respectively. In the two bottom panels, we show the relative deviation from the analytic solution. Given that quite coarse timestepping was used for this problem (about $`32`$ timesteps for 20 sec), the match between the theoretical and numerical solutions is quite good. Even at a radial distance of $`20\text{cm}`$, where the Gaussian profile has dropped to less than one hundredth of its central amplitude, relative deviations are at around 10 percent for $`t=20\text{sec}`$ and below 3 percent at $`t=30\text{sec}`$. Note however that at all times the numerical solution maintains a nice Gaussian shape so that the deviations can be interpreted as a small modification in the effective conductivity. We then see that at late times, when the temperature gradients are resolved by more particles, the SPH estimate of the conductivity becomes ever more accurate.

## 4 Spherical models for clusters of galaxies

Observed temperature profiles of clusters of galaxies are often characterised by a central decline of temperature, while otherwise appearing fairly isothermal over the measurable radial range. Provided the conductivity is non-negligible, there should hence be a conductive heat flux into the inner parts, which would then counteract central cooling losses. Motivated by this observation, Zakamska & Narayan (2003) (ZN henceforth) have constructed simple analytic models for the structure of clusters, invoking as a key assumption a local equilibrium between conduction and cooling. Combined with the assumption of hydrostatic equilibrium and spherical symmetry, they were able to quite successfully reproduce the temperature profiles of a number of observed clusters.

We here use the model of ZN as a test-bed to check the validity of our conduction modelling in a realistic cosmological situation. In addition, the question for how long the ZN solution can be maintained is of immediate interest, and we address this with our simulations of spherical clusters as well.

For definiteness, we briefly summarise the method of ZN for constructing a cluster equilibrium model. The cluster is assumed to be spherically symmetric, with the gas of density $`\rho (r)`$ and pressure $`P(r)`$ being in hydrostatic equilibrium, i.e.

$$\frac{1}{\rho }\frac{\text{d}P}{\text{d}r}=\frac{\text{d}\mathrm{\Phi }}{\text{d}r},$$
(28)
where the gravitational potential $`\mathrm{\Phi }`$ is given by

$$\frac{1}{r^2}\frac{\text{d}}{\text{d}r}\left(r^2\frac{\text{d}\mathrm{\Phi }}{\text{d}r}\right)=4\pi G\left(\rho _{\mathrm{DM}}+\rho \right)\text{.}$$
(29)
The dark matter density $`\rho _{\mathrm{DM}}(r)`$ is described by a standard NFW (Navarro, Frenk & White, 1996) halo, optionally modified by ZN with the introduction of a small softened dark matter core. The structure of the dark halo can then be fully specified by the virial radius $`r_{200}`$, defined as the radius enclosing a mean overdensity of 200 with respect to the critical density, and a scale length $`r_s`$.

The radial heat flux $`F`$ due to electron conduction is given by

$$F=\kappa \frac{\text{d}T}{\text{d}r},$$
(30)
where the conductivity $`\kappa `$ is taken to be a constant fraction of the Spitzer value. In order to balance local cooling losses, we require

$$\frac{1}{r^2}\frac{\text{d}}{\text{d}r}\left(r^2F\right)=j,$$
(31)
where $`j`$ is describing the local gas cooling rate. If cooling is dominated by thermal bremsstrahlung, $`j`$ can be approximated as

$$j=2.1\times 10^{27}n_e^2T^{1/2}\text{ergs cm}^3\text{s}^1.$$
(32)
The electron number density $`n_e`$ herein correlates with the gas density $`\rho `$ and depends on the hydrogen mass fraction $`X`$ we assume for the primordial matter. With a proton mass of $`m_p`$, it equals $`n_e=\rho (X+1)/(2m_p)`$. Following ZN, we use this simplification in setting up our cluster models.

Equations (28) to (32) form a system of differential equations that can be integrated from inside out once appropriate boundary conditions are specified. We adopt the values for central gas density and temperature, total dark matter mass, and scale radius determined by ZN for the cluster A2390 such that the resulting temperature profile matches the observed one well in the range $`0rR_{\mathrm{out}}`$, where $`R_{\mathrm{out}}`$ was taken to be $`2r_s`$.

Integrating the system of equations out to $`R_{\mathrm{out}}`$, we obtain a result that very well matches that reported by ZN. However, we also need to specify the structure of the cluster in its outer parts in order to be able to simulate it as an isolated system. ZN simply assumed this part to be isothermal at the temperature $`T_{\mathrm{out}}=T(R_{\mathrm{out}})`$, thereby implicitly assuming that the equilibrium condition between cooling and conduction invoked for the inner parts is not valid any more. It is a bit unclear why such a sudden transition should occur, but the alternative assumption, that equation (31) holds out to the virial radius, clearly leads to an unrealistic global temperature profile. In this case, the temperature would have to monotonically increase out to the outermost radius. Also, since the cumulative radiative losses out to a radius $`r`$ have to be balanced by a conductive heat flux of equal magnitude at that radius, the heat flux also monotonically increases, such that all the energy radiated by the cluster would have to be supplied to the cluster at the virial radius.

We nevertheless examined both approaches, i.e. we construct cluster models where we solve equations (28)-(32) for the whole cluster out to a radius of $`r_{200}`$, and secondly, we construct initial conditions where we follow ZN by truncating the solution at $`2r_\mathrm{s}`$, continuing into the outer parts with an isothermal solution that is obtained by dropping equations (30) and (31). Since conduction may still be important at radii somewhat larger than $`2r_s`$, these two models may hence be viewed as bracketing the expected real behaviour of the cluster in a model where cooling is balanced by conduction.

For the plasma conductivity, we assumed a value of $`0.3\kappa _{\mathrm{sp}}`$, as proposed by ZN in their best-fit solution for A2390. Note that this value implies that a certain degree of suppression of conduction by magnetic fields is present, but that this effect is (perhaps optimistically) weak, corresponding to what is expected for chaotically tangled magnetic fields (Narayan & Medvedev, 2001).

Having obtained a solution for the static cluster model, we realized it as 3D initial conditions for GADGET. We used $`2\times 10^5`$ gas particles with a total baryonic mass of around $`8.6\times 10^{14}\text{M}_{}`$. For simplicity, we described the NFW dark matter halo of mass $`2.4\times 10^{15}\mathrm{M}_{}`$ as a static potential. We then simulated the evolution of the gas subject to self-gravity (with a gravitational softening length of $`3\text{kpc}`$), radiative cooling, and thermal conduction, but without allowing for star formation and associated feedback processes.

By construction, we expect that thermal conduction will be able to offset radiative cooling for these initial conditions, at least for some time. This can be verified in Figure 4, where we plot the local cooling rate and conductive heating rate as a function of radius. Indeed, at time $`0.15\mathrm{Gyr}`$, after a brief initial relaxation period, the two energy transfer rates exhibit the same magnitude and cancel each other with good accuracy. This represents a nice validation of our numerical implementation of thermal conduction with the temperature-dependent Spitzer-rate.

As the simulation continues, we see that the core temperature of the cluster slowly drifts to somewhat higher temperature. A secular evolution of some kind is probably unavoidable, since the balance between cooling and conductive heating cannot be perfectly static. The cluster of course still loses all the energy it radiates, which eventually must give rise to a slow quasi-static inflow of gas, and a corresponding change of the inner structure of the cluster. Because of the different temperature dependences of cooling and conductive heating, it is also not clear that the balance between cooling and conductive heating represents a stable dynamical state. While a stability analysis by ZN and Kim & Narayan (2003) suggests that thermal instability is sufficiently suppressed in models with conduction, Soker (2003) argues that local perturbations in the cooling flow region will grow to the non-linear regime rather quickly and that a steady solution with a constant heat conduction may therefore not exist.

In any case, it is clear that the inclusion of conduction strongly reduces central mass drop out due to cooling. This is demonstrated explicitly in Figure 5, where we show the temperature and cumulative baryonic mass profiles of the cluster after a time of $`0.6\mathrm{Gyr}`$. We compare these profiles also to an identical simulation where conduction was not included. Unlike the conduction run, this cooling-only simulation shows a very substantial modification of its temperature profile in the inner parts, where the core region inside $`0.1\mathrm{Mpc}`$ becomes much colder than the initial state. We can also see that the cooling-only simulation has seen substantial baryonic inflow as a result of central mass drop-out. Within $`0.6\mathrm{Gyr}`$, the baryonic mass enclosed in a sphere with radius $`50\mathrm{kpc}`$ has increased by 86 percent in this simulation, clearly forming a cooling flow. This contrasts strongly with the run that includes conduction, where beyond a cental distance of $`0.02\mathrm{Mpc}`$, there is no significant difference in the cumulative mass profile between $`t=0`$ and $`t=0.6\mathrm{Gyr}`$.

We also checked that runs with a reduced conductivity show a likewise reduced effect, suppressing the core cooling and matter inflow to a lesser extent than with our default model with $`0.3\kappa _{\mathrm{sp}}`$.

Finally, we note that the results discussed here are quite insensitive to whether we set-up the initial conditions with an isothermal outer part, or whether we continue the equilibrium solution to the virial radius. There is only a small difference in the long term evolution of the solution.

The model presented here helps us to verify the robustness of our conduction implementation when a temperature dependent conductivity is used, but it does not account for the time dependence of conduction expected in the cosmological context. The static, spherically symmetric construction does not reflect the hierarchical growth of clusters that is a central element of currently favoured cosmologies. To be able to make reliable statements about the effects of thermal conduction on galaxy clusters, one therefore has to trace their evolution from high redshift to the present, which is best done in full cosmological simulations. Here, the hierarchical growth of clusters implies that the conductivity is becoming low early on, when the ICM temperature is low, and becomes only large at late times when the cluster forms.

## 5 Cosmological Cluster Simulations

In this section, we apply our new numerical scheme for the treatment of thermal conduction in fully self-consistent cosmological simulations of cluster formation. We will in particular compare the results of simulations with and without conduction, both for runs that follow only adiabatic gas-dynamics, and runs that also include radiative cooling of gas. While it is beyond the scope of this work to present a comprehensive analysis of the effects of conduction in cosmological simulations, we here want to investigate a set of small fiducial runs in order to further validate the stability of our method for real-world applications, and secondly, to give a first flavour of the expected effects of conduction in simulated clusters. A more detailed analysis of cosmological implications of conduction is presented in a companion paper (Dolag et al., 2004).

We focus on a single cluster, extracted from the ‘GIF’ simulation (Kauffmann et al., 1999) of the $`\mathrm{\Lambda }`$CDM model, and resimulated using the “zoomed initial conditions” technique (Tormen et al., 1997). To this end, the particles that make up the cluster at the present time are traced back to their original coordinates in the initial conditions. The Lagrangian region of the cluster identified in this way is then resampled with particles of smaller mass, and additional small-scale perturbations from the CDM power spectrum are added appropriately. Far away from the cluster, the resolution is progressively degraded by using particles of ever larger mass.

The cluster we selected has a virial mass of $`1.1\times 10^{15}\mathrm{M}_{}`$. It is the same cluster considered in the high-resolution study of Springel et al. (2001a), with our mass resolution corresponding to their S1 simulation, except that we split each of the $`450000`$ high resolution dark matter particles into a gas and a dark matter particle (assuming $`\mathrm{\Omega }_b=0.04`$), yielding a gas mass resolution of $`8.4\times 10^9\mathrm{M}_{}`$. The boundary region was sampled with an additional 3 million dark matter particles. We then evolved the cluster forward in time from a starting redshift $`z=30`$ to the present, $`z=0`$, using a comoving gravitational softening length of $`20\mathrm{kpc}`$.

We consider 4 different simulations of the cluster, each with different physical models for the gas: (1) Adiabatic gasdynamics only, (2) adiabatic gas and conduction, (3) radiative cooling and star formation without conduction, and finally, (4) radiative cooling, star formation and conduction. The simulations with cooling and star formation use the sub-resolution model of Springel & Hernquist (2003) for the multi-phase structure of the ISM (without including the optional feedback by galactic winds offered by the model). In the two runs with conduction, we adopt the full Spitzer rate for the conductivity. While this value is unrealistically large for magnetised clusters, it serves our purpose here in highlighting the effects of conduction, given also that our cluster is not particularly hot, so that effects of conductivity can be expected to be weaker than in very rich clusters.

In Figure 6, we show projections of the mass-weighted temperature, X-ray emissivity, and gas mass density for all four simulations. Each panel displays the gas contained in a box of side-length $`8.6\mathrm{Mpc}`$ centred on the cluster. Comparing the simulations with and without conduction, it is nicely seen how conduction tends to wipe out small-scale temperature fluctuations. It is also seen that the outer parts of the cluster become hotter when conduction is included.

These trends are also borne out quantitatively when studying radial profiles of cluster properties in more detail. In Figure 7, we compare the cumulative baryonic mass profile, the temperature profile, and the radial profile of X-ray emission for all four simulations. Note that the innermost bins, for $`R<30\mathrm{kpc}`$, may be affected by numerical resolution effects. Interestingly, the temperature profiles of the runs with conduction are close to being perfectly isothermal in the inner parts of the cluster. While this does not represent a large change for the adiabatic simulation, which is close to isothermal anyway, the simulation with radiative cooling is changed significantly. Without conduction, the radiative run actually shows a pronounced rise in the temperature profile in the range of $`100200\mathrm{kpc}`$, as a result of compressional heating when gas flows in to replace gas that is cooling out of the ICM in a cooling flow. Only in the innermost regions, where cooling becomes rapid, we see a distinct drop of the temperature. Interestingly, conduction eliminates this feature in the temperature profile, by transporting the corresponding heat energy from the maximum both to parts of the cluster further out and to the innermost parts. The latter effect is probably small, however, because a smooth decline in the temperature profile in the inner parts of the cluster, as present in the ZN model, does not appear in the simulation. As a consequence, a strong conductive heat flow from outside to inside cannot develop.

Conduction may also induce changes in the X-ray emission of the clusters, which we show in the bottom panel of Figure 7. Interestingly, the inclusion of conduction in the adiabatic simulation has a negligible effect on the X-ray luminosity. This is because in contrast to previous suggestions (Loeb, 2002), the cluster does not lose a significant fraction of its thermal energy content to the outside intergalactic medium, and the changes in the relevant part of the gas and temperature profile are rather modest. We do note however that the redistribution of thermal energy within the cluster leads to a substantial increase of the temperature of the outer parts of the cluster.

For the simulations with cooling, the changes of the X-ray properties are more significant. Interestingly, we find that allowing for thermal conduction leads to a net increase of the bolometric luminosity of our simulated cluster. The panel with the cumulative baryon mass profile reveals that conduction is also ineffective in significantly suppressing the condensation of mass in the core regions of the cluster. In fact, it may even lead to the opposite effect. We think this behaviour simply occurs because a temperature profile with a smooth decline towards the centre, which would allow the conductive heating of this part of the cluster, is not forming in the simulation. Instead, the conductive heat flux is pointing primarily from the inside to the outside, which may then be viewed as an additional “cooling” process for the inner cluster regions.

It is interesting to note that in spite of the structural effects that thermal conduction has on the ICM of the cluster, it does not affect its star formation history significantly. In the two simulations that include cooling and star formation, the stellar component of the mass profile (Figure 7) does not show any sizeable difference between the run including thermal conduction and the one without it. This does not come as surprise as $`90\%`$ of the stellar content of the cluster has formed before a redshift of $`z=0.85`$. At these early times, the temperature of the gas in the protocluster was much lower, such that conduction was unimportant. In fact, for that reason, the stellar mass profiles for both simulation runs coincide.

In summary, our initial results for this cluster suggest that conduction can be important for the ICM, provided the effective conductivity is a sizable fraction of the Spitzer value. However, the interplay between radiative cooling and conduction is clearly complex, and it is presently unclear whether temperature profiles like those observed can arise in self-consistent cosmological simulations. We caution that one should not infer too much from the single object we examined here. A much larger set of cluster simulations will be required to understand this topic better.

## 6 Conclusions

Hot plasmas like those found in clusters of galaxies are efficiently conducting heat, unless electron thermal conduction is heavily suppressed by magnetic fields. Provided the latter is not the case, heat conduction should therefore be included in hydrodynamical cosmological simulations, given, in particular, that conduction could play a decisive role in moderating cooling flows in clusters of galaxies. Such simulations are then an ideal tool to make reliable predictions of the complex interplay between the nonlinear processes of cooling and conduction during structure formation.

In this paper, we have presented a detailed numerical methodology for the treatment of conduction in cosmological SPH simulations. By construction, our method manifestly conserves thermal energy, and we have formulated it such that it is robust against the presence of small-scale temperature noise. We have implemented this method in a modern parallel code, capable of carrying out large, high-resolution cosmological simulations.

Using various test problems, we have demonstrated the accuracy and robustness of our numerical scheme for conduction. We then applied our code to a first set of cosmological cluster formation simulations, comparing in particular simulations with and without conduction. While these results are preliminary, they already hint that the phenomenology of the coupled dynamics of radiative cooling and conduction is complex, and may give rise to results that were perhaps not anticipated by earlier analytic modelling of static cluster configurations.

For example, we found that conduction does not necessarily reduce a central cooling flow in our simulations; the required smoothly declining temperature profile in the inner cluster regions does not readily form in our cosmological simulations. Instead, the profiles we find are either flat, or have a tendency to slightly rise towards the centre, akin to what is seen in the cooling-only simulations. In this situation, conduction may in fact lead to additional cooling in certain situations, by either transporting thermal energy to the outer parts, or by modifying the temperature and density structure in the relevant parts of the cluster such that cooling is enhanced. This can then manifest itself in an increase of the bolometric X-ray luminosity at certain times, which is actually the case for our model cluster at $`z=0`$. Interestingly, we do not find that our cluster loses a significant fraction of its thermal heat content by conducting it to the external intergalactic medium.

In a companion paper (Dolag et al., 2004), we analyse a larger set of cosmological cluster simulations, computed with much higher resolution and with more realistic sub-Spitzer conductivities. This set of cluster allows us to investigate, e.g., conduction effects as a function of cluster temperature and the influence of conduction on cluster scaling relations. While our first results of this work suggest that conduction by itself may not resolve the cooling-flow puzzle, it also shows that conduction has a very strong influence on the thermodynamic state of rich clusters if the effective conductivity is a small fraction of the Spitzer value or more. In future work, it will hence be very interesting and important to understand the rich phenomenology of conduction in clusters in more detail.

## Acknowledgements

The simulations of this paper were carried out at the Rechenzentrum der Max-Planck-Gesellschaft, Garching. K. Dolag acknowledges support by a Marie Curie Fellowship of the European Community program “Human Potential” under contract number MCFI-2001-01227.


format: text

text:
# Possible fractal structure of exact universal exchange-correlation potential

## Acknowledgements

The author would like to thank VEGA grant 1/9179/02 for the financial support.


text:
# Spectroscopic imaging of single atoms within a bulk solid

## Abstract

The ability to localize, identify and measure the electronic environment of individual atoms will provide fundamental insights into many issues in materials science, physics and nanotechnology. We demonstrate, using an aberration-corrected scanning transmission electron microscope, the spectroscopic imaging of single La atoms inside CaTiO3. Dynamical simulations confirm that the spectroscopic information is spatially confined around the scattering atom. Furthermore we show how the depth of the atom within the crystal may be estimated.

Detection and measurement of the response of individual atoms has become a challenging issue to provide new insight into many fields in condensed matter and nanoscale sciences. Distributions of isolated atoms deeply modify the physical properties of many of the technologically most relevant and scientifically most interesting materials. Therefore, analytical techniques capable of probing single-atom identity and location are increasingly in demand. Here we show how the aberration-corrected scanning transmission electron microscope (STEM) allows not only the imaging of individual atoms inside a crystal, but their spectroscopic identification, with spatial resolution at the atomic scale. This substantial improvement in sensitivity opens up the possibility of probing the electronic environment of a single atom. Furthermore, by comparing signals from columns adjacent to that containing the single atom with dynamical simulations, the depth of the atom in the crystal can be estimated.

Much of the reported work on the identification and imaging of single atoms has been achieved by indirect techniques such as image simulation, including complex reconstruction of phase images with through focal series restoration Jia et al. (2003); Meyer et al. (2000). The STEM provides direct images of individual atoms as first demonstrated many years ago by the imaging of individual U and Th atoms on a thin carbon film Crewe et al. (1970). More recently, at a resolution of $`0.13`$ nm, the STEM has successfully imaged single atoms on and within a variety of materials, including catalysts Nellist and Pennycook (1996) and semiconductors Crewe et al. (1970); Voyles et al. (2002); Lupini and Pennycook (2003). Although the optics of the STEM is ideally suited to simultaneous imaging and spectroscopy, very few results have been reported that combine direct imaging with spectroscopic identification of individual atoms. Electron energy loss spectroscopy (EELS) signals also display high spatial resolution. Sub-unit cell analysis was first achieved by Spence and Lynch Spence and Lynch (1982) at 1.1 nm resolution. Single Gd atoms were for the first time identified by Suenaga and coworkers Suenaga et al. (2000) by combining phase contrast TEM imaging with EELS mapping in the STEM. However, their EELS resolution was only 0.6 nm, sufficient to resolve their widely spaced Gd atoms but insufficient to resolve the atomic structure of most materials. Scanning probe microscopies have also succeeded in achieving atomic-resolution spectroscopic identification of single atoms on surfaces, but they cannot probe individual atoms within the bulk environment.

Recent advances in aberration-corrected STEM provide the perfect scenario to address such a problem Dellby et al. (2001). A Z-contrast image is obtained by collecting, point by point, the Rutherford-scattered electrons at high angles using a high angle annular dark field (HAADF) detector. Energy loss electrons passing through the central hole in the detector can be collected simultaneously with the HAADF image. This correlation allows EELS to be performed with atomic resolution Browning et al. (1993); Duscher et al. (1998); Kaiser et al. (2003), limited only by the size of the probe Rafferty and Pennycook (1999); Allen et al. (2003a). Recently, by correcting the aberrations of the STEM probe-forming lens, a sub-Ångstrom beam has been demonstrated on a VG Microscopes HB501UX field emission STEM operating at 120 kV which is equipped with a Nion aberration corrector Batson et al. (2002). In a similar microscope, operated at 100 kV, a spatial resolution close to 0.11 nm is routinely achieved, and in this work we show how this enables us not only to image an individual atom within a crystal but to identify it spectroscopically.

The superlattice sample used for this study was a stack of CaTiO3 and LaxCa1-xTiO3 layers with a very low concentration of La dopants, prepared by pulsed laser deposition at $`800`$ C, in an O2 pressure of $`75`$ mTorr, with a laser energy of $`3`$ J/cm2 at $`248`$ nm and a repetition rate of $`10`$ Hz. One commercial CaTiO3 target and a specially prepared LaxCa1-xTiO3 target with $`x=0.04\pm 0.002`$ were used. The sample is epitaxial and highly ordered, and contains reference layers of La0.04Ca0.96TiO3 approximately 2.3 nm thick grown using $`100`$ shots of the La-doped target. Layers with correspondingly lower La concentrations were grown by using fewer shots, and were correspondingly thinner. From this calibration we calculate that a layer grown using just a single shot results in a La concentration of 1 atom per $`60`$ nm2. Assuming a typical TEM cross-section specimen to be about $`10`$ nm thick, on average, $`1`$ La atom should be found every 16 $`\pm `$ 4 unit-cell columns. Simple statistical arguments yield a probability of only 3% that any given La-containing column contains more than one La atom. Such cross-section specimens were prepared by conventional methods. To minimize surface damage the final cleaning was performed by ion milling at 0.5 kV. Plasma cleaning was also used to prevent sample contamination.

Figure 1 shows a high resolution Z-contrast image of a La0.002Ca0.998TiO3 layer adjacent to a thicker La0.04Ca0.96TiO3 layer. The higher atomic number of La, as opposed to Ca, O, or Ti, makes single La atoms appear as bright spots in this environment, although this brightness is dependent on atom depth Nellist and Pennycook (2000); Nakamura et al. (1997). Figure 1(b) shows the intensity traces across the atom. While we cannot completely exclude the possibility that any particular column contains two La atoms, many other single La atoms were imaged at different positions of the sample, in order to get better statistics, most giving an increase in columnar intensity by $`10\%`$. They are seen to sit on the Ca site of the LaxCa1-xTiO3 perovskite structure, easily identifiable due to the fact that Ca columns appear darker than TiO columns.

Identification of the atoms was carried out using a McMullan design charge-coupled-device parallel-detection EELS system McMullan et al. (1990). This system does not provide optimized coupling from the aberration-corrected beam, having a collection efficiency of only 8%. \[Collection efficiency is calculated as $`\beta ^2/(\alpha ^2+\theta _\mathrm{E}^2)`$ where $`\beta `$ = 7 mrad is the spectrometer entrance aperture, $`\alpha `$ = 25 mrad is the probe forming aperture and $`\theta _E=\mathrm{\Delta }E/2E_0`$ is the characteristic angle of inelastic scattering for an energy loss $`\mathrm{\Delta }E`$ and incident energy E0.\] Nevertheless, by scanning the STEM probe directly over a highly magnified image of an isolated La atom (0.1 nm $`\times `$ 0.1 nm in size as shown in figure 2), a clear EELS signature is obtained. In figure 2, spectrum 3 shows two distinct peaks that correspond to the La M4,5 lines at 832 and 849 eV. The intensity is low, as one would expect from the small number of counts associated with the excitations of a single atom, but clearly observable. The total number of M4,5 excitations for this spectrum is of the order of 104. To be certain that a single La atom can be detected, ten different bright spots were analyzed and each showed characteristic La lines with similar intensity. No meaningful beam damage was observed during the acquisition times employed, which were of the order of tenths of seconds.

There has been much discussion about detection limits and the role that the delocalization of the underlying ionization process plays in limiting the resolution of the STEM image Spence and Lynch (1982); Rafferty and Pennycook (1999); Allen et al. (2003a); Kohl and Rose (1985); Muller and Silcox (1995); Holbrook and Bird (1995). To analyze the localization of the EELS signal we have taken spectra with the probe positioned on nearby columns as indicated on figure 2(a). Given the fact that neighboring atomic columns (2 and 4) correspond to different atomic species (TiO and O columns) the channeling and therefore the localization of the signal would be different. As depicted in figure 2(b), the EELS signal attributed to the La single atom is clearly localized on the atomic column containing this ion. When the electron beam is placed on adjacent columns the signal is substantially reduced, to only $`10\pm 5`$% on the neighboring TiO column (labeled 2), and to $`20\pm 5`$% on the neighboring O site (labeled 4), only 0.19 nm away. This clearly indicates that the La can be localized within significantly better than 0.2 nm, which is much larger than the impact parameter for La-M excitation Rafferty and Pennycook (1999) and therefore suggests that the residual intensity from neighboring columns is almost certainly the effect of channeling Allen et al. (2003b).

To explore the role of channeling, dynamical diffraction calculations were performed. Figure 3(a) shows the simulated La signal strength for an aberration-free probe situated on the Ca column, O column, and TiO column as a function of the impurity depth. A simplified cubic model for the structure suffices for this estimate, but subsequent calculations are based upon the distorted perovskite structure Buttner and Maslen (1992). The contribution on the Ca column initially increases due to the electrostatic attraction to the column (channeling), but then decreases as elastic and inelastic scattering mechanisms broaden the probe. The contribution on the O column rises gradually; simulations of the electron density make clear that this is due to probe spreading. From this figure it is estimated that a $`20\pm 5\%`$ signal on the O column is obtained with an impurity depth of around 10 nm. The lack of any appreciable signal from the TiO column is an expression of the channeling strength of this column. To explain the $`10\pm 5`$% signal measured it must be appreciated that in the experiment the probe was not situated upon a single point but raster scanned across a highly magnified image. This introduces a finite, square scan region, of approximately 0.1 nm side length. Figures 3(b) and 3(c) show the HAADF images simulated Allen et al. (2003b); Findlay et al. (2003) for a 15 nm thick crystal using an aberration-free probe and an aberration-balanced probe respectively. The latter uses parameters measured from autotuning which, though recorded on a different occasion, can be regarded as typical. The latter image has been convolved with a Gaussian of half-width 0.04 nm, accounting for position instability of the probe and finite source size Nellist and Rodenburg (1994). EELS images were simulated using mixed dynamic form factor theory Findlay et al. (2003); Oxley and Allen (1998).

Figures 3(d) and 3(e) show EELS images corresponding to the parameters used in figures 3(b) and 3(c). The asymmetry of the structure, particularly whereby O columns distorted towards the Ca column provide a stronger signal than those distorted away from it, is evident. Figure 3(f) shows the scan regions and the signal strengths therein as a percentage of that on the Ca column. In these simulations the impurity was set at a depth of 12.2 nm. Note that the O column signal strength is on average higher than that of the experiment, and the TiO column signal strength is on average lower. For the second probe there are columns which agree with the experimental results to within the error bars. While the parameters used may not correspond precisely to those characterizing the experiment, the simulations have elucidated the mechanism by which the signals are obtained. The simulations support the localization of the ionization interaction and demonstrate that the dynamical propagation, in transferring some electron intensity to adjacent columns, is responsible for the signals obtained. The depth dependence seen in figure 3(a) is suggestive. It will be interesting to investigate how accurately the depth of an impurity atom may be determined in this fashion, particularly using larger probe-forming apertures and variable defocus.

In summary, an individual atom within a bulk solid has been identified spectroscopically with spatial resolution at the atomic level. This not only allows the detection and identification of single atoms, but also measurement of their electronic properties. Through dynamical simulations of probe spreading, an indication of atom depth is also obtained. Furthermore, the local formal oxidation state can be determined from the ratio of white line features present in L edges giving direct information about the magnetic moment and the valence of the atom. The present statistics are sufficient for such an analysis, although in this case the La ion is known to have a fixed 3+ valence state. However, optimum coupling of the aberration-corrected beam into the spectrometer would not only improve collection efficiency by over an order of magnitude, but also allow correction of spectrometer aberrations to give improved energy resolution Dellby et al. (2003). This would allow better comparison of local electronic structure with first-principles calculations Buczko et al. (2000).

There are many situations where the valence state of localized impurity atoms is variable, and it is very often the critical factor controlling macroscopic properties. Examples include band-bending effects at grain boundaries Kim et al. (2001); Hilgenkamp and Mannhart (2002), charge ordering in perovskite-based oxides Renner et al. (2002), charge-modulation superlattices Ohtomo et al. (2002), and molecular species in nanotubes Smith et al. (2000); Suenaga et al. (2003). We anticipate single atom analysis will open up many other areas of science where isolated atoms control macroscopic phenomena, including the origin of ductility or embrittlement in structural alloys, the nature of trap states at semiconductor interfaces, luminescent quantum efficiency in nanocrystals and active sites in heterogeneous catalysts.

The single atom represents the smallest quantum of matter that retains the characteristics of the material through its local electronic environment. The ability to probe the electronic structure of a material at this level represents the ultimate advance for understanding the atomic origins of materials properties. This work opens up this possibility.

M. Varela acknowledges fruitful discussions with R. Sanchez. Research of M. Varela is performed as a Eugene P. Wigner Fellow and staff member at ORNL. This research was sponsored by the Laboratory Directed Research and Development Program of ORNL, managed by UT-Battelle, LLC, for the U.S. Department of Energy under Contract No. DE-AC05-00OR22725 and by appointment to the ORNL Postdoctoral Research Program administered jointly by ORNL and ORISE. M. P. Oxley and L. J. Allen acknowledge support by the Australian Research Council.


text:
# Search for Southern-Hemisphere Giant Radio Galaxies with SALT

## 1. Introduction

Radio galaxies (RGs) represent the largest single objects in the Universe. Powered by an active galactic nucleus (AGN) jets emerge from the central engine, which very likely is a super-massive black hole accreting matter surrounding it. There is a huge range of linear extent of the RGs: from less than $`10^2`$ pc – Gigahertz-Peaked Spectrum (GPS), through $`10^2`$$`10^4`$ pc – Compact Steep Spectrum (CSS), $`10^4`$$`10^6`$ pc – normal-size sources, up to greater than 1 Mpc – Giant Radio Galaxies (GRG). The three largest GRGs, recognized up to now, are shown in Fig. 1.

Although giant-size radio sources are very rare among RGs, from many years they have been of a special interest for several reasons. Their very large angular size on the sky give an excellent opportunity for the study of radio source physics. They are also very useful to study the density and evolution of the intergalactic and intracluster environment. One of the key issues of the current research is attempt to trace an individual evolution of RGs. Is there a single evolutionary scheme governing the linear extent of radio sources, or do small and large sources evolve in a different way? To answer this question, in a number of papers, both theoretical and observational, attempts were undertaken to recognize factors which may differentiate giants from normal-size sources. It seems that there is no a single factor responsible for the size of classical radio sources; the large size of GRGs probably results from a combination of different factors like: age of a source, jet power, density of environment, etc. Still very limited number of well studied GRGs is a reason of that uncertainty. Therefore the phenomenon of GRG is still open for a further research.

## 2. Search for New Giants

During the IAU Symposium No. 199 (December 1999) Machalski & Jamrozy (2002) presented an evidence that only a very small fraction of expected faint GRGs of Fanaroff-Riley (1974) type II (FRII) was detected at that time. In order to find those missed giant sources we inspected the radio maps available from the large radio surveys: NVSS (Condon et al., 1998) and the first part of FIRST (Becker et al., 1995). The maps of these surveys, made with two different angular resolution (45′′ and 5′′, respectively) at the same observing frequency of 1.4 GHz, allowed (i) an effective removal of confusing sources, (ii) a reliable determination of morphological type of the giant candidate, and (iii) a detection of the compact radio core necessary for the proper identification of the source with its parent optical object. As the result we selected a sample of 36 GRG candidates (cf. Machalski et al., 2001). In order to identify their host galaxy, to determine its distance and other physical properties, we have carried out several radio and optical observations of the sample sources. In particular, we already made optical spectroscopy and got redshift for 17 out of 36 galaxies (spectroscopic redshifts of 5 sample galaxies were available prior our research). Out of 22 galaxies, 19 host giant radio sources. In the meantime, similar efforts have been undertaken by Schoenmakers et al. (2001) and Lara et al. (2001). Owing to the above studies, the statistics of giant radio galaxies is enlarged. The numbers of FRII-type GRGs, expected from our population analysis (Machalski & Jamrozy 2002), are recalled in Table 1 and compared with the observed numbers. The observed numbers denoted by an asterisk refer to the data available in 1999, while other are from the beginning of the year 2003.

Two examples of GRGs from our sample are shown in Fig. 2. J1343+3758 with the linear size of 3.14 Mpc has appeared to be the third largest source in the universe (Machalski & Jamrozy 2000), while J1604+3438 represents a very rare type of AGN – a so-called double-double RG (cf. Schoenmakers et al. 2000) which shows two pairs of lobes likely originating from an old and a new cycle of activity. Low-resolution optical spectra of host galaxies of these two giant radio sources are shown in Fig. 3. Some of the above data are used to constrain the existing analytical models of a dynamical evolution of FRII-type radio sources (Machalski et al. 2003; Chyży et al. 2003).

## 3. Why Southern Hemisphere

Our investigations of the new giant radio sources are in progress. However, we would like to extend them on GRGs on the southern sky. There are several scientific reasons for such a project, and the main of them are:

* All of the recent systematic search for new giants (Lara et al. 2001, Machalski et al. 2001, Schoenmakers et al. 2001) were performed on the northern sky. Furthermore, only about 17% of the presently known GRGs have negative declinations, and gross of them are high flux density ($`S>`$0.5 Jy) nearby objects. Therefore, one can expect a large number of undetected yet GRGs on the southern hemisphere very useful for a further enlargement of their still unsatisfactory statistics.
* The development of astronomical high-technology facilities i.e. the existing and planned large optical telescopes on the South (VLT, SALT) is very rapid. Therefore, it should be easy to attain the redshift of new GRG hosting galaxies which is the crucial observational parameter for determination of all physical parameters of the radio sources like their distance, projected linear size, volume of their lobes or cocoon, luminosity, etc. The above needs low-resolution spectroscopic observations of usually faint optical counterparts (which in many cases have very low apparent magnitudes $`R21^m23^m`$) in a reasonably short time.
* There is a high probability that the planned powerful radio interferometer – the Square Kilometer Area (SKA) will be located in the South (there are two southern site candidates for its possible location, i.e. Australia or South Africa). It would give the opportunity to detect and study low surface-brightness GRGs with a very high angular resolution and sensitivity in the future. This, in turn, would allow to recognize a very last stage of the dynamical evolution of classical double radio galaxies and learn a typical lifetime and end of that physical process.

## References

Becker, R.H., White, R.L., & Helfand, D.J. 1995, ApJ, 450, 559

Chyży, K., Jamrozy, M., Kleinman, S.J., Krzesiński, J., Machalski, J., McMillan, Nitta, A., Serafimovich, N., & Zola, S. 2003, Baltic Astronomy, (Proceedings of the JENAM 2003, astro-ph/0310606)

Condon, J.J., Cotton, W.D., Greisen, E.W., Yin, Q.F., Perley, R.A., Taylor, G.B., & Broderick, J.J. 1998, AJ, 115, 1693

Fanaroff, B.L., & Riley, J.M. 1974, MNRAS, 167, 31

Lara, L., Cotton, W.D., Feretti, L., Giovannini, G., Marcaide, J.M., Marquez, I., & Venturi, T. 2001, A&A, 370, 409

Machalski, J., Chyży, K., & Jamrozy, M. 2003, submitted for MNRAS, (astro-ph/0210546)

Machalski, J., & Jamrozy, M. 2002, in IAU Symposium 199, The Universe at Low Radio Frequencies, ed. A. P. Rao, G. Swarup, & Gopal-Krishna (San Francisco: ASP), 203

Machalski, J., Jamrozy, M., & Zola, S. 2001, A&A, 371, 445

Machalski, J., & Jamrozy, M. 2000, A&A, 363, L17

Schoenmakers, A.P., de Bruyn, A.G., Röttgering, H.J.A., & van der Laan, H. 2001, A&A, 374, 861

Schoenmakers, A.P., de Bruyn, A.G., Röttgering, H.J.A., van der Laan, H., & Kaiser, C.R. 2000, MNRAS, 315, 371


text:
# A Lie Group without universal covering

## Abstract.

We present an example of a disconnected Lie group for which there is no universal covering (as Lie group).

A connected and locally pathwise connected topological space $`X`$ with a base point $`x`$ admits a universal covering $`\stackrel{~}{X}`$ which is constructed as space of equivalence classes of curves starting at the given base point. The constant curve with value $`x`$ defines a base point $`\stackrel{~}{x}\stackrel{~}{X}`$. Continuous maps between two such spaces mapping the base point of the source space to the base point of the target space lift uniquely to a continuous map between the respective universal coverings mapping basepoint to basepoint.

Using this, for every connected Lie group the universal covering admits a structure as Lie group, because the maps defining the group structure lift uniquely. The neutral element of the Lie group is the canonical choice for the basepoint.

However, the pictures changes if we regard disconnected spaces. In order to construct a universal covering for a disconnected locally pathwise connected topological space we need a basepoint for each connected component, because we need to treat each connected component separately. But for a Lie group only the connected component containing the neutral element $`e`$ has a canonical choice of a basepoint, namely $`e`$. For the other connected components there is no natural choice. As a consequence, we may construct a universal covering as a manifold for a disconnected Lie group, but there is no natural way to put a group structure on this manifold, because there is no canonical way to lift the map defining the group structure.

In fact, there is an example, where there does not exist any compatible group structure.

###### Proposition 1.

Let $`\stackrel{~}{X}=_2\times _2\times `$ and $`X=_2\times _2\times (/)`$ as manifolds and $`\pi :\stackrel{~}{X}X`$ the natural projection induced by the natural map $`/`$. (Here $`_2=/2`$.)

Let $`X`$ be endowed with the following group structure:

$$(a,b,v)(c,d,w)=(a+c,b+d,v+w+\frac{1}{2}(ad))$$

with $`a,b,c,d_2`$ and $`v,w/`$.

Then there does not exist any Lie group structure on $`\stackrel{~}{X}`$ for which $`\pi `$ is a group homomorphism.

###### Proof.

Assume the contrary. We begin by observing that the connected component $`X^0`$ of $`X`$ containing $`e`$ is central in $`X`$. It follows that $`ghg^1h^1\pi ^1(e)`$ for all $`g\stackrel{~}{X}`$, $`h(\stackrel{~}{X})^0`$. Since $`\zeta :(g,h)ghg^1h^1`$ is continuous, $`(\stackrel{~}{X})^0`$ is connected, $`\pi ^1(e)`$ is discrete and $`\zeta (g,e)=e`$ it follows that $`\zeta (g,h)=e`$ for all $`g\stackrel{~}{X}`$, $`h(\stackrel{~}{X})^0`$. In other words, $`(\stackrel{~}{X})^0`$ is central.

Now let $`g\stackrel{~}{X}(\stackrel{~}{X})^0`$. Then $`g^2(\stackrel{~}{X})^0`$, because $`\stackrel{~}{X}/(\stackrel{~}{X})^0X/X^0_2\times _2`$. Since $`(\stackrel{~}{X})^0`$ is central and $`(\stackrel{~}{X})^0(,+)`$, it follows that every connected component of $`\stackrel{~}{X}`$ contains a unique element of order two.

Let $`g\stackrel{~}{X}`$. Then conjugation by $`g`$ stabilizes each of the connected components, because $`\stackrel{~}{X}/(\stackrel{~}{X})^0`$ is commutative. Therefore conjugation by $`g`$ must stabilize the unique element of order two for every connected component. Furthermore, conjugation by $`g`$ acts trivially on the central subgroup $`(\stackrel{~}{X})^0`$. The elements of order two together $`(\stackrel{~}{X})^0`$ generate the whole group $`\stackrel{~}{X}`$. Hence conjugation by $`g`$ is entirely trivial. Since $`g`$ was arbitrary, we have deduced that $`\stackrel{~}{X}`$ is commutative. But this is in contradiction to the fact that $`X`$ is not commutative. ∎