text:
# Cold filaments in galaxy clusters: effects of heat conduction

## 1 Introduction

Intergalactic space within clusters of galaxies contains enormous quantities of gas that is at the virial temperature. In about three quarters of all rich clusters, the temperature falls by a factor of 2 to 3 within a few tens of kiloparsecs of the cluster centre – such clusters are said to possess a ‘cooling flow’. Nearby examples are the Virgo, Hydra and Perseus clusters.

Optical emission-line filaments are observed around NGC 1275 at the centre of the Perseus cluster (Minkowski 1957; Lynds 1970; Conselice, Gallagher & Wyse 2001, and reference therein). Similar structures of extended $`\mathrm{H}\alpha `$ emission are seen near the centres of other cooling-flow clusters (e.g., Virgo, Hydra, A2597, A2052, A1795; see Heckman 1981; Hu, Cowie & Wang 1985; Johnstone, Fabian & Nulsen 1987; Heckman et al. 1989). The origin of these filaments and the mechanism responsible for their ionization have been extensively studied, but definitive conclusions have not been reached. Two formation scenarios have been proposed. In the first, the filaments form as the intracluster medium (ICM) cools (e.g., Fabian & Nulsen 1977, Cowie, Fabian & Nulsen 1980). In the other scenario, filaments are produced through cosmic infall of either cold gas or gas-rich galaxies (e.g., Soker, Bregman & Sarazin 1991; Baum 1992; Sparks 1992). The fact that filaments are observed only in clusters with cooling flows has often been considered evidence for the first scenario, but by no means all cooling-flow clusters have detectable optical emission lines (e.g., A2029; see Fabian 1994). Some observational data favour the infall scenario. For example, Sparks, Macchetto & Golombek (1989) showed that filaments contain dust with normal Galactic extinction properties, which filaments formed by condensation of the ICM should not (see also Donahue & Voit 1993). In addition, the filaments appear to be dynamically disturbed systems, suggesting that they form through violent processes, such as merging, rather through the quiescent condensation of the X-ray emitting ICM.

An important related question is, what mechanism is responsible for heating the filaments? Many possible ionization sources have been considered: Active Galactic Nuclei (AGN); shocks; massive stars; radiation from the ICM; magnetic reconnection (e.g., Heckman et al. 1989). None of these ionization mechanisms, individually, can satisfactory account for the observed emission-line ratios and intensities, so possible combinations of these processes have been explored (e.g., Sabra, Shields & Filippenko 2000). However, there are indications that the heating source is physically associated with the filaments (Conselice et. al 2001 and references therein). Thus, it is likely that stellar photoionization makes an important contribution to powering their optical emission-lines, since young stars are systematically found in the filaments (e.g., McNamara, O’Connell & Sarazin 1996).

The existence of cool filaments embedded in hot ambient plasma has implications for the problem of the plasma’s thermal conductivity, which must be moderated by the intracluster magnetic field. If the conductivity were high, one would expect the filaments to evaporate rapidly. Thus, in principle, their very existence limits the conductivity in the same way that of cold fronts in the ICM does (Markevitch et al. 2000, 2003). Given our very limited understanding of the connectedness of the intergalactic magnetic field, and the possibility that thermal conductivity contributes significantly to the global dynamics of cooling flows (Narayan & Medvedev 2001, Voigt & Fabian 2003 and references therein), observational constraints on thermal conductivity are valuable.

In this paper we address two main questions:

1. Why are cold filaments observed only in cooling-flow clusters? and why only in the cluster core?
2. To what degree does the existence of these filaments constrain the thermal conductivity of the ICM?

The competition between evaporation and condensation when a cool body of plasma is immersed in a hotter ambient medium has been studied extensively in the literature (Graham & Langer 1973; Cowie & McKee 1977; McKee & Cowie 1977; Balbus & McKee 1982; Giuliani 1984; Draine & Giuliani 1984; Balbus 1986; Böhringer & Hartquist 1987; Böhringer & Fabian 1989; McKee & Begelman 1990; Ferrara & Shchekinov 1993). Most effort has been devoted either to spherical clouds or to planar inhomogeneities, although Cowie & Songaila (1977) considered spheroidal clouds. Since the morphology of observed filaments suggests that cylindrical symmetry should provide a useful idealization, we revisit the problem for this case, and discuss the implications of our results for real systems.

In any geometry, clouds below some critical size are evaporated, while clouds above this size grow by condensation of the hot medium. We determine the critical size that separates these two regimes by finding the steady-state solution in which conductivity exactly balances radiative losses. The typical size of observed filaments is likely to be comparable to the critical size: much smaller filaments are rapidly evaporated away, while filaments above this critical value grow slowly, at a rate determined by the cooling time. Significantly larger filaments are expected to be very rare: in principle, if they formed they could fundamentally change the nature of the system by eliminating the cooling hot ambient plasma.

In Section 2 we derive the critical cylindrical solution. In Section 3 we discuss the application of the model to cluster filaments. Our conclusions are in Section 4.

## 2 Radiatively stabilized cylindrical clouds

### 2.1 Model and equations

We consider a long cylindrically symmetric cloud in which there is perfect balance between conductive heat flow towards the symmetry axis and radiative losses, so the plasma flows neither inwards nor outwards (see McKee & Cowie 1977). The temperature on the axis is assumed to be $`T10^4\mathrm{K}`$, while far from the origin $`T`$ asymptotes to a value $`10^610^8\mathrm{K}`$. With these assumptions, and neglecting the effects of gravity, the energy equation reads

$$𝐪=n_e^2\mathrm{\Lambda }(T),$$
(1)
where $`𝐪`$ is the heat flux, $`n_e`$ is the electron number density, and $`\mathrm{\Lambda }(T)`$ is the normalized cooling function. In the unsaturated thermal conduction regime111In the Appendix we discuss the possible effects of saturation., the heat flux is given by

$$𝐪=\kappa (T)T,$$
(2)
with $`\kappa (T)`$ thermal conductivity. For a hydrogen plasma, in the absence of magnetic fields, Spitzer (1962) derives $`\kappa (T)=\kappa _0T^{5/2}`$, with $`\kappa _01.84\times 10^5(\mathrm{ln}\mathrm{\Lambda })^1\mathrm{erg}\mathrm{s}^1\mathrm{cm}^1\mathrm{K}^{7/2}`$, where $`\mathrm{ln}\mathrm{\Lambda }`$ is the Coulomb logarithm, which is only weakly dependent on $`n_e`$ and $`T`$ (in the following we assume $`\mathrm{ln}\mathrm{\Lambda }=30`$). Confinement of thermal electrons by the intracluster magnetic field reduces the actual conductivity to

$$\kappa (T)=f\kappa _0T^{5/2},$$
(3)
where the parameter $`f1`$ is the heat conduction suppression factor (e.g., Binney & Cowie 1981, Böhringer & Fabian 1989).

At least in some clusters, observations indicate that emission-line nebulae are in pressure equilibrium with the surrounding hot plasma (e.g., Heckman et al. 1989). So we require that the pressure $`P=n_eT`$ is constant throughout the interface, and we use this relation to eliminate $`n_e`$ in the right-hand side of equation (1). Thus, under the assumption of cylindrical symmetry, eliminating $`𝐪`$ with the help of equations (2) and (3), the energy equation (1) is reduced to the differential equation

$$\frac{f\kappa _0}{R}\frac{d}{dR}\left(RT^{5/2}\frac{dT}{dR}\right)=P^2\frac{\mathrm{\Lambda }(T)}{T^2},$$
(4)
where $`R`$ is the cylindrical radius.

### 2.2 Boundary conditions

Equation (4) applies outside the edge at $`R_{\mathrm{in}}`$ of the $`T10^4\mathrm{K}`$ cloud: we assume that the cloud contains its own heat sources, such as photoionizing stars, so on its surface the temperature gradient and heat flux vanish. We require $`T`$ to attain some specified value on $`R=R_{\mathrm{out}}`$. Thus we solve equation (4) subject to three boundary conditions, $`T(R_{\mathrm{in}})=T_{\mathrm{in}}=10^4\mathrm{K}`$, $`dT/dR(R_{\mathrm{in}})=0`$, and $`T(R_{\mathrm{out}})=T_{\mathrm{out}}`$, the temperature of the ambient hot plasma, which varies with distance from the cluster centre. For completeness we will explore the wide range $`10^6\begin{array}{c}<\hfill \\ \hfill \end{array}T_{\mathrm{out}}\begin{array}{c}<\hfill \\ \hfill \end{array}10^8\mathrm{K}`$.

Introducing the dimensionless quantities $`\xi =R/R_{\mathrm{out}}`$, $`\tau =T/T_{\mathrm{out}}`$, and $`\stackrel{~}{\mathrm{\Lambda }}_{T_{\mathrm{out}}}(\tau )\mathrm{\Lambda }(\tau T_{\mathrm{out}})/\mathrm{\Lambda }_{23}`$, with $`\mathrm{\Lambda }_{23}=10^{23}\mathrm{erg}\mathrm{s}^1\mathrm{cm}^3`$, equation (4) reads

$$\frac{1}{\xi }\frac{d}{d\xi }\left(\xi \tau ^{5/2}\frac{d\tau }{d\xi }\right)=\lambda \frac{\stackrel{~}{\mathrm{\Lambda }}_{T_{\mathrm{out}}}(\tau )}{\tau ^2},$$
(5)
where

$$\lambda =\frac{P^2R_{\mathrm{out}}^2\mathrm{\Lambda }_{23}}{f\kappa _0T_{\mathrm{out}}^{11/2}}$$
(6)
is an unknown dimensionless parameter that is the eigenvalue of equation (5). Expressed in dimensionless variables, the boundary conditions are $`\tau (\xi _{\mathrm{in}})=\tau _{\mathrm{in}}`$, $`d\tau /d\xi (\xi _{\mathrm{in}})=0`$, and $`\tau (1)=1`$ with $`\xi _{\mathrm{in}}R_{\mathrm{in}}/R_{\mathrm{out}}`$ and $`\tau _{\mathrm{in}}T_{\mathrm{in}}/T_{\mathrm{out}}`$. Due to the presence of $`\stackrel{~}{\mathrm{\Lambda }}_{T_{\mathrm{out}}}`$, the eigenfunctions $`\tau (\xi )`$ and eigenvalues $`\lambda `$ depend on the outer temperature $`T_{\mathrm{out}}`$ in addition to the dimensionless numbers $`\xi _{\mathrm{in}}`$ and $`\tau _{\mathrm{in}}`$.

The inner and outer radii $`R_{\mathrm{in}}`$ and $`R_{\mathrm{out}}`$ cannot be exactly identified with observational quantities. However, we can assume that $`R_{\mathrm{in}}`$ corresponds approximately to the width of the filament, and $`R_{\mathrm{out}}`$ to its length $`l`$. The latter assumption is consistent with the hypothesis of cylindrical geometry, which does not hold at $`R\begin{array}{c}>\hfill \\ \hfill \end{array}l`$. Observed filaments are highly elongated, with ratios length/width $`\begin{array}{c}>\hfill \\ \hfill \end{array}100`$ (Conselice et al. 2001). Thus, we explore a range of values of $`\xi _{\mathrm{in}}`$ that satisfy $`\xi _{\mathrm{in}}1`$.

### 2.3 Numerical solutions

Equation (5) is non–linear and would require numerical integration even if $`\stackrel{~}{\mathrm{\Lambda }}_{T_{\mathrm{out}}}(\tau )`$ were approximated by a simple analytical function. Thus, use of the tabulated cooling function of Sutherland & Dopita (1993) for metallicity $`[\mathrm{Fe}/\mathrm{H}]=0.5`$ introduces no additional computational complexity. We use a shooting method, with Newton iteration, to calculate the solutions.

Numerical integration of equation (5) shows that, for any fixed $`T_{\mathrm{out}}`$ in the range considered, the value of $`\lambda `$ is virtually independent of $`\xi _{\mathrm{in}}`$, provided that $`\xi _{\mathrm{in}}`$ is sufficiently small (i.e., $`\xi _{\mathrm{in}}\begin{array}{c}<\hfill \\ \hfill \end{array}10^2`$). Thus, in what follows we can simplify our treatment by assuming that the eigenvalue $`\lambda `$ depends only on $`T_{\mathrm{out}}`$.

Although the details of the inner temperature profile $`\tau (\xi )`$ depend significantly on $`\xi _{\mathrm{in}}`$, the behaviour of the profile at $`\xi \begin{array}{c}>\hfill \\ \hfill \end{array}10^2`$ is virtually independent of $`\xi _{\mathrm{in}}`$ provided $`\xi _{\mathrm{in}}\begin{array}{c}<\hfill \\ \hfill \end{array}10^3`$. Henceforth we set $`\xi _{\mathrm{in}}\begin{array}{c}<\hfill \\ \hfill \end{array}10^3`$ so that both $`\lambda `$ and the observationally significant part of $`\tau (\xi )`$ are independent on the exact radial range considered.

Fig. 1 shows $`\lambda `$ as a function of $`T_{\mathrm{out}}`$, while Fig. 2 shows a number of the corresponding eigenfunctions $`T(R)`$. For values of $`T_{\mathrm{out}}`$ in the relevant range, $`\lambda `$ spans a relatively small range ($`0.08\begin{array}{c}<\hfill \\ \hfill \end{array}\lambda \begin{array}{c}<\hfill \\ \hfill \end{array}0.37`$), which is reduced to $`0.25\begin{array}{c}<\hfill \\ \hfill \end{array}\lambda \begin{array}{c}<\hfill \\ \hfill \end{array}0.37`$ for $`T_{\mathrm{out}}\begin{array}{c}>\hfill \\ \hfill \end{array}4\times 10^6\mathrm{K}`$.

From Fig. 2 is apparent that over a significant fraction of the radial range the temperature profiles are roughly power laws, with similar indices. This suggests that an approximate analytical solution could be obtained by exploring the asymptotic behaviour of the eigenfunctions of equation (5) for $`\xi 1`$ ($`\tau 1`$). When we assume $`\tau =\xi ^n`$, and approximate the cooling function near $`T_{\mathrm{out}}`$ by $`\stackrel{~}{\mathrm{\Lambda }}_{T_{\mathrm{out}}}(\tau )[\mathrm{\Lambda }(T_{\mathrm{out}})/\mathrm{\Lambda }_{23}]\tau ^\alpha `$, equation (5) requires $`n=4/(112\alpha )`$, and $`\lambda =\frac{7}{2}n^2/[\mathrm{\Lambda }(T_{\mathrm{out}})/\mathrm{\Lambda }_{23}]`$. Since $`1/2\begin{array}{c}<\hfill \\ \hfill \end{array}\alpha \begin{array}{c}<\hfill \\ \hfill \end{array}1/2`$ for the relevant temperatures, the outer slope of the eigenfunction is expected to lie in the small range $`1/3\begin{array}{c}<\hfill \\ \hfill \end{array}n\begin{array}{c}<\hfill \\ \hfill \end{array}2/5`$ (cf. Fig. 2). The run of the eigenvalue as a function of $`T_{\mathrm{out}}`$ (Fig. 1) reflects the proportionality to $`1/\mathrm{\Lambda }(T_{\mathrm{out}})`$, modulated by $`n^2`$, which accounts for the dependence on the local slope of the cooling function $`\alpha `$.

## 3 Applications

### 3.1 Critical length for filaments

We now discuss the physical implications of the results presented in the previous section. Equation (6) relates the numerical value of $`\lambda `$ to the dimensional quantities associated with a filament and its surroundings. If, as explained above, we identify $`R_{\mathrm{out}}`$ with the length of a filament, then this equation yields for a given environment the critical length $`l_{\mathrm{crit}}`$ of filaments, such that smaller filaments evaporate while larger ones grow. Numerically we have

$`l_{\mathrm{crit}}`$ $``$ $`R_{\mathrm{out}}`$ (7)
$``$ $`1.4f^{1/2}\lambda ^{1/2}\left({\displaystyle \frac{T_{\mathrm{out}}}{10^7\mathrm{K}}}\right)^{11/4}\left({\displaystyle \frac{P}{10^6\mathrm{K}\mathrm{cm}^3}}\right)^1\mathrm{kpc}.`$
This equation indicates that $`l_{\mathrm{crit}}`$ has a strong dependence on the outer temperature ($`l_{\mathrm{crit}}T_{\mathrm{out}}^{11/4}`$, neglecting the weak dependence of $`\lambda `$ on $`T_{\mathrm{out}}`$; see Fig. 1), and also is inversely proportional to the pressure $`P`$. Even though the dependence of $`l_{\mathrm{crit}}`$ on the thermal conduction suppression factor $`f`$ is not very strong ($`l_{\mathrm{crit}}f^{1/2}`$), $`f`$ can significantly affect $`l_{\mathrm{crit}}`$ because the value of $`f`$ is very uncertain: $`f1`$ seems to be excluded, but different authors gives values differing by two orders of magnitude, from a few $`10^1`$ to $`10^3`$ – see Narayan & Medvedev (2001) for a discussion. We consider values in the range $`10^3\begin{array}{c}<\hfill \\ \hfill \end{array}f\begin{array}{c}<\hfill \\ \hfill \end{array}1`$ and discuss the possibility of constraining $`f`$ from observations of filaments in Section 3.4.

Consider the implications of the dependence of $`l_{\mathrm{crit}}`$ on $`T_{\mathrm{out}}`$ at fixed pressure. The left panel of Fig. 3 shows $`l_{\mathrm{crit}}(T_{\mathrm{out}})`$ from equation (7), for pressure $`P=10^6\mathrm{K}\mathrm{cm}^3`$, which is a typical value in the core of a cooling-flow galaxy cluster \[e.g., for the Perseus cluster, Conselice et al. (2001), Schmidt, Fabian & Sanders (2002); see also Heckman et al. (1989)\]. A few values of the suppression factor ($`f=1,0.1,0.01,0.001`$) are considered. For given $`f`$, the relation $`l_{\mathrm{crit}}(T_{\mathrm{out}})`$ divides the the diagram in two regions: filaments longer than $`l_{\mathrm{crit}}`$ are condensing, while smaller ones are evaporating. It is apparent from the steep slope of the curves in the diagram that the process of evaporation or condensation of filaments is very sensitive to the temperature of the ambient medium. For a plausible value of the suppression factor, $`f0.01`$, at fixed pressure $`10^6\mathrm{K}\mathrm{cm}^3`$, any filament longer than $`0.1\mathrm{kpc}`$ is condensing in a cluster with ambient temperature $`T_{\mathrm{out}}=10^7\mathrm{K}`$, while filaments as long as $`10\mathrm{kpc}`$ are evaporating if $`T_{\mathrm{out}}=6\times 10^7\mathrm{K}`$.

Although $`l_{\mathrm{crit}}`$ depends on pressure more weakly than on temperature, pressure also plays an important role in determining the fate of filaments. This is apparent in the right diagram of Fig. 3, where we fix $`T_{\mathrm{out}}=3\times 10^7`$ and plot $`l_{\mathrm{crit}}`$ as a function of $`P`$, for the set of values of $`f`$ considered above. The diagram shows that, at fixed temperature, in low-pressure clusters virtually all filaments are evaporating (e.g., $`l_{\mathrm{crit}}100\mathrm{kpc}`$ for $`P\mathrm{few}\times 10^4\mathrm{K}\mathrm{cm}^3`$ and $`f=0.01`$), but condensation is effective in the highest-pressure ambient gas ($`l_{\mathrm{crit}}1\mathrm{kpc}`$ for $`P\begin{array}{c}>\hfill \\ \hfill \end{array}10^6\mathrm{K}\mathrm{cm}^3`$ and $`f=0.01`$).

### 3.2 Time-scales

What are the characteristic time-scales for evaporation and condensation of filaments? A detailed answer to this question would require the solution of the full hydrodynamic equations for non-vanishing mass flow. However, for the purpose of our investigation we are interested in an order of magnitude estimate, which can be obtained by considering that the time variation of temperature is described by a differential equation of the form

$$\frac{dT}{dt}=\frac{T}{t_\kappa }\frac{T}{t_{\mathrm{cool}}},$$
(8)
where $`t_\kappa `$ and $`t_{\mathrm{cool}}`$ are the thermal conduction and cooling time-scales, respectively. The relation between these time-scales is given by

$$\frac{t_\kappa }{t_{\mathrm{cool}}}=\frac{n_e^2\mathrm{\Lambda }(T)}{|𝐪|}\frac{n_e^2\mathrm{\Lambda }(T)l^2}{\kappa (T)T},$$
(9)
because $`|𝐪|\kappa (T)T/l^2`$, where $`l`$ is the temperature variation length scale (Begelman & McKee 1990). In the considered case, $`l`$ is the radial extension of the interface, which is of the order of the filament length (see Section 2). Considering that the equilibrium solution ($`ł=l_{\mathrm{crit}}`$) corresponds to $`t_\kappa t_{\mathrm{cool}}`$, from equation (9) we get

$$t_\kappa \left(\frac{l}{l_{\mathrm{crit}}}\right)^2t_{\mathrm{cool}},$$
(10)
and equation (8) can be approximated as

$$\frac{dT}{dt}=\frac{T}{t_{\mathrm{cool}}}\left[\left(\frac{l_{\mathrm{crit}}}{ł}\right)^21\right].$$
(11)
Consequently, filaments shorter than $`\frac{1}{2}l_{\mathrm{crit}}`$ evaporate in a fraction $`(l/l_{\mathrm{crit}})^2`$ of the cooling time, which for example is $`\begin{array}{c}<\hfill \\ \hfill \end{array}300\mathrm{Myr}`$ within $`4\mathrm{kpc}`$ of the centre of the Hydra cluster, (David et al. 2000). Filaments smaller than $`l_{\mathrm{crit}}`$ by a factor $`\begin{array}{c}>\hfill \\ \hfill \end{array}8`$ have evaporation times $`\begin{array}{c}<\hfill \\ \hfill \end{array}5\mathrm{Myr}`$ that are smaller than their likely dynamical times, $`l/v`$, where $`v100\mathrm{km}\mathrm{s}^1`$ is the typical gas velocity estimated from line widths (Sabra et al. 2000 and references therein). On the other hand, equation (11) in the limit $`ll_{\mathrm{crit}}`$ indicates that, as expected, very long filaments condense on time-scales of the order of the cooling time of the ambient medium.

These considerations strongly support the hypothesis that only filaments with lengths comparable to the critical length or longer are likely to survive for long enough to be observed.

### 3.3 The fate of filaments in cooling-flow and non–cooling-flow clusters

The results presented above have interesting implications, when compared with the properties of observed galaxy clusters.

1. A first simplistic implication is that, due to the strong dependence of $`l_{\mathrm{crit}}`$ on temperature, in massive (high-temperature) clusters most cold gas structures are evaporating, while even relatively small clouds can survive (and condense) in small (low-temperature) clusters. On the other hand one needs to bear in mind that clusters with higher temperatures have larger virial radii $`R_\mathrm{V}`$, and it is natural to expect the typical size of clouds to increase proportional to $`R_\mathrm{V}`$. The left panel of Fig. 3 shows the relation between $`R_\mathrm{V}`$ and the virial temperature ($`R_\mathrm{V}T_\mathrm{V}^{1/2}`$) predicted for galaxy clusters under the assumption of self-similarity (the slope of this scaling relation is consistent with that derived from X-ray observations of galaxy clusters; see, e.g., Ettori, De Grandi & Molendi 2002, and references therein). Thus, as is apparent from the diagram, the relation $`l_{\mathrm{crit}}(T_{\mathrm{out}})`$ is much steeper than the relation between cluster size and temperature. Hence, in the highest-temperature clusters any cold filament is evaporating, while in the lowest-temperature clusters, any resolvable filament is condensing.
2. The scenario presented in the previous point would hold under the assumption that the physical properties (gas density and pressure) of clusters of given temperature would not differ significantly from cluster to cluster. In fact, we know this is not the case: for fixed gas temperature, non–cooling-flow clusters (i.e., ones with long cooling times) have, by definition, lower central density (and pressure) than cooling-flow clusters. Thus, as shown by the right panel of Fig. 3, since $`l_{\mathrm{crit}}1/P`$, a filament long enough to be condensing in the core of a cooling-flow cluster with central temperature $`T`$ might evaporate if it were at the centre of a non–cooling-flow cluster with the same temperature (but lower density). This argument suggests that filaments of cold gas are not observed in non–cooling-flow clusters because in these clusters they are evaporated very efficiently, even at the cluster centre. On the other hand, whatever their formation mechanism, filaments are more likely to survive and grow in the denser, higher pressure cores of cooling-flow clusters.
3. In a given cooling-flow cluster, $`T`$ increases as one moves away from the centre, while $`P`$ falls. Consequently, equation (7) predicts that $`l_{\mathrm{crit}}`$ increases rapidly with distance from the cluster centre. This result explains why, even in the infall or merging scenario for the formation of the filaments, these are observed only in the cooling cores of clusters, and not in the outer regions. For example, in the Perseus cluster, the gas temperature at a distance $`100\mathrm{kpc}`$ from the center is roughly twice the central temperature. Also, in the same radial range the pressure decreases by a factor of $`2.5`$ (e.g., Schimdt et al. 2002; Churazov et al. 2003). Thus, equation (7) predicts that $`l_{\mathrm{crit}}`$ grows by over a factor of 10 between the centre and $`100\mathrm{kpc}`$. Interestingly, even in the core, observations indicate that the mean $`\mathrm{H}\alpha `$ emission is a strongly decreasing function of the cluster radius (e.g., Conselice et al. 2001). This result may arise because as one moves out, $`\mathrm{H}\alpha `$ emission from smaller filaments contributes less and less to the overall signal.

We can summarize the points above by considering the relation between $`l_{\mathrm{crit}}`$ and the cooling time of ambient gas. For simplicity, we use the approximate isobaric cooling time $`t_{\mathrm{cool}}10^4n_e^1T^{1/2}\mathrm{yr}`$ (where $`n_e`$ is in $`\mathrm{cm}^3`$ and $`T`$ in $`\mathrm{K}`$; e.g., Sarazin 1986), considering the temperature range $`10^710^8\mathrm{K}`$. The dotted curves in Fig. 4 represent the loci of constant cooling time in density-temperature space, while the solid curves show the loci of constant $`l_{\mathrm{crit}}`$, for suppression factor $`f=0.01`$. We see that if the cooling time is long (as in the cores of non cooling-flow clusters, and in the outer regions of any cluster), $`l_{\mathrm{crit}}`$ is large, so evaporation is very likely. On the other and, small $`l_{\mathrm{crit}}`$ corresponds to a short cooling time. Hence it is precisely the physical conditions found in the cores of cooling-flow clusters that favour condensation.

### 3.4 Suppression of thermal conductivity in the intracluster medium

As outlined in the Introduction, the existence of cold filaments could be used to constrain the value of the thermal conductivity in the ICM. Indeed, we could derive the value of the suppression factor $`f`$ from equation (7) if we knew all the other quantities involved. Unfortunately, while $`T_{\mathrm{out}}`$ and $`P`$ are well constrained by observations, it is quite difficult to put strong observational constraints on $`l_{\mathrm{crit}}`$. A detailed morphological study of the filaments, available only for the Perseus cluster (Conselice et al. 2001), reveals complex structures on different scales. However, we note that even a quite conservative estimate of $`l_{\mathrm{crit}}`$ can lead to interesting constraints on $`f`$. For instance, the data suggest that in the Perseus cluster filaments $`10\mathrm{kpc}`$ long are present at least up to radii $`4050\mathrm{kpc}`$. This means $`l_{\mathrm{crit}}\begin{array}{c}<\hfill \\ \hfill \end{array}10\mathrm{kpc}`$, where the ICM temperature and pressure are $`4.5\times 10^7\mathrm{K}`$ and $`10^6\mathrm{K}\mathrm{cm}^3`$, respectively (Schmidt et al. 2002). According to equation (7), for these physical properties of the ambient medium, this condition is translated into the upper limit $`f\begin{array}{c}<\hfill \\ \hfill \end{array}0.04`$ for the thermal conduction suppression factor. This result is compatible with that found by Markevitch et al. (2003) on the basis of observations of cold fronts in A754. We note that a substantially higher efficiency would be required in order for the thermal conductivity to significantly affect cooling flows in the centres of clusters (e.g., Narayan & Medvedev 2001).

In the discussion above we considered the thermal conduction suppression factor as a global parameter characterizing the ICM in clusters. In fact, the reduction of thermal conductivity could vary from cluster to cluster, or within a given cluster, depending on the detailed properties of the magnetic field. However, there are indications that the effect of the reconnection of the magnetic field lines between the filaments and the ICM would be to enhance the efficiency of heat conduction in the centres of cooling-flow clusters (Soker, Blanton & Sarazin 2004). Thus, at least according to this scenario, it seems reasonable to extend the derived upper limit on $`f`$ to the ICM in general.

### 3.5 X-ray emission from conduction fronts of filaments

Here we explore the possibility that the conduction fronts described above can be detected in soft X-rays. This is particularly interesting since Fabian et al. (2003) find soft X-ray emission associated with the brightest $`\mathrm{H}\alpha `$ emission filaments around NGC1275, suggesting that it is produced in the conductive evaporation or condensation of the filaments. We quantify the emission properties of the conduction fronts for the considered cylindrically symmetric clouds, computing the differential emission measure

$$EM(T)n_e^2\frac{dV}{d\mathrm{log}T},$$
(12)
where $`n_e`$ is the electron number density and $`dV`$ is the volume element. Under the hypothesis of constant pressure $`P=n_eT`$, the differential emission measure can be related directly to the temperature profile $`T(R)`$: using $`dV=2\pi lRdR`$, with $`l`$ length of the filament, we get

$$EM(T)=2\pi lP^2\frac{R}{T^2}\frac{dR}{d\mathrm{log}T}.$$
(13)
As discussed in Section 2, the temperature profile, for given $`T_{\mathrm{out}}`$, can be considered independent of $`\xi _{\mathrm{in}}`$ at large enough radii ($`\xi \begin{array}{c}>\hfill \\ \hfill \end{array}10^2`$); correspondingly we find that the differential emission measure at $`T\begin{array}{c}>\hfill \\ \hfill \end{array}0.1T_{\mathrm{out}}`$ is not sensitive to the value of $`\xi _{\mathrm{in}}`$. Fig. 5 plots the emission measure as a function of $`T`$ for a subset of values of the outer temperature ($`T_{\mathrm{out}}=1,2,5,10`$ in units of $`10^7\mathrm{K}`$). For the normalization of $`EM(T)`$ we assume $`P=10^6\mathrm{K}\mathrm{cm}^3`$, and $`l=l_{\mathrm{crit}}=R_{\mathrm{out}}`$, in accordance with the discussion in Section 2. Due to the dependence of $`l_{\mathrm{crit}}`$ on the thermal conduction suppression factor $`f`$ (equation 7), the normalization of $`EM(T)`$ is proportional to $`f^{3/2}`$. Thus it is difficult to make predictions on the absolute value of the emission measure, also because of the strong dependence on pressure. On the other hand, the shape of $`EM(T)`$ is not affected by these parameters, being determined by the temperature profile, which depends only on the temperature of the ambient medium. As apparent from Fig. 5, in the considered temperature range the differential emission measure can be roughly approximated with a power-law $`EM(T)T^\beta `$. We find that $`EM(T)`$ is steeper for lower $`T_{\mathrm{out}}`$, with best–fitting slope $`\beta `$ ranging from $`\beta 3.4`$ for $`T_{\mathrm{out}}=10^8\mathrm{K}`$ to $`\beta 4.6`$ for $`T_{\mathrm{out}}=10^7\mathrm{K}`$. Accurate measurements of the soft X-ray spectra around filaments could be interestingly compared with the slopes of these emission measure profiles.

## 4 Conclusions

Embedded within cooling flows, one frequently detects filaments of cold gas. Thermal conduction will cause heat to flow into these filaments from the hot ambient gas, with the result that filaments below a critical size will be evaporated. On the other hand radiative cooling is very efficient in cold, dense gas, so larger filaments can survive, and indeed nucleate the cooling of the ambient medium. We have determined the critical filamentary size that divides the two regimes.

The absolute value of this length scale, $`l_{\mathrm{crit}}`$, is uncertain because we do not know how effectively the magnetic field suppresses thermal conduction. By contrast, the way in which $`l_{\mathrm{crit}}`$ scales with the temperature and pressure of the ambient medium is well determined and leads to suggestive results.

The critical length scale rises sharply with ambient temperature: $`l_{\mathrm{crit}}T^{11/4}`$, and decreases as $`P^1`$ with ambient pressure. Consequently, in the Perseus cluster $`l_{\mathrm{crit}}`$ increases by over an order of magnitude between the cluster centre and a radius of $`100\mathrm{kpc}`$, where it is probably larger than $`10\mathrm{kpc}`$. Within a few kiloparsecs of the cluster centre, filaments only a kiloparsec long may either grow, or at least take a long time to evaporate because their lengths are just below the critical value, while such small filaments are rapidly evaporated in the main body of the cluster.

Where do filaments come from in the first place? They cannot emerge spontaneously from the ambient medium because this is thermally stable,222Balbus (1991) shows that thermal instability is possible in the presence of even a weak magnetic field. However, a weak field can only control low-amplitude perturbations, and buoyant stabilization will reassert itself once the perturbation amplitude is greater than $`P_{\mathrm{mag}}/P_{\mathrm{gas}}10^3`$ in clusters. It is also the case that for technical reasons Balbus had to insist on a highly artificial background field for which field-line closure is problematic. so they must be injected. Simulations of structure formation predict that significant amounts of cold gas fall into clusters (Katz et al. 2003), and the observed filaments could represent just the fraction of this infall that reaches the centre. There is also some indication that extremely close to the cluster centre a small amount of gas cools to filamentary temperatures, and some of this gas could be ejected by the AGN to the radii where filaments are observed.

Whether infall or AGN are responsible for injecting filamentary gas, there is no reason to expect the lengths of the injected filaments to have a characteristic value. It is natural to think in terms of a spectrum of filamentary lengths $`l`$, in which the number of filaments decreases with increasing filament length, probably as a power law. Observations of filaments would then be dominated by filaments with $`ll_{\mathrm{crit}}`$ because shorter filaments would be rapidly evaporated, and larger ones would be rare. The local number density of filaments would be a decreasing function of $`l_{\mathrm{crit}}`$, explaining why observed filaments are concentrated near the centres of cooling-flow clusters.

Since filament growth depends on the efficiency of radiative cooling, $`l_{\mathrm{crit}}`$ tends to be small when the cooling time of the ambient medium is short. Hence clusters with short central cooling times, have smaller values of $`l_{\mathrm{crit}}`$ and more filaments than clusters with similar temperatures but longer cooling times. Consequently, it is natural that filaments are only detected in cooling-flow clusters.

For the case of a filament of critical length, we have calculated the temperature profile throughout the interface between gas at $`10^4\mathrm{K}`$ that emits optical emission lines, and the ambient X-ray emitting plasma. This profile yields the run of differential emission measure with temperature, which can be tested observationally.

As cosmic clustering proceeds, the temperature of gas that is trapped in gravitational potential wells rises as longer and longer length-scale perturbations collapse. This temperature rise causes $`l_{\mathrm{crit}}`$ to rise much faster than the virial radius of the potential well. Hence, while at early times filaments that are much smaller than the system can grow, growth is later restricted to filaments that are comparable in size to the entire system. Consequently, the ability of infalling cold gas to nucleate star formation declines as clustering proceeds.

Semi-analytic models of galaxy formation that can account for the observed slope of the galaxy luminosity function at the faint end need efficient feedback from star formation, and over-produce luminous galaxies (Benson et al. 2003). Binney (2003) has argued that galaxy formation is dominated by infall of cold gas, and that it is the efficient evaporation of filaments in massive systems that suppresses the formation of massive galaxies. Our results make a step towards quantifying this picture.

## Acknowledgments

We thank Luca Ciotti for helpful comments on the manuscript, and the referee, Noam Soker, for interesting suggestions that helped improve the paper. C.N. was partially supported by a grant from Università di Bologna.

## Appendix A Thermal conduction saturation parameter

We have hitherto assumed that the thermal conduction is not saturated, so the heat flux is given by the classical equation (2), with the thermal conductivity given by equation (3). The latter formula is derived under the assumption that $`l_e/l_T1`$, where $`l_e`$ is the electron mean free path and $`l_T=T/|T|`$ is the temperature scale height. When this condition is violated, equation (2) overestimates the heat flux, which is limited by the thermal speed of electrons, and the modulus of $`𝐪`$ is given by the saturated heat flux (Cowie & McKee 1977)

$$q_{\mathrm{sat}}=0.4\left(\frac{2k_\mathrm{B}T}{\pi m_e}\right)^{1/2}n_ek_\mathrm{B}T,$$
(14)
where $`k_\mathrm{B}`$ is the Boltzmann constant, and $`m_e`$ is the electron mass. Cowie & McKee introduced a discriminant of whether the heat flux in a plasma is saturated: this is the saturation parameter $`\sigma q_{\mathrm{cl}}/q_{\mathrm{sat}}`$, where $`q_{\mathrm{cl}}=\kappa (T)|T|`$ is the modulus of the classical heat flux (equation 2). The critical value that separates the two regimes is $`\sigma 1`$: for $`\sigma 1`$ the heat flux is unsaturated, while for $`\sigma 1`$, equation (14) applies.

For given $`T_{\mathrm{out}}`$ and corresponding critical temperature profile $`T(R)`$, the value of the saturation parameter as a function radius in the interface can be obtained combining equations (2), (3), and (14), and using the relations $`n_e=P/T`$ and $`R_{\mathrm{out}}=l_{\mathrm{crit}}`$ (equation 7). Numerically we get

$$\sigma (R)8\times 10^3f^{1/2}\lambda ^{1/2}\left(\frac{T_{\mathrm{out}}}{10^7\mathrm{K}}\right)^{1/4}\tau ^2\frac{d\tau }{d\xi },$$
(15)
where $`\xi R/R_{\mathrm{out}}`$, $`\tau =\tau (\xi )`$ is the normalized temperature profile, and $`\lambda =\lambda (T_{\mathrm{out}})`$ the corresponding eigenvalue (see Section 2). We note that $`\sigma (R)`$ is independent of the pressure and scales as the square root of the suppression factor $`f`$. Applying equation (15) to conduction fronts with ambient gas temperature in the considered range $`10^6\begin{array}{c}<\hfill \\ \hfill \end{array}T_{\mathrm{out}}\begin{array}{c}<\hfill \\ \hfill \end{array}10^8\mathrm{K}`$, we find in all cases $`\sigma (R)\begin{array}{c}<\hfill \\ \hfill \end{array}10^2f^{1/2}`$. Thus, the condition for unsaturated heat conduction ($`\sigma 1`$) is satisfied throughout the interface, even for the largest values of the suppression factor $`f1`$.


format: text

text:
# Tailoring ferromagnetic chalcopyrites

If magnetic semiconductors are ever to find wide application in real spintronic devices, their magnetic and electronic properties will require tailoring in much the same way that band gaps are engineered in conventional semiconductors. Unfortunately, no systematic understanding yet exists of how, or even whether, properties such as Curie temperatures and band gaps are related in magnetic semiconductors. Here we explore theoretically these and other relationships within 64 members of a single materials class, the Mn-doped II-IV-V2 chalcopyrites, three of which are already known experimentally to be ferromagnetic semiconductors. Our first-principles results reveal a variation of magnetic properties across different materials that cannot be explained by either of the two dominant models of ferromagnetism in semiconductors. Based on our results for structural, electronic, and magnetic properties, we identify a small number of new stable chalcopyrites with excellent prospects for ferromagnetism.

One of the tantalizing promises of dilute magnetic semiconductors is to conjoin, in a single material, the advantages of nonvolatility and bandgap engineeringohno\_science\_1998a ; maekawa2002a ; zutic\_rev\_mod\_phys\_2004a . In normal (nonmagnetic) semiconductors, precise tailoring of band structure is possible due to our detailed understanding—both empirical and theoretical—of the mechanisms underlying band formation. Unlike the situation for normal semiconductors, however, there does not yet exist any comprehensive understanding of how the magnetic properties of magnetic semiconductors are related to their structural and electronic properties.

For example, there is no firm understanding of the relationship between the properties of the host semiconductor and the Curie temperatures attainable by doping with magnetic impurities such as Mn. Experimentally, a number of materials issues—including Mn content, compensation, and phase purity—are not yet sufficiently under control to permit a systematic description across different semiconductors. Theoretically, there is an emerging consensus that while the ferromagnetism originates from the interaction between itinerant electrons (or holes) and localized electrons, there is no single model that can systematically relate the physical and electronic properties of the host semiconductor to its resulting magnetic properties when dopeddietl\_nature\_mater\_2003a . Models based on the weak- and strong-coupling limits of this interaction (the Zenerdietl\_science\_2000a ; konig\_phys\_rev\_lett\_2000a  and double-exchangechattopadhyay\_phys\_rev\_lett\_2001a ; nagaev2002a  models, respectively) do make specific predictions, but it is increasingly evident that real materials are not sufficiently well described by either limit for these predictions to be reliable across different hosts.

Here we address this problem by exploring theoretically these relationships within a single materials class: the II-IV-V2 chalcopyrites consisting of the 64 possible combinations of II=(Be,Mg,Zn,Cd), IV=(C,Si,Ge,Sn), and V=(N,P,As,Sb). Three of these—CdGeP2, ZnGeP2, and ZnSnAs2—were recently shown to become ferromagnetic upon Mn doping, with remarkably high Curie temperaturesmedvedkin\_japanese\_journal\_of\_applied\_physics\_2000a ; medvedkin\_journal\_of\_crystal\_growth\_2002a ; cho\_phys\_rev\_lett\_2002a ; choi\_solid\_state\_comm\_2002a . It is an open question whether other chalcopyrites can be made ferromagnetic; nor is it known which properties of the host chalcopyrite are important for ferromagnetism and, in particular, for high Curie temperatures. These are the central issues we address in this work.

For each chalcopyrite, we use density-functional theory (DFT) in the generalized-gradient approximationperdew\_phys\_rev\_b\_1992a ; kresse\_phys\_rev\_b\_1993a ; kresse\_phys\_rev\_b\_1996a  to compute several properties of the undoped hosts: equilibrium crystal structure, i.e. lattice parameters $`a`$ and $`c`$ and the single internal coordinate; phase stability, as determined by the enthalpy of formation; and electronic band gap, which may be direct or indirect. A comprehensive survey of the magnetic properties within such a large class of materials is impractical, so we consider here only a few key quantities that control magnetism. Since Mn is isoelectronic with the group-II elements, it is believed that carrier-mediated ferromagnetism in the II-IV-V2 chalcopyrites is favored when Mn occupies the group-IV site, where it is likely to be electrically activemahadevan\_phys\_rev\_lett\_2002a . For this reason we first explore the competition between substitutional Mn doping on the II and IV sites by computing the Mn impurity formation energy on each site. Second, we address whether occupation of the group-IV sites is indeed likely to lead to a ferromagnetic phase; in particular, for each host we compute the Heisenberg spin coupling for an experimentally reasonable Mn concentration. The magnitude (and sign) of this coupling should serve as a useful roadmap for the experimental exploration of new ferromagnetic chalcopyrites.

Of the 64 chalcopyrites considered here, about 1/3 are known to exist and have been studied experimentally. A few crystallize in both the chalcopyrite and sphalerite structure (e.g. CdGeP2 and ZnSnSb2), while two take the wurzite structure (BeSiN2 and ZnGeN2); for simplicity we treat all 64 in the chalcopyrite structure. For those with measured lattice parameters, we can compare our theoretical predictions. As expected from DFT calculations, $`a`$ is given very accurately (within $``$1%), while $`c`$ is less so ($``$15%). As Fig. 1 shows, nearly half of the 64 chalcopyrites enjoy a close lattice match to an important mainstream semiconductor.

As with all semiconductors, DFT in the generalized-gradient approximation significantly underestimates band gaps. For chalcopyrites whose measured gaps are in the range 1–2 eV, the DFT values are approximately 1 eV smaller than measured, consistent with previous findingscontinenza\_phys\_rev\_b\_1992a . Hence the theoretical gaps shown in Fig. 1 must be viewed with caution: the true gaps are likely to be at least 1 eV larger than predicted. Likewise, hosts that have no gap within DFT may in reality have a gap (which is likely to be 1 eV or smaller). These caveats aside, the trend in Fig. 1 is standard for covalently bonded semiconductors: smaller lattice constants are associated with larger band gaps. In particular, the nitrides have the smallest lattice constants and largest gaps, while the antimonides have the largest lattice constants and smallest gaps. From the perspective of our survey, the range and distribution of lattice constants and band gaps is satisfyingly large, guaranteeing our roadmap good coverage of this unexplored territory.

Not all of the chalcopyrites considered here are thermodynamically stable, even in the absence of Mn. Fig. 2 shows the calculated enthalpies of formation, $`\mathrm{\Delta }H_f(\text{II-IV-V}_2)=E(\text{II-IV-V}_2)E(\text{II})E(\text{IV})2E(\text{V})`$, where the latter terms refer to the ground-state elemental phases. Of the 64 chalcopyrites, 3/4 have negative enthalpies, suggesting many potentially stable host materials. (We have not considered phase stability with respect to disproportionation into compound products, and so the terms “stable” and “unstable” should be understood here with this caveat.) The trend in Fig. 2 is similar, and indeed intimately related, to that of Fig. 1: smaller lattice constants are usually associated with higher stability. Hence, as expected, the nitrides tend to be the most stable chalcopyrites and the antimonides the least stable. Note that there are important deviations from this trend: for example, nearly all the carbides are unstable.

The solubility of Mn in a host chalcopyrite is an extremely important issue for two reasons. The first pertains to the well-known tendency, within Mn-doped DMS materials, toward phase separation with increasing Mn contentohno\_science\_1998a ; maekawa2002a . This tendency is reflected in the very low equilibrium solubility limit of Mn in, for example, III-V hosts—typically of order parts per thousand or less. These limits can be greatly exceeded, by 2–3 orders of magnitude, by relying on kinetic barriers to maintain the metastable phase. It is not feasible to address this issue theoretically for the chalcopyrites. Instead, we propose using the theoretical solubility limit of Mn as a rough proxy for the metastability doping limit. The solubility limit of Mn, $`x_{\mathrm{sol}}`$, is determined by its impurity formation energy, $`\mathrm{\Delta }H_f(\text{Mn})`$, simply according to $`x_{\mathrm{sol}}=\mathrm{exp}(\mathrm{\Delta }H_f(\text{Mn})/kT)`$. In this formulation, a negative formation energy implies that spontaneous incorporation of Mn impurities will be limited only by kineticszhang\_j\_phys\_cond\_mat\_2002a .

The second reason concerns ferromagnetism per se. In general, the solubility limits for Mn substituting on the group-II and -IV sites will be different. Mahadevan and Zunger have shown that, in the case of CdGeP2, substitutional MnII leads to antiferromagnetic interactions while MnIV leads to ferromagnetic interactionsmahadevan\_phys\_rev\_lett\_2002a . Hence the II-IV solubility difference will greatly influence the stability, or even existence, of the ferromagnetic phase, as well as the attainable Curie temperatures. There are two other important issues, which we do not address here. One is the possibility of Mn occupying interstitial sites (which is increasingly likely for hosts with small lattice constantserwin\_phys\_rev\_b\_2003a ) or forming complexes with other Mn or native defects; if electrically active, these can have consequences for ferromagnetism. The other is the possibility that even substitutional Mn may be subject to a Jahn-Teller distortion, which can affect the magnetic interactions as discussed in Ref. \[kacman\_semicond\_sci\_technol\_2001a, \].

It is important to realize that the impurity formation energy (and hence the solubility limit) in a multicomponent host is not uniquely defined, but depends on the growth environment—in particular, on the energy $`\mu `$ required to transfer one atom from a reservoir into the host compoundvandewalle\_phys\_rev\_b\_2001a . For the ternary chalcopyrites there are three such chemical potentials, {$`\mu _{\mathrm{II}}`$, $`\mu _{\mathrm{IV}}`$, $`\mu _\mathrm{V}`$}, which are related by $`\mu _{\mathrm{II}}+\mu _{\mathrm{IV}}+2\mu _\mathrm{V}=\mathrm{\Delta }H_f(\text{II-IV-V}_2).`$ The formation energies shown in Fig. 3 were evaluated for $`\mu _{\mathrm{II}}`$ and $`\mu _{\mathrm{IV}}`$ at their highest allowed values, corresponding to a simultaneously II-rich and IV-rich growth environment. This condition can often be realized experimentally, unless a more stable compound intervenes. With this choice, half of the chalcopyrite hosts (including all of the carbides) prefer Mn on the group-II site, while half (including most of the silicides) prefer the group-IV site. For both the II and IV sites, there is the expected trend of larger lattice constants leading to smaller impurity formation energies—and hence to higher solubilities. This trend begins to disappear for lattice constants larger than $``$6 Å; in this regime there is enough room available for the Mn impurity that the formation energies become essentially independent of both site and host.

In the remainder of this work we assume that the prospects for ferromagnetism are best served when Mn substitutes on the group-IV sitemahadevan\_phys\_rev\_lett\_2002a . The doping can be biased in this direction by changing the growth environment to the simultaneously II-rich and IV-poor condition, which reduces the relative MnIV formation energy, $`\mathrm{\Delta }E=\mathrm{\Delta }H_f(\text{Mn}\text{IV})\mathrm{\Delta }H_f(\text{Mn}\text{II})`$, to its minimum physically allowed value. The magnitude of the reduction depends on the nature of other possible intervening phases IIxIVyVz, which we do not pursue here (see Ref. \[zhao\_phys\_rev\_b\_2004a, \] for a related discussion). In the best-case scenario, the limiting IV-poor condition will lead to a change in $`\mathrm{\Delta }E`$ by an amount $`\mathrm{\Delta }H_f(\text{II-IV-V}_2)`$, relative to the values shown in Fig. 3. Using our calculated enthalpies, we find this limiting IV-poor condition favors Mn on the group-IV site in 50 of the 64 hosts (the carbides are again the dominant exception). Thus we conclude that under suitable growth conditions, substitutional Mn can be forced to preferentially occupy the group-IV site in the majority of chalcopyrite hosts.

Even in such favorable cases, one expects some fractional occupation of Mn on group-II sites. Since this doping is isoelectronic it is unlikely to produce either holes or electrons, and thus will probably require additional doping to order ferromagnetically. A more pressing question is whether Mn doping on the group-IV sites necessarily leads to ferromagnetism in every II-IV-V2 chalcopyrite. And, for those hosts that can be made ferromagnetic, it is of great interest to know what Curie temperatures one can expect.

These are very difficult questions to address comprehensively, in a materials-specific framework, for any random magnetic alloy. Instead, we adopt again our earlier view and focus on a simple proxy for ferromagnetism, namely the nearest-neighbor Heisenberg coupling, $`J`$, for a specific distribution of Mn in each host. In the mean-field theory of Heisenberg ferromagnetism, the Curie temperature is simply proportional to $`|J|`$, with a proportionality constant depending on the spatial distribution of Mn. Nevertheless, we caution that while qualitative trends in $`J`$ are likely to be meaningful within this materials class, the quantitative prediction of Curie temperatures in ferromagnetic semiconductors remains a difficult problem. Here we choose a fixed, experimentally reasonable Mn concentration of 12.5% (see, for example, Ref. \[cho\_phys\_rev\_lett\_2002a, \]), and examine how $`J`$ varies across different hosts. In a related study of the I-III-VI2 chalcopyrite CuGaSe2, Picozzi et al. showed that $`|J|`$ increases with Mn concentrationpicozzi\_phys\_rev\_b\_2002a . Our method for computing $`J`$ is the same as in Ref. \[picozzi\_phys\_rev\_b\_2002a, \], and gives nearly identical results for CuGaSe2 with 12.5% Mn.

Surprisingly, we find that ferromagnetic alignment is favored ($`J`$ is negative) in only half of the 64 chalcopyrites. Three main factors determine the behavior of $`J`$ across different hosts. The first is the group-II element: the majority of Mg- and Cd-containing hosts favor ferromagnetic alignment, while the majority of Be- and Zn-containing hosts favor antiferromagnetic alignment. The second factor is the group-IV element: for many combinations of II and V elements, we find the ordering $`J(\mathrm{C})<J(\mathrm{Si})<J(\mathrm{Ge})<J(\mathrm{Sn})`$, except in the nitrides (see Supplementary Information, Fig. S1). The group-V element plays the third—and lesser—role, with the ordering $`J(\mathrm{P})<J(\mathrm{As})<J(\mathrm{Sb})`$ obtained for most cases except the carbides (see Supplementary Information, Fig. S2). None of these trends can be attributed to the approximate scaling, $`Ja^3`$, predicted by the mean-field solution of the Zener modeldietl\_phys\_rev\_b\_2001a ; indeed, across this class of chalcopyrites we find no universal dependence of $`J`$ on lattice constant (see Supplementary Information, Fig. S3). Nor do we find any simple correlation between $`J`$ and the bandgap; such a correlation might have been expected if the small-gap hosts reflected Zener physics (with modest Curie temperatures) and the large-gap hosts double-exchange physics (with higher Curie temperatures)pearton\_j\_appl\_phys\_2003a .

Fig. 4 shows the computed values of $`J`$ plotted against the relative MnIV impurity formation energy, $`\mathrm{\Delta }E`$. We do not imply any correlation between these two quantities, but rather display them together because each is centrally important for designing new ferromagnetic chalcopyrites. In particular, the lower-left quadrant contains chalcopyrites for which occupation by, and ferromagnetic alignment of, Mn on the group-IV site is favorable under experimentally plausible II- and IV-rich growth conditions. This conservatively defined sweet spot contains just 19 chalcopyrites out of the full set of 64. Of these, only eight are thermodynamically stable: four are lattice-matched to a mainstream semiconductor (BeSnN2, MgGeP2, MgSiAs2, MgGeAs2) and four are not (BeGeN2, MgSiSb2, MgGeN2, MgGeSb2). Seven of these eight have values of $`J`$ comparable to, or larger than, that of CdGeP2 and thus can be expected to have similar Curie temperatures—of order room temperature or larger. This small set of materials should provide a propitious starting point for more detailed experimental and theoretical scrutiny.

One can also consider the best-case scenario described earlier and consider $`\mathrm{\Delta }E`$ in the limiting IV-poor growth condition. This leads to horizontal shifts, by $`\mathrm{\Delta }H_f`$, of the points in Fig. 4—leftward for negative enthalpies and rightward for positive—and thereby moves 13 more chalcopyrites into the lower left sweet spot. Ten of these are lattice-matched to mainstream semiconductors (MgSiN2, MgSiP2, MgSnP2, MgSnAs2, ZnSiP2, ZnGeN2, CdSiP2, CdSiAs2, CdGeAs2, CdSnP2), and three are not (BeSiP2, CdSiN2, CdGeP2); all are thermodynamically stable.

Experimentally, CdGeP2 and ZnSnAs2 have ferromagnetic ground states with comparable Curie temperatures, while ZnGeP2 has an antiferromagnetic ground state (ferromagnetism disappears below 47 K \[cho\_phys\_rev\_lett\_2002a, \]). Our results, which presume Mn doping on the group-IV site, predict CdGeP2 to be ferromagnetic but ZnSnAs2 and ZnGeP2 to be antiferromagnetic. Reconciling these findings is problematic, because the experimental group-II and -IV chemical potentials—and thereby the favored Mn site—are in general not known. The available evidence points to MnIV substitution in CdGeP2 \[mahadevan\_phys\_rev\_lett\_2002a, \], MnII in ZnGeP2 \[cho\_phys\_rev\_lett\_2002a, \], and mixed occupation in ZnSnAs2 \[choi\_solid\_state\_comm\_2002a, \]. Future work addressing the magnetic interactions between Mn on group-II sites is in progress.

Finally, we address the prospects for using ferromagnetic chalcopyrites as sources of spin-polarized current for injection into nonmagnetic semiconductors. Many issues play a role in this process, including the spin polarization of states at the Fermi level, the band velocities of those states, native defects and other structural inhomogeneities, and spin scattering at the interfacezutic\_rev\_mod\_phys\_2004a . Here we focus on the first of these. For dilute Mn doping (12.5% and 6.25%, and presumably below this as well), all of the 21 stable ferromagnetic chalcopyrites identified above are “idealized half-metals,” by which we mean at that at zero temperature and for an ordered arrangement of Mn, there are states at the Fermi level only in one spin channel. The effects of both finite temperature and site disorder will generally reduce the spin polarization from this idealized case. At higher Mn concentration (50%), half-metallicity persists for all but two of these 21 hosts (BeSnN2, CdSiN2). Even when Mn completely occupies the group-IV sublattice, the resulting stoichiometric compounds are half-metallic in all cases except for the nitrides. Thus the prospects for using ferromagnetic chalcopyrites as spin-polarized sources appear quite favorable.

In summary, we have examined theoretically the prospects for ferromagnetism within the class of all possible II-IV-V2 chalcopyrites having constituents chosen from the first four rows of the Periodic Table. As expected for semiconductors, we find that the electronic properties of the host materials (band gap and enthalpy of formation) are closely related to their structural properties (lattice constant, $`a`$). Contrary to conventional wisdom, we do not find that Mn is generally more soluble on the group-II site than on the group-IV site; instead, the relative site solubility depends on both the nature of the host (e.g. carbide vs. silicide) and the choice of growth condition, while the absolute Mn solubility is largely determined by the host lattice constant. Our results for the Heisenberg coupling between Mn show that the Curie temperatures expected in Mn-doped chalcopyrites do not follow the scaling with lattice constant predicted by the Zener model; nor do they show any systematic variation with band gap, as the double-exchange model would predict. By identifying those chalcopyrites that simultaneously exhibit thermodynamic stability, favorable Mn doping site, and ferromagnetic Mn interactions, we identify two small sets of chalcopyrites that show excellent prospects for stable ferromagnetism under realistic and attainable experimental conditions.

This work was supported by ONR and the DARPA SpinS program. I.Ž. acknowledges financial support from the National Research Council. Discussions with G.A. Medvedkin are gratefully acknowledged. Computations were performed at the DoD Major Shared Resource Center at ASC.

Methods. All the numerical results reported here are based on density-functional theory in the generalized-gradient approximationperdew\_phys\_rev\_b\_1992a , using ultrasoft pseudopotentials as implemented in vasp kresse\_phys\_rev\_b\_1993a ; kresse\_phys\_rev\_b\_1996a . For each host chalcopyrite, the plane-wave cutoff was separately determined by the three constituent elements (four if Mn was included). For the undoped hosts, equilibrium lattice parameters and the internal coordinate were optimized using a 4$`\times `$4$`\times `$4 Monkhorst-Pack sampling of the Brillouin zone. Band gaps were evaluated by identifying band edges based on a 12$`\times `$12$`\times `$12 sampling that included the zone center and high-symmetry points.

For enthalpies of formation, the total energies of the elemental phases were calculated using the following structures, which are either ground-state phases or energetically very close thereto: hcp (Be, Mg, Cd, Zn); diamond (C, Si, Ge); white-tin (Sn); molecular dimer (N); black phosphorous (P); $`\alpha `$As (As, Sb).

Mn impurity formation energies were computed using simple cubic (or, for non-ideal $`c/a`$, tetragonal) supercells of 64 host atoms. To properly represent the dilute impurity limit, all atomic positions were relaxed with fixed ideal supercell lattice parameters, within the constraint that Mn remain on-center. A single Monkhorst-Pack $`k`$-point was used. Chemical potentials were defined with respect to their thermodynamic upper limit determined by the elemental phases listed above. For the Mn chemical potential, a non-magnetic fcc structure was used to approximate the more complicated $`\alpha `$Mn ground-state phase.

Heisenberg spin couplings for 12.5% Mn-doped chalcopyrites were computed using the method described in Ref. picozzi\_phys\_rev\_b\_2002a, . The lattice parameters and internal coordinate of each supercell were assumed to depend linearly on the Mn content (Vegard’s law), and were interpolated using lattice parameters and internal coordinate calculated for the fully MnIV-substituted host. Two Mn per 64-atom supercell were then arranged on an uniform lattice and full atomic relaxation was performed. Total energies for two spin configurations were computed—ferromagnetic, and antiferromagnetic with $`[`$100$`]`$ ordering wavevector—whereupon the energy difference gives $`J`$ within the nearest-neighbor Heisenberg model.

Supplementary Information


text:
# Bianchi IX Chaoticity: BKL Map and Continuous Flow

## Abstract

We analyze the Bianchi IX dynamics (Mixmaster) in view of its stochastic properties; in the present paper we address either the original approach due to Belinski, Khalatnikov and Lifshitz (BKL) as well as a Hamiltonian one relying on the Arnowitt–Deser–Misner (ADM) reduction.
We compare these two frameworks and show how the BKL map is related to the geodesic flow associated with the ADM dynamics. In particular, the link existing between the anisotropy parameters and the Kasner indices is outlined.

The Bianchi classification characterizes all the admissible three-dimensional homogeneous spaces and determines 9 non-equivalent types. The Bianchi models were widely studied in cosmology as the simplest generalization of the Friedmann–Lemaitre–Robertson–Walker Universe, the Einstein equations reducing to an ordinary differential system with respect to time .
Among the Bianchi classification, the types VIII and IX (Mixmaster) were particularly attractive from a cosmological point of view, being the most general ones allowed within the homogeneity constraint and characterized by a chaotic dynamics near the Big Bang; here we discuss the type IX in some detail, yet all considerations hold also for the type VIII.

The Bianchi IX model is summarized by the line element

$$ds^2=N^2(t)dt^2\gamma _{ij}\sigma ^i\sigma ^j,i,j=1,2,3$$
(1)
where $`N`$ denotes the lapse function and the $`\sigma ^i`$ are the 1-forms defined by the homogeneity property. In vacuum we can take

$$\gamma _{ij}(t)=\text{diag}\left(e^{h_i(t)}\right)$$
(2)
and the Einstein equations provide

$`h_i^{\prime \prime }=\left(e^{h_j}e^{h_k}\right)^2e^{2h_i},ijk`$ (3a)
$`{\displaystyle \underset{i}{}}h_i^{\prime \prime }={\displaystyle \frac{1}{2}}{\displaystyle \underset{jk}{}}h_j^{}h_k^{}`$ (3b)
where $`()^{}=d/d\eta `$ and the time $`\eta `$ corresponds to the choice $`N=\text{exp}\left(_i\frac{h_i}{2}\right)`$; equation (3b) plays the role of a first integral for the system (3a).
Belinski, Khalatnikov and Lifshitz (BKL)  showed that the asymptotic approach to the Big Bang of the type IX model is characterized by an infinite sequence of Kasner epochs, i.e. intervals of time during which the right-hand side of (3a) is negligible. Thus, in the limit $`\eta \mathrm{}`$, we have a piecewise solution, whose single steps take the Kasner form

$$h_ip_i\eta +h_i^{},(p_i,h_i^{})=\text{const.},$$
(4)
where the Kasner indices $`p_i`$ satisfy the conditions $`_ip_i=_i(p_i)^2=1`$.
The Mixmaster dynamics is then reduced to a map on these indices: adopting the parametrization which orders the $`p_i`$’s as

$$p_1=\frac{u}{1+u+u^2},p_2=\frac{1+u}{1+u+u^2},p_3=\frac{u(1+u)}{1+u+u^2},u(1,\mathrm{})$$
(5)
then we get the following discrete map on the parameter $`u`$

$$u^{}=\{\begin{array}{c}u1u>2\hfill \\ \frac{1}{u1}1<u2\hfill \end{array}$$
(6)
The main achievement of the BKL approach is the discovery that such a map has stochastic properties and $`u`$ admits the probability distribution $`w`$

$$w(u)=\frac{1}{u(u+1)\mathrm{ln2}}.$$
(7)
In  the above analysis was extended, point by point in space, to a generic inhomogeneous cosmological model.
To get a continuous description, let us introduce the Misner variables $`(\alpha ,\beta _+,\beta _{})`$ and then the Misner–Chitre-like ones $`(\tau ,\xi ,\theta )`$, via the transformations

$`h_1=\beta _++\sqrt{3}\beta _{},h_2=\beta _+\sqrt{3}\beta _{},h_3=2\beta _+`$ (8a)
$`\alpha =e^\tau \xi ,\beta _\pm =e^\tau \sqrt{\xi ^21}\times \{\begin{array}{c}\mathrm{cos}\theta \hfill \\ \mathrm{sin}\theta \hfill \end{array}`$ (8d)
with $`\xi (1,\mathrm{})`$ and $`\theta [0,2\pi )`$.
By an Arnowitt–Deser–Misner (ADM) reduction of the Einstein–Hilbert action, it can be shown  that the Bianchi IX dynamics is summarized by

$$S_{\text{IX}}=𝑑\tau \left(p_\xi \frac{d\xi }{d\tau }+p_\theta \frac{d\theta }{d\tau }\sqrt{\epsilon ^2+U}\right),$$
(9)
where $`\epsilon ^2=\left(\xi ^21\right)p_\xi ^2+\frac{p_\theta ^2}{\xi ^21}`$.
Near the Big Bang $`(\tau \mathrm{})`$, the potential is modelled by the infinite walls

$$U=\underset{i}{}\mathrm{\Theta }_{\mathrm{}}(H_i),\mathrm{\Theta }_{\mathrm{}}(x)=\{\begin{array}{c}0x>0\hfill \\ \mathrm{}x0\hfill \end{array}$$
(10)
where the quantities $`H_ih_i/(_kh_k)`$, $`(_iH_i=1)`$ are the anisotropy parameters which in the Misner–Chitre-like variables read as

$`H_1={\displaystyle \frac{1}{3}}{\displaystyle \frac{\sqrt{\xi ^21}}{3\xi }}\left(\mathrm{cos}\theta +\sqrt{3}\mathrm{sin}\theta \right)`$ (11a)
$`H_2={\displaystyle \frac{1}{3}}{\displaystyle \frac{\sqrt{\xi ^21}}{3\xi }}\left(\mathrm{cos}\theta \sqrt{3}\mathrm{sin}\theta \right)`$ (11b)
$`H_3={\displaystyle \frac{1}{3}}+2{\displaystyle \frac{\sqrt{\xi ^21}}{3\xi }}\mathrm{cos}\theta ;`$ (11c)
the functions $`H_i`$ expressed in (11) do not depend on the (time) variable $`\tau `$.
The domain $`\mathrm{\Gamma }_H`$ where $`U`$ vanishes is dynamically closed and, within it, $`\epsilon `$ behaves as a constant of motion, i.e. $`\frac{d\epsilon }{d\tau }=\frac{\epsilon }{\tau }=0`$, which implies $`\epsilon =E=\text{const.}`$.
Thus the Bianchi IX dynamics is isomorphic to a billiard on a Lobachevski plane whose line element reduces to

$$dl^2=E^2\left[\frac{d\xi ^2}{\xi ^21}+(\xi ^21)d\theta ^2\right].$$
(12)
This system is stochastic  and, upon the transformation of the conjugate momenta as

$$p_\xi =\epsilon \frac{\mathrm{cos}\varphi }{\sqrt{\xi ^21}},p_\theta =\epsilon \sqrt{\xi ^21}\mathrm{sin}\varphi ,0\varphi <2\pi ,$$
(13)
the dynamics admits a uniform invariant measure like

$$d\mu =\frac{1}{8\pi ^2}d\xi d\theta d\varphi .$$
(14)
The link between this approach based on continuous variables and the one based on the discrete BKL one can be achieved considering the following set of coordinates transformations

$`\stackrel{}{y}`$ $`={\displaystyle \frac{1+\xi }{\sqrt{\xi ^21}}}(\mathrm{cos}\theta ,\mathrm{sin}\theta ),`$ (15a)
$`\stackrel{}{\iota }`$ $`={\displaystyle \frac{\stackrel{}{y}+\stackrel{}{b}}{(\stackrel{}{y}+\stackrel{}{b})^2}}\stackrel{}{b},`$ (15b)
where $`\stackrel{}{\iota }`$ denotes the Poincaré variables and $`\stackrel{}{b}`$ belongs to the boundary of the Lobachevsky plane, i.e. $`b^2=1`$; let us consider afterwards the transformation to $`(u,v)`$

$`\stackrel{}{\iota }={\displaystyle \frac{2}{\sqrt{3}}}`$ $`\left[\left(u+{\displaystyle \frac{1}{2}}\right)\stackrel{}{b}^{}+v\stackrel{}{b}\right],`$ (16)
$`v0,\mathrm{}<u<\mathrm{},\stackrel{}{b}=(0,1),\stackrel{}{b}^{}=(1,0).`$
In terms of the $`(u,v)`$ variables, the anisotropy parameters (11) read

$`Q_1(u,v)={\displaystyle \frac{u}{\delta }},Q_2(u,v)={\displaystyle \frac{1+u}{\delta }},Q_3(u,v)={\displaystyle \frac{u(u+1)+v^2}{\delta }},`$ (17)
with $`\delta (u,v)(u^2+u+1+v^2)`$. The quantitities in (17) reduce to the Kasner exponents (5) when evaluated at $`v=0`$.
This correspondence provides the link between the dynamics representation in continuous variables and the BKL map .


text:
# Volcanic forcing improves Atmosphere-Ocean Coupled General Circulation Model scaling performance

## Abstract

Recent Atmosphere-Ocean Coupled General Circulation Model (AOGCM) simulations of the twentieth century climate, which account for anthropogenic and natural forcings, make it possible to study the origin of long-term temperature correlations found in the observed records. We study ensemble experiments performed with the NCAR PCM for 10 different historical scenarios, including no forcings, greenhouse gas, sulfate aerosol, ozone, solar, volcanic forcing and various combinations, such as natural, anthropogenic and all forcings. We compare the scaling exponents characterizing the long-term correlations of the observed and simulated model data for 16 representative land stations and 16 sites in the Atlantic Ocean for these scenarios. We find that inclusion of volcanic forcing in the AOGCM considerably improves the PCM scaling behavior. The scenarios containing volcanic forcing are able to reproduce quite well the observed scaling exponents for the land with exponents around 0.65 independent of the station distance from the ocean. For the Atlantic Ocean, scenarios with the volcanic forcing slightly underestimate the observed persistence exhibiting an average exponent 0.74 instead of 0.85 for reconstructed data.

While many modelers and climatologists focus their studies on trends caused by natural and anthropogenic forcings during the twentieth century Stott ; Mitchell ; Meehl ; Ammann ; Santer ; Storch  we here focus on another important aspect of temperature anomalies - long-term correlations. One of the first studies on the question how general circulation models reproduce observed climate variability was performed by Manabe . Using power spectrum analysis, which is affected by nonstationarities in time series, they argue that the GFDL AOGCM reproduced the natural climate variability on decadal and centennial scales correctly for a 1000 year control run integration. Several recent studies of KBunde1 ; KBunde2 ; Pelletier2 ; Talkner ; Pelletier1 ; Caballero ; Eichner  clearly demonstrate that surface air temperature (SAT) anomalies are long-term correlated with a fluctuation exponent $`\alpha `$ close to 0.65. On the other hand, results of Bell ; Govindan1 ; Govindan2 ; Syroka ; Vjushin  indicate that AOGCMs underestimate surface air temperture (SAT) persistence for control run, greenhouse gas forcing only and greenhouse gas forcing plus sulfate aerosols scenarios. In contrast, recent studies of reconstructed data Fraedrich ; Blender  claim that inner continental regions do not show long-term correlations, and thus AOGCMs successfully reproduce the natural persistence for the control run and greenhouse gas forcing only scenarios. Recently, this claim was tested on observed records by Bunde1  with the finding that the SAT fluctuation exponents for continental sites do show long-term correlation and the $`\alpha `$ values do not depend on the distance of the site from the ocean.

In this Letter we show that the recent PCM model simulations properly reproduce the observed long-term correlations for SAT on land only for those scenarios that include volcanic forcing. These scenarios also show better scaling performance over the ocean than the other scenarios.

In order to present the land-surface temperature profile for the last century, 16 observed daily maximum temperature time series are considered. They have been collected from different representative weather stations around the globe for the following sites: Vancouver, Tucson, Cheyenne, Luling, Brookings, Albany, Oxford, Prague, Kasan, Tashkent, Surgut, Chita, Seoul, Jakutsk, Melbourne and Sydney. We also analyze the gridded monthly mean sea surface temperature (SST) for 16 sites in the Atlantic ocean with a spatial resolution of 2.5ox2.5o for the period of $`19002002`$ from the Kaplan Extended SSTA data set (see also Monetti ). For $`19001981`$ this is the analysis of Kaplan  which uses optimal estimation in the space of 80 empirical orthogonal functions (EOFs) in order to interpolate ship observations of the U.K. Meteorological Office database Parker . The data after $`1981`$ consists of gridded data from the National Center for Environmental Prediction optimal interpolation analysis, which combines ship observations with remotely sensed data Reynolds . This analysis is performed on the same set of 80 EOFs as used in Kaplan  in order to provide enhanced data quality.

The model considered in our study is the Parallel Climate Model (PCM), which was developed at the National Center for Atmospheric Research (NCAR). It is a fully coupled global ocean-atmosphere-sea ice-land surface model that produces a stable climate without flux adjustment. The horizontal resolution of the atmosphere is equivalent to 2.8ox2.8o, with 18 levels in the vertical. Resolution of the ocean is roughly 2/3o, increasing to 1/2o at the equator, with 32 levels. The detailed description of the model and results from experiments using various forcings and their combinations may be found in Washington ; Dai ; Meehl  and Ammann .

Here we study 10 forcing combinations: no forcings, greenhouse gas, sulfate aerosol, ozone, solar, volcanic, solar + volcanic, ozone + solar + volcanic, greenhouse gas + sulfate aerosol + ozone, all forcings. Greenhouse gas forcing is based on historical observations of CO2, N2O, CH4, CFC-11, CFC-12, and ozone Dai . Evolution of direct forcing from tropospheric sulfate aerosol is reported by Kiehl . Historical changes of solar irradiance were reconstructed by Hoyt  and volcanic forcing by Ammann . The period of all experiments is 1890-1999.

For the no forcings and solar+volcanic scenarios we analyze the available 3-member ensembles, whereas for other scenarios 4-member ensembles are available. For each scenario, we selected the temperature records of the 4 grid points closest to each site, and bilinearly interpolated the data to the location of the observed site.

For each record, we analyse daily (or monthly) temperature anomalies $`\mathrm{\Delta }T_i`$. The $`\mathrm{\Delta }T_i`$ are called long-term correlated if their autocorrelation function $`C(s)`$ decays with time lag $`s`$ by a power law

$$C(s)s^\gamma ,0<\gamma <1.$$
(1)
To overcome possible nonstationarities in the data, we do not calculate $`C(s)`$ directly. Instead we construct the “profile” $`Y_n=_{i=1}^n\mathrm{\Delta }T_i`$ and study the fluctuation function $`F(s)`$ of the profile in segments of length $`s`$ by using the second order detrended fluctuation analysis (DFA2) Peng ; Jan . In DFA2 we determine in each segment the best second-order polynomial fit of the profile. The standard deviation of the profile from these polynomials represents the square of the fluctuations in each segment.

The fluctuation function $`F(s)`$ is the root mean square of the fluctuations in all segments. For the relevant case of long-term power-law correlations given by Eq. (1), with $`0<\gamma <1`$, the fluctuation function $`F(s)`$ increases according to a power law,

$$F(s)s^\alpha ,\alpha =1\frac{\gamma }{2}.$$
(2)
For uncorrelated data (as well as for short-range correlations represented by $`\gamma 1`$ or exponentially decaying correlation functions), we have $`\alpha =\frac{1}{2}`$. For long-term correlations we have $`\alpha >\frac{1}{2}`$.

First we plot the results of DFA2 (DFA curves of higher order show the same performance) for the observed daily maximum temperature (Figure 1a) and NCAR PCM simulations from the B06.61 run (Figure 1b). Run B06.61 represents one of the runs from the all forcings ensemble. All curves are shown in a double logarithmic presentation. We plot 10 typical DFA curves chosen from our 16 sites over land. The sites chosen include coastal, near coastal, and inland locations.

The approximate period of the observed records is 1880-1990, with the maximum length for Prague (1775-1992) and the minimum for Seoul (1908-1993). The period of the B06.61 run is 1890-1999. The slopes in Figure 1a correspond to fluctuation exponents of the observed SAT anomalies, and vary from 0.62 to 0.68, with an average close to 0.65. Figure 1a demonstrates that SAT anomalies for all sites studied obey long-term power-law correlations independent of the distance from the nearest ocean. The slopes in Figure 1b range from 0.62 to 0.69. Comparing Figures 1a and 1b shows that the scaling of the NCAR PCM output agrees quite well with the scaling of the observed data over land.

Figure 2 shows DFA2 curves for the 10 sites in the Atlantic ocean for the Kaplan reconstructed monthly SST anomalies (Figure 2a) and for the NCAR PCM monthly averaged SST anomalies from the all forcings B06.61 run (Figure 2b). The slopes for the reconstructed SSTA vary from 0.71 in the equatorial part of the Atlantic to 1.0 in the Northern Atlantic, with an average of 0.85. The SSTA exponents characterizing the memory effect on decadal and centennial scales seem to depend on complex ocean circulation dynamics. The variation of the scaling exponents over the Atlantic Ocean is significantly larger than on land, which is probably due to different ocean circulation patterns in equatorial, mid-latitude, and high-latitude regions. In a double logarithmic presentation, the slopes of the DFA2 curves for the simulated ocean records have an average of 0.72 which is noticeably lower than the observed average of 0.85.

Figure 3 presents the fluctuation exponent distribution for the observed data and for the 10 NCAR PCM scenarios considered for the 16 land locations. In each panel the “grey” column in the range \[0.62,0.68) represents the distribution of the fluctuation exponents for the observed data. As seen from the figure, scenarios containing volcanic forcing best reproduce the observations since they have a peak at the same range \[0.62-0.68) as the observed data. Their average fluctuation exponent is close to 0.65. In contrast, the other six scenarios that do not contain volcanic forcing, have an average fluctuation exponent for land less than 0.6 (see also Govindan2 ).

Similar behavior is found over the Atlantic Ocean. Only those scenarios containing the volcanic forcing exhibit an average fluctuation exponent greater than 0.7 with the largest value equal to 0.76 for the volcanic forcing only scenario. Figure 4 shows the fluctuation exponent distribution for the Atlantic Ocean for the Kaplan data (in grey) and for 10 studied scenarios (in black). Thus the PCM underestimates the fluctuation exponents obtained for reconstructed data for the Atlantic Ocean by 10-15%.

Therefore, we can conclude that for the NCAR PCM addition of volcanic forcing to any other forcing combination immediately improves its scaling behavior both for the land and the ocean. This fact suggests that (besides the atmosphere-ocean coupling) the volcanic forcing is mostly responsible for the presence of the long-term correlations in the NCAR PCM over land on annual and decadal scales. For the ocean, the addition of volcanic forcing leads to stronger memory and consequently higher fluctuation exponents comparing to those for the land. However the NCAR PCM still underestimates the observed persistence of the oceans.

The main conclusion from our research is that the NCAR PCM is able to reproduce the scaling behavior of the observed land SAT records for the last century after taking into account all historically based natural and anthropogenic forcings. However, even the best scenario for the land slightly underestimates natural SST persistence in the ocean, possibly due to errors in the simulations of deep ocean circulation, the atmosphere-ocean interaction, and/or an insufficiently long spin up period of the ocean component of the AOGCM.

The results presented in this letter may also help to clarify the controversy about the values of SAT fluctuation exponents for inner continental regions (see Fraedrich ). Our study indicates that not only the observed records for inner continental regions are long-term correlated in agreement with Bunde1 , but also the recent PCM simulations for these regions show similar fluctuation exponents, characteristic of long-term persistence.

Finally, this letter also supports the suggestion of Govindan2  that the inability of the seven leading AOGCMs, for their control runs, greenhouse gas forcing only, and greenhouse gas plus aerosols scenarios, to mimic the observed SAT persistence is caused by the absence of natural forcings, in particular volcanic forcing. The results of our detrended fluctuation analysis for the 16 land stations around the globe and the 16 sites in the Atlantic ocean suggest that volcanic forcing has by far the largest impact on the AOGCM long-memory persistence.

###### Acknowledgements.

This work has been supported by the Deutsche Forschungsgemeinschaft and the Israel Science Foundation. We are grateful to W. M. Washington, M. Wehner, and G. Strand for providing access to the NCAR PCM simulations data. The authors thank J. Eichner, J. Kantelhardt and E. Koscielny-Bunde for fruitful discussions.


text:
# Existence of Dynamical Scaling in the Temporal Signal of Time Projection Chamber

\[

## Abstract

The temporal signals from a large gas detector may show dynamical scaling due to many correlated space points created by the charged particles while passing through the tracking medium. This has been demonstrated through simulation using realistic parameters of a Time Projection Chamber (TPC) being fabricated to be used in ALICE collider experiment at CERN. An interesting aspect of this dynamical behavior is the existence of an universal scaling which does not depend on the multiplicity of the collision. This aspect can be utilised further to study physics at the device level and also for the online monitoring of certain physical observables including electronics noise which are a few crucial parameters for the optimal TPC performance.

\]

Dynamical scaling refers here to powerlaw distribution or a powerlaw with an exponential cut-off that describes certain correlation phenomena in the dynamical systems distinctly being different from the random statistical processes. In the thermodynamical context, it describes a critical phenomena associated with the phase transition  while in many other complex systems it corresponds to the so called self organised criticality . Although a powerlaw distribution $`(fx^\beta )`$ is a common dynamical feature, a powerlaw with an additional exponential $`(fx^\beta e^{\alpha x})`$ makes the distribution normalizable for all values of $`\beta `$ (where $`\alpha `$ and $`\beta `$ are positive constants) and also many real world systems like World Wide Web and social networks show this cut-off . In this paper, we show the presence of both type of scalings, a powerlaw and a powerlaw with exponential in the temporal signals (comprising of time gap and bunch length distributions) of a large Time Projection Chamber (TPC) which is a type of gas detector used for three dimensional tracking of the charged particles passing through it during a particle physics experiment. The time gap distribution along the drift direction shows a dynamical scaling which is independent of the multiplicity of the collisions above a critical value. While this scaling behavior itself is an interesting aspect to investigate the correlation phenomena in a gas detector at the device level, it can be further utilised for online monitoring of certain important TPC parameters including the constant rise in electronics noise level due to long exposure to radiation.

The TPC is the main tracking device of the ALICE (A Large Ion Collider Experiment)  at the Large Hadron Collider (LHC) at CERN optimized for the study of heavy ion collisions at a centre of mass energy $`5.5`$ ATeV. It is a large gas filled detector of cylindrical design with an inner radius of about 80 cm, an outer radius of about $`250`$ cm, and an over all length in the beam direction of $`500`$ cm. A charged particle passing through the gas volume creates electrons by ionization. The electrons drift in the electric field towards the read out chambers (multiwire proportional counters with more than $`550000`$ cathode pads read out located at the two end-caps of the TPC cylinder) where they are amplified in the field of the sense wires by a factor of several $`10^4`$. This signal is coupled to the read out pads which are on ground potential and at a few milimeters distance behind the sense wire-plane. The detail aspect of the design and simulation results can be found in ref . For simulation, the charged particles (mostly pions) with different multiplicities are generated using HIJING parametrization corresponding to $`Pb+Pb`$ collisons at 5.5 ATeV. The simulation is carried out using a microscopic simulator  incorporated in ALIROOT which is a GEANT3.21 and ROOT based simulation package used by the ALICE collaboration .

While we refer to  for deatil, in the following we briefly mention the salient features of the simulation. The ionization in the gas proceeds in two stages. Firstly, the electromagnetic interactions of the primary particles with the TPC gas (90$`\%`$ Ne and 10$`\%`$ CO2) lead to the release of primary electrons with a statistics that follows a Poisson distribution. Thus, the distance $`S`$ between two successive collisons leading to primary ionization can be simulated through an exponential distribution given by $`\mathrm{exp}(S/D)/D`$ where $`D=(N_{prim}f)^1`$ is the mean distance between primary ionizations, $`N_{prim}`$ is the number of primary electrons per cm produced by a Minimum Ionizing Particle (MIP) and $`f`$ is the Bethe-Bloch curve. At sufficient kinetic energy, the primary electrons produce secondaries creating an electron cluster with a total number of electrons given by $`N_{tot}=(E_{tot}I_{pot})W_i+1`$ where $`E_{tot}`$ is the energy loss in a given collision, $`W_i`$ is the effective energy required to produce an electron-ion pair and $`I_{pot}`$ is the first ionization potential. The electron cluster which is assumed to be point like undergoes diffusion while drifting towards the end-cap which is described by a three dimensional Gaussian with widths $`\delta _T=D_T\sqrt{L_d}`$ and $`\delta _L=D_L\sqrt{L_d}`$ where $`D_T`$ and $`D_L`$ are transverse and longitudinal diffusion constants and $`L_d`$ is the total drift distance. An electron arriving at the anode wire creates an avalanche which induces a charge on the pad plane. The time signal is obtained by folding the avalanche with the shaping function of the pre-amplifier/shaper with a shaping time $`200`$ ns which is a compromise between the need for achieving a high signal to noise ratio and for avoiding overlap of successive signals. This signal is sampled with a frequency $`5.66MHz`$ which divides the total drift time of $`88\mu s`$ into about $`500`$ time bins. The microscopic simulator also takes into account the electron loss in the drift gas due to presence of electron negative gas like O2 and also $`\mathrm{𝐄𝐱𝐁}`$ effect near the anode wires. The signal is digitized by a $`10`$ bit A/D converter that generates Gaussian random noise with r.m.s. about $`1000e`$. Finally, the digitized data is processed and formatted by an Application Specific Integrated Circuit (ASIC) called ALTRO (ALICE TPC Read Out) . A few typical parameters which are used in the simulations are given in table I. Although, we concentrate only on the TPC data, all component of the ALICE detectors as well as all passive materials are included in the simulation so as to create a situation close to the real experiment.

In ALTRO data format, zero suppressed data is recorded for each pad over all the time bins. This means, if we call $`bunch`$ a group of adjacent over threshold samples coming from one pad, the signal can be represented bunch by bunch. Figure 1 shows a typical plot of time bin versus amplitude for a given pad. Since the data is zero suppressed, it is sufficient to record three types of data, the sample amplitude in a given bunch, the bunch length and the time gap between two consecutive bunches. The amplitude distribution gives the energy loss spectrum with a long Landau tail and is not illuminating for the present purpose. Therefore, in the following, we will consider only the time gap and bunch length distributions built over all the pads.

The time gap distribution corresponds to the distance between two trajectories along the drift direction (say z-direction). Figure 2 shows the plot of time gap distribution at different multiplicities ($`M=20000`$, $`40000`$, $`60000`$ and $`80000`$). Although, the tail of the distribution is linear (in the semi-log scale) and depends on $`M`$ as expected, it deviates from the linearity at the shorter distances. The above behavior can be described by a distribution of the type,

$$f(z)=Az^\beta e^{\alpha z}forz1$$
(1)
The constant $`A`$ is fixed by the requirement of normalization which gives $`A=[Li_\beta (e^\alpha )]^1`$ where $`Li_n(x)`$ is the $`n^{th}`$ polylogarithm of $`x`$ given by,

$$Li_n(x)=\underset{i=1}{\overset{\mathrm{}}{}}\frac{x^i}{i^n}$$
(2)
Defining $`\overline{z}=_{z=1}^{\mathrm{}}zf(z)/_{z=1}^{\mathrm{}}f(z)`$, the average of the above distribution is

$$\overline{z}=\frac{Li_{\beta 1}(e^\alpha )}{Li_\beta (e^\alpha )}$$
(3)
Note that for $`\beta 0`$, the above distribution becomes a pure exponential with $`\overline{z}=(1e^\alpha )^1`$.

Figure 3 shows the parameters $`\beta `$ (filled circles in the left) and $`\alpha `$ (filled circles on the right) extracted from fitting the simulated data points. Note that $`\alpha `$ is a statistical parameter which has a near linear dependence on the multiplicty $`M`$. However, the index $`\beta `$ increases with $`M`$ and becomes a (nearly) constant at higer $`M`$ values. We will show below that this index is responsible for the dynamical scaling which does not depend on external parameters like multiplicity $`M`$ or magnetic field $`B`$ above a critical value but depends on certain intrinsic TPC parameters. Another interesting observation is that the product $`\omega =\overline{z}\alpha `$ also changes very slowly with $`M`$ above a critical value of $`M30000`$. In order to appreciate the effect of dynamical scaling, we can plot figure 2 in a reduced scale by dividing $`z`$ by $`\overline{z}`$ and multiplying $`f`$ by $`\overline{z}/A`$ in Eq.(1). By this rescaling, the index $`\beta `$ will not change, but $`\alpha `$ will be rescaled to say $`\omega =\alpha \overline{z}`$. Since $`\beta `$ and $`\omega `$ are (approximately) independent of $`M`$ above a certain value, the rescaled distributions will merge with each other as shown in figure 4. The solid curve is the fit of the type $`f(z)z^\beta e^{\omega z}`$ with $`\beta 0.45`$ and $`\omega 0.67`$. Recall here that this scaling is quite similar to the famous KNO scaling  which says that at high energies $`s`$, the probability distributions $`P_n(s)`$ of producing $`n`$ particles in a certain collision process should exhibit the scaling relation $`\overline{n}P_n(s)=f(n/\overline{n})`$. This scaling hypothesis asserts that if we rescale $`P_n(s)`$ measured at different energies via streching (shrinking) the vertical (horizontal) axes by $`\overline{n}`$, the rescaled curves will coincide with each other. In the present case, the observed scaling is identical to KNO scaling if we replace energy $`s`$ by the multiplicity $`M`$ and $`\overline{n}`$ by average distance $`\overline{z}`$ although both corresponds to two different physical situations.

Apart from this striking similarity with KNO, this rescaling is an elegant way to remove the statistical dependency from the dynamical behavior. It may be mentioned here that such type of scaling is also expected in case of an exponential distribution when the exponent $`\alpha `$ is small . Under such limit ($`\alpha <<1`$), since the average of an exponential distribution $`\overline{z}=(1e^\alpha )^1\alpha ^1`$, $`\omega =\alpha \overline{z}1`$ (The $`\omega `$ value may become more than unity if $`\alpha `$ is large). However, the interesting aspect of the present scaling is the deviation from unity and also having same $`\omega `$ at all $`M`$. This means $`\omega `$ deviates from unity due to presence of the dynamical exponent $`\beta `$. In absence of $`\beta `$, the scaling would have followed the linear behaviour in the semi-log scale as shown by the dashed curve in figure 4. Therefore, both $`\beta `$ and $`\eta =1\omega `$ are dynamical exponents whose values do not depend on $`M`$.

In the following, we investigate the origin of this dynamical behavior and also the parameters that affect its values. The origin of this dynamical phenomena can be associated with the set of measurements which are strongly correlated. Although, at a given interaction point, the creation of primary electrons due to ionization is a random statistical process, the set of interaction points created by a single charged particles passing through the TPC are well correlated. Obviously, this correlation will depend on the multiplicty $`M`$ and also on the number of measurements (number of pads, size etc). Depending on these parameters, a critical value is reached beyond which the correlation becomes independent of $`M`$.

Although, the dynamical behavior may depend on other geometrical TPC parameters, in this study, we do not intend to change any such parameters as they have been optimized based on certain physics criteria. However, we can see the effect of other relatively softer parameters like diffusion constants, drift velocity, noise etc which are likely to change during operation. It is noticed that out of all dynamical parameters, $`\eta `$ or $`(1\omega )`$ has some dependence on $`D_T`$ and strong dependence on electronic noise. Since the produced electron clouds are broadened in transverse direction due to increased $`D_T`$, the correlation effect is also enhanced. It is found that increasing $`D_T`$ by a factor of two, $`\omega `$ reduces only from $`0.67`$ to $`0.62`$. Further increase in $`D_T`$ has very little effect on $`\omega `$.

On the other hand, noise above the threshold reduces the correlation effect as it acts as an additional random source of electrons. As shown in figure 5, the middle curve is obtained with the parameters as given in the tablei (reffered as default parameters), while the upper curve corresponds to $`D_T`$ twice the default value and the lower curve is due to an increased noise level from $`1000e`$ to $`2000e`$. The correlation effect is lost due to noisy signal for which $`\omega 1`$. With increased noise level $`\omega `$ rises above unity (slowly) which is typical of an exponential behavior. Note that the dynamical parameters like $`D_T`$, $`D_L`$ and drift velocity etc are expected to vary within $`<4\%`$ during operation. This small fluctuation does not effect the $`\omega `$ value. On the otherhand, noise is quite unpredictable and also may go up due to the exposure of the electronics to the radiation for an extended period. Therefore, the parameter $`\omega `$ can serve as an excellent online tool to monitor the noise level.

So far, we have discussed only about the time gap distribution. The tail of the bunch length distribution also shows a powerlaw behaviour $`An^\gamma `$ as shown in figure 6(a). Since the tail of the bunch length distribution corresponds to low energy electrons and delta rays, the exponent $`\gamma `$ is quite sensitive to the applied magnetic field (up to some limit). Figure 6(b) shows the bunch length distributions at different $`B`$ values. May be this exponent can be used to monitor the magnetic field setting during the operation. Apart from this, we have not found the dependency of $`\gamma `$ on any other dynamical TPC parameters which could have been utilized for online monitoring.

We would like to add here that the powerlaw exponent that affect the time gap distribution corresponds to a region of short distances. Therefore, it is very difficult to extract the $`\beta `$ parameters acurately and also fitting the time gap distribution with other functional forms like two exponentials can not be ruled out. However, we have seen that the quality of fit is much better with powerlaw exponent. Further, as we have argued before, dynamical phenomena is expected in side a TPC due to many correlated interaction points. The tail of the bunch length distribution having a perfect powerlaw behavior reflects this aspect rather unambiguously. Due to the same dynamical origin, it is reasonable to assume a powerlaw with exponential cutoff for time gap distribution as well. We would also like to add here that the KNO type of scaling found in case of time gap distribution is independent of any fitting procedure and the deviation of the slope from unity is an indication that the noise level has remained below the threshold. This we consider as an important observation of this analyses irrespective of the fitting procedures.

In conclusion, the time gap distribution shows an universal scaling behavior (KNO type) (at sufficiently large multiplicity) with an exponent deviating from unity. This is an interesting aspect to study the correlation phenomena in gas detector and also to understand physics (more explicitly, the physics that has gone into simulation) at the device level. An important pratical utility of this phenomena is the utilization of the above scaling exponent to monitor the noise level above a given threshold without any computational complexity and also without building any rigorous models. In that sense, this analyses provide a model independent way to monitor the quality of the data what is being recorded.