format: text

text:
# On the Angular Correlation Function of SZ Clusters : Extracting cosmological information from a 2D catalog

## 1 Introduction

The Sunyaev–Zel’dovich (SZ) effect (Sunyaev & Zel’dovich 1970, 1972) has become a practical observational tool for studying galaxy clusters and cosmology (for recent reviews see Birkinshaw 1999; Carlstrom et al. 2002). Current observations of individual clusters, when combined with X–ray observations, constrain cosmological parameters via gas mass fractions (Grego et al. 2002) and angular–diameter distance determinations (Mason et al. 2001; Jones et al. 2001; Reese et al. 2002). Multi–band millimeter observations of a handful of clusters have already been used to set limits on peculiar velocities (Holzapfel et al. 1997; Benson et al. 2003), and theoretical studies of this technique show its promise for the future (Aghanim et al. 2002a, 2002b; Holder 2002). A new generation of optimized, dedicated instruments, both large–format bolometer arrays and interferometers with high sensitivity receivers, will qualitatively improve these studies. And the arrival of these instruments within the next few years, in addition to the Planck mission111http://astro.estec.esa.nl/SA-general/Projects/Planck/, will move the field forward to its next important step: surveying. This will open a new observational window onto large–scale structure and its evolution out to large redshifts.

The ultimate goal of these SZ surveys is the construction of large cluster catalogs with multi–wavelength follow–up observations in order to perform cosmological studies; for example, constraining cosmological parameters with the counts and redshift distribution (e.g., Barbosa et al. 1996; Eke et al. 1996; Colafrancesco et al. 1997; Haiman et al. 2001; Holder et al. 2001; Kneissl et al. 2001; Weller et al. 2002; Benson et al. 2002). Driving this effort are the particular advantages of SZ–based cluster catalogs (Bartlett 2000; Bartlett 2001): Firstly, SZ surveys are intrinsically efficient at finding clusters at large redshift, due to the surface brightness constancy of the SZ effect222for a cluster of fixed properties. The thermal SZ spectrum is furthermore universal, the same for all clusters at any redshift333in the non–relativistic limit. Other emission mechanisms, in contrast, suffer from cosmological dimming and the need for accurate $`k`$–corrections. Secondly, SZ surveys select clusters based on their thermal energy. Since the spectrum is the same for all clusters, the total observable SZ flux from a cluster can be expressed in a frequency independent manner as the integrated Compton $`y`$–parameter, $`Y=𝑑\mathrm{\Omega }y(\widehat{n})`$, where the integral is over the cluster profile (see Equation 4 below). The $`y`$ parameter being the pressure integrated along the line–of–sight ($`ydlnT)`$, this then implies $`YMT`$, i.e., the thermal energy of the intracluster medium (ICM). This is important, because the total thermal energy of the ICM is given by energy re-partition during cluster collapse and is independent of any thermal or spatial structure in the gas. It is hence a more robust quantity than, for example, the X–ray emission measure that depends in a more complicated fashion on both the ICM density and temperature. Hydrodynamical simulations confirm this expectation by showing a tight SZ flux–mass relationship with little scatter (da Silva et al. 2003). Object selection in a flux–limited SZ catalog is therefore relatively easy to interpret in terms of cluster mass and redshift. For instance, it is easy to show that the minimum detectable cluster mass is almost independent of redshift. This is particularly advantageous for evolutionary studies, because one is able to follow the evolution of the same kind of object over redshift, instead of comparing massive objects at high redshift to less massive ones at low redshift, as is the case with X–ray samples.

Detailed follow–up of SZ surveys will, however, be time–consuming, and an enormous effort for the more than $`10^4`$ clusters expected from Planck. Large–area photometric surveys in the optical and infrared (e.g., SDSS) will help (Bartelmann & White 2002), but it is important to identify the kind of science that may be done directly with the two dimensional catalog of cluster positions and SZ fluxes, what we will refer to in the following as the SZ photometric catalog. This will certainly be the first science to be performed. Source counts represent the primary avenue of 2D study that has been discussed extensively in the literature. In this paper, we examine the next higher order catalog statistic, namely, the angular correlation function $`w(\theta )`$ of SZ–detected clusters (Diaferio et al. 2003). We quantify its information content and study its potential use as a cosmological probe. The angular function samples the catalog redshift distribution, because it is a projection of the real–space correlation function along the line–of–sight. With an appropriate model for the real–space correlation function of the catalog, we may gain some insight on this distribution, and hence on the underlying cosmological model.

Moscardini et al. (2002) studied the 3D clustering properties of SZ detected clusters, taking the Planck survey as an example and accounting for evolutionary effects along the past light–cone. Diaferio et al. (2003) examined the angular function of SZ clusters as a means of identifying probable physical cluster pairs (in 3D) and superclusters. Our modeling is similar to theirs, although we focus instead on the angular function as a cosmological probe permitting the extraction of cosmological information from a photometric SZ catalog. The idea is of course not new, and has been applied in the past to, for example, optical and X–ray cluster catalogs; but we reiterate the advantages of an SZ catalog in this context: the cluster selection function is relatively easy to model (compared to other observing bands), and it extends out to large redshift, giving a longer base–line for viewing evolutionary effects.

We should distinguish at the outset the difference between the angular power spectrum, $`C_l^{\mathrm{sz}}`$, of SZ–induced temperature fluctuations (secondary anisotropies in the cosmic microwave background \[CMB\]) and the angular correlation function of detected clusters in a SZ survey, $`w(\theta )`$. The angular power spectrum $`C_l^{\mathrm{sz}}`$ is a two–point statistic quantifying the integrated contribution of the entire cluster population to the CMB sky brightness fluctuations. It is dominated by the Poisson term and its overall shape is determined by the mean SZ profile of clusters. Cluster–cluster correlations add additional power on the order of $`2030`$% of the pure Poisson term (Komatsu & Kitayama 1999). Since it is defined relative to the mean cosmic microwave background temperature, we expect the SZ fluctuation power to increase with the surface density of clusters on the sky. This is quantitatively confirmed by both numerical simulations and analytical calculations that indicate $`C_l^{\mathrm{sz}}\sigma _8^7`$ (all other factors held constant), where $`\sigma _8`$ is the amplitude the density perturbations, the quantity most directly influencing cluster abundance. The fluctuation power spectrum is an analysis method appropriate in a low signal–to–noise context (the current situation) where individual source identification is not possible444Either for low signal–to–noise observations or when pushing constraints on source counts below the detection threshold. It is the SZ equivalent of the “background fluctuation analysis” of radio and X–ray Astronomy.. Fluctuations induced by the SZ effect have been invoked as a possible explanation for the excess power at high multipole $`l`$ reported by the CBI collaboration (Bond et al. 2002) and consistent with new VSA data (Grainge et al. 2002). If this were entirely due to the cluster population, it would imply a surprisingly large value for $`\sigma _8`$ ($`>1`$; Bond et al. 2002; Komatsu & Seljak 2002; Holder 2002), although an important contribution from heated gas at reionization may also be expected (Oh et al. 2003).

The angular function $`w(\theta )`$ quantifies the projected clustering of a 2D catalog of individually detected clusters. It refers to the object positions and makes no reference to the mean background sky brightness. There are several ways to imagine using the information contained in the SZ cluster angular function (e.g., Diaferio et al. 2003). In the following we choose to illustrate its use by examining constraints on the matter density $`\mathrm{\Omega }_\mathrm{M}`$ and $`\sigma _8`$ in the context of flat CDM–like models. The SZ counts provide one constraint on a combination of these parameters. To extract additional information from the 2D catalog using the angular function, we are forced to model the real–space cluster correlation function. In CDM scenarios, clusters form from peaks in the density field whose clustering may be analytically calculated. We adopt the approach proposed by Mo & White (1996; 2002). Any conclusions that we draw are, therefore, unavoidably dependent on this clustering model (Moscardini et al. 2002; Diaferio et al. 2003); however, it is well founded in the context of CDM cosmogonies and compares well with the results of numerical simulations, at least at redshifts lower than $`10`$ (Reed et al. 2003; Jenkins et al. 2001). Another important issue concerns the modeling of the intracluster medium. Moscardini et al. (2002), for example, find non–negligible dependence of the 3D correlation function on ICM properties. We discuss the issue below for the angular function, but we note here that it tends to be dominated by massive clusters (depending of course on the exact flux limit; see below) that follow relatively well observed scaling laws. As we show, adding constraints from the local X–ray cluster abundance permits us to simultaneously constrain the cosmological parameters and the cluster baryon content.

In the following section, we give our master equation for calculating $`w(\theta )`$ and identify the necessary modeling ingredients. We outline our cluster model in Section 3. Results for the angular function are presented in Section 4, where we apply the results to constraining the cosmological parameters $`\sigma _8`$ and $`\mathrm{\Omega }_\mathrm{M}`$. We discuss the influence of ICM physics on the results and how to use additional information from the local cluster abundance to simultaneously constrain this physics and the cosmology. Our fiducial example is the Planck mission. Section 5 closes with a final discussion and summary.

## 2 The Angular Correlation Function

In this section we relate the angular correlation function to the real–space correlation function in the context of SZ observations. The 3–dimensional (auto)correlation function $`\xi `$ quantifies the 2–point clustering of a population in terms of the probability in excess of Poisson of finding two objects at a separation $`r`$. The angular correlation $`w(\theta )`$ of the same population is the projection of $`\xi `$ onto the sky:

$`d\mathrm{\Omega }_1d\mathrm{\Omega }_2\mathrm{\Sigma }^2[1+w(\theta )]=`$ (1)
$`d\mathrm{\Omega }_1d\mathrm{\Omega }_2{\displaystyle }dz_1dz_2{\displaystyle \frac{dV_1}{d\mathrm{\Omega }_1dz_1}}{\displaystyle \frac{dV_2}{d\mathrm{\Omega }_2dz_2}}\times `$
$`{\displaystyle 𝑑M_1𝑑M_2\frac{dn_1}{dM_1}\frac{dn_2}{dM_2}(1+\xi )}`$
where the integrals concern two lines–of–sight (los) of solid angles $`d\mathrm{\Omega }_1`$ and $`d\mathrm{\Omega }_2`$ separated on the sky by angle $`\theta `$. In this expression, $`dn/dM`$ represents the cluster mass function, and we assume the small angle approximation here and throughout. The surface density of sources (the counts) with integrated Compton parameter larger than $`Y_{\mathrm{lim}}`$ (the SZ ‘flux’ limit; see the following section, Equation 4) is given by

$$\mathrm{\Sigma }(Y_{\mathrm{lim}})=_0^{\mathrm{}}𝑑z\frac{dV}{d\mathrm{\Omega }dz}_{M_{\mathrm{lim}}(Y_{\mathrm{lim}},z)}^{\mathrm{}}𝑑M\frac{dn}{dM}(M,z)$$
(2)
The correlation function $`\xi `$ depends on cluster mass, redshift and, according to statistical isotropy, on physical separation $`r(\theta )=r_1^2(z_1)+r_2^2(z_2)2r_1r_2\mathrm{cos}(\theta )`$, where $`r_1`$ and $`r_2`$ are angular–diameter distances. Assuming that correlations fall off sufficiently rapidly with distance, as is observed, we may take the two clusters to be at approximately the same redshift and write $`\xi [M_1,M_2,z,r(\theta )]`$. We furthermore adopt a linear biasing scheme in which $`\xi (M_1,M_2,z,r)=b(M_1,z)b(M_2,z)\xi _{\mathrm{dm}}(z,r)`$, where $`\xi _{\mathrm{dm}}(z,r)`$ is the correlation function of the underlying cold dark matter and $`b(M,z)`$ is the bias factor for clusters (see below). Then, using the short–hand notation $`\varphi (M,z,Y_{\mathrm{lim}})=(1/\mathrm{\Sigma })b(M,z)(dV/dzd\mathrm{\Omega })dn/dM`$ for the joint distribution of clusters in mass and redshift, weighted by the bias factor, we arrive at the expression

$`w(\theta ,Y_{\mathrm{lim}})=`$
$`{\displaystyle _0^{\mathrm{}}}dz{\displaystyle }{\displaystyle _{M_{\mathrm{lim}}(Y_{\mathrm{lim}},z)}^{\mathrm{}}}dM_1dM_2\varphi (M_1,z)\varphi (M_2,z)\times `$
$`{\displaystyle 𝑑r\frac{dz}{dr}(z)\xi _{\mathrm{dm}}[z,r(\theta )]}`$
$`{\displaystyle _0^{\mathrm{}}}𝑑z\mathrm{\Phi }^2(z,Y_{\mathrm{lim}}){\displaystyle 𝑑r\frac{dz}{dr}(z)\xi _{\mathrm{dm}}[z,r(\theta )]}`$
a sort of Limber’s equation appropriate for SZ sources (Limber 1953; Peebles 1993; Diaferio et al. 2003) in which we explicitly show the dependence on limiting flux $`Y_{\mathrm{lim}}`$. From this equation we clearly see that the three key ingredients are the the mass–limit function $`M_{\mathrm{lim}}(z,Y_{\mathrm{lim}})`$, the distribution function $`\varphi `$ and the correlation function $`\xi _{\mathrm{dm}}`$. We now discuss our modeling of each.

## 3 The SZ Population

The SZ cluster population inherits its properties from two sources: its constituent dark matter halos, whose properties are the sole result of gravitational evolution, and the relationship between observable SZ flux and these halos, governed by more difficult to model baryonic physics. It is reasonable to characterize dark halos by their mass and redshift, and we will apply the results of N–body simulations that give both their abundance (i.e., the mass function) and spatial correlations (i.e., clustering bias $`b(M,z)`$ and the correlation function $`\xi _{\mathrm{dm}}`$) as a function of these two fundamental descriptors. More difficult to model, the baryonic physics of the cluster gas requires particular attention to various (and at times contradictory) observational constraints and theoretical scaling laws.

### 3.1 Halo properties

The abundance of galaxy clusters is given by the mass function of collapsed objects, which is completely specified once the linear power spectrum of dark matter perturbations is specified. For the latter we adopt the BBKS (Bardeen et al. 1986) transfer function (see also below), while for the mass function we employ the fitting formula (improved Press–Schechter) given by Sheth & Tormen (1999):

$`{\displaystyle \frac{dn}{dM}}(M,z)dM=A\left(1+{\displaystyle \frac{1}{\nu ^{2q}}}\right)\sqrt{{\displaystyle \frac{2}{\pi }}}{\displaystyle \frac{\overline{\rho }}{M}}{\displaystyle \frac{d\nu ^{}}{dM}}\mathrm{exp}\left({\displaystyle \frac{\nu ^2}{2}}\right)dM`$ (3)
where $`\overline{\rho }`$ is the universal mean mass density and the constants $`A0.322`$ and $`q=0.3`$; the parameter $`\nu ^{}=\sqrt{a}\nu `$, where $`\nu \frac{\delta _c}{D(z)\sigma (M)}`$ is the usual critical peak height ($`\delta _c1.69`$) normalized to the mass–density perturbation variance $`\sigma (M)`$ in spheres containing mass $`M`$, and the constant $`a=0.707`$. The expression for the growth factor for flat models with $`\mathrm{\Lambda }>0`$ is taken from Carroll et al. (1992):

$`D(z,\mathrm{\Omega }_\mathrm{M},\mathrm{\Omega }_\mathrm{\Lambda })=`$ $`{\displaystyle \frac{g(z)}{g(0)(1+z)}}`$
$`g(z,\mathrm{\Omega }_\mathrm{M},\mathrm{\Omega }_\mathrm{\Lambda })`$ $`{\displaystyle \frac{5}{2}}\mathrm{\Omega }_\mathrm{M}(z)[\mathrm{\Omega }_\mathrm{M}(z)^{4/7}`$
$`\mathrm{\Omega }_\mathrm{\Lambda }(z)`$ $`+(1+\mathrm{\Omega }_\mathrm{M}(z)/2)(1+\mathrm{\Omega }_\mathrm{\Lambda }(z)/70)]^1`$
with the definitions $`\mathrm{\Omega }_\mathrm{M}(z)\mathrm{\Omega }_\mathrm{M}(1+z)^3/E^2(z)`$, $`\mathrm{\Omega }_\mathrm{\Lambda }(z)\mathrm{\Omega }_\mathrm{\Lambda }/E^2(z)`$, and $`E^2(z)=[\mathrm{\Omega }_\mathrm{\Lambda }+(1\mathrm{\Omega }_\mathrm{M}\mathrm{\Omega }_\mathrm{\Lambda })(1+z)^2+\mathrm{\Omega }_\mathrm{M}(1+z)^3]`$; $`\mathrm{\Omega }_M`$ and $`\mathrm{\Omega }_\mathrm{\Lambda }`$ written without an explicit redshift dependence will indicate present–day values ($`z=0`$).

As mentioned above, we use a linear bias scheme to relate the cluster–cluster correlation function to that of the dark matter and employ the analytic fitting formula for $`b(M,z)`$ given by Sheth et al. (2001); the formula includes corrections for ellipsoidal perturbation collapse:

$`b(M,z)=1+{\displaystyle \frac{1}{\sqrt{a}\delta _c(z)}}\left[\sqrt{a}(a\nu ^2)+\sqrt{a}b(a\nu ^2)^{1c}\right]`$
$`{\displaystyle \frac{(a\nu ^2)^c}{(a\nu ^2)^c+b(1c)(1c/2)}}`$
where $`\delta _c`$, $`\nu `$ and $`a`$ are given above, and $`b=0.5`$ and $`c=0.6`$. We model the linear dark matter perturbation spectrum with the BBKS transfer function with shape parameter fixed at $`\mathrm{\Gamma }=0.25`$ and scale–invariant primordial density fluctuations ($`n=1`$). This seems to provide a good fit to galaxy clustering data (Percival et al. 2001) and is consistent with constraints on $`n`$ from CMB anisotropies (Spergel et al. 2003, and references therein). The resulting linear theory $`\xi _{\mathrm{dm}}`$ is adequate on most scales ($`\theta >10\mathrm{}`$), although we also include non–linear corrections according to the fitting formula developed by Peacock & Dodds (1996).

### 3.2 Intracluster Medium

With the abundance and clustering of halos now specified, we next relate the observable SZ flux to cluster mass and redshift. This relation is particularly robust from a theoretical viewpoint, contrary to, for example, the situation for X–ray luminosity. In the non–relativistic regime, the surface brightness of the thermal SZ effect – measured relative to the mean sky intensity – at position $`\widehat{\mathrm{\Omega }}`$ on a cluster image is

$$\mathrm{\Delta }i_\nu (\widehat{\mathrm{\Omega }})=y(\widehat{\mathrm{\Omega }})j_\nu $$

where the Compton $`y`$ parameter is an integral of the pressure along the line–of–sight

$$y(\widehat{\mathrm{\Omega }})=_{los\widehat{\mathrm{\Omega }}}𝑑l\frac{kT}{mc^2}n\sigma _T$$

with $`T`$ the gas temperature (strictly speaking, that of the electrons), $`\sigma _T`$ the Thompson cross section, $`k`$ and $`m`$ the Boltzmann constant and the electron mass, respectively, and where $`j_\nu `$ is a universal spectral function that is the same for all clusters, independent of their properties (Birkinshaw 1999). Since the thermal SZ spectrum is the same for all clusters, we may express the total flux in a frequency independent manner using the integrated Compton parameter $`Y=𝑑\mathrm{\Omega }y(\widehat{\mathrm{\Omega }})`$, where the integral is over the cluster profile. The total flux density (e.g., in Jy) is then $`S_\nu =Yj_\nu `$. This is the total flux (density) measured by an experiment with low angular resolution in which clusters are simply point sources. From these definitions, we find that the observable flux is directly proportional to the total thermal energy of the ICM (e.g., Barbosa et al. 1996):

$$Y(M,z)=\frac{k\sigma _T}{mc^2}\frac{N_\mathrm{e}T}{D_{\mathrm{ang}}^2(z)}\frac{f_{\mathrm{gas}}(M,z)T(M,z)M}{D_{\mathrm{ang}}^2(z)}$$
(4)
where $`N_\mathrm{e}`$ is the total number of electrons, $`f_{\mathrm{gas}}`$ is the ICM mass fraction, $`D_{\mathrm{ang}}(z)=H_\mathrm{o}^1d_{\mathrm{ang}}(z)`$ is the angular diameter distance, and $`T`$ is to be understood as the mean (particle, and not emission–weighted) electron temperature; note that with this understanding, the relation does not depend on any assumption of isothermality. It is this direct relation between observable SZ flux and thermal gas energy that lies at the heart of some of the advantages of SZ over X–ray surveys (Bartlett 2001).
We have taken care to write the quantities $`f_{\mathrm{gas}}`$ and $`T`$ as general functions of mass and redshift. In the absence of efficient heating/cooling and gas reprocessing, the cluster population will be fully self–similar. Simple theoretical arguments based on energetics during collapse suggest in this case the existence of a scaling law between cluster temperature and mass:

$$T(M,z)=T_{}\left(M_{15}h\right)^{2/3}\left(\mathrm{\Delta }(z)E(z)^2\right)^{1/3}\left[12\frac{\mathrm{\Omega }_\mathrm{\Lambda }(z)}{\mathrm{\Delta }(z)}\right]$$
(5)
adopting the notation of Pierpaoli et al. (2002). In this expression, the mass $`M_{15}`$ is measured in units of $`10^{15}`$M, $`\mathrm{\Delta }(z)`$ is the full non–linear overdensity inside the virial radius relative to the critical density ($`178`$). In the self–similar model, the gas mass fraction $`f_{\mathrm{gas}}`$ is constant, essentially proportional to the universal ratio $`\mathrm{\Omega }_\mathrm{B}/\mathrm{\Omega }_\mathrm{M}`$.
These scaling expectations are indeed manifest in hydrodynamical simulations that neglect cooling, and supported by observations of the more massive clusters with cooling timescales longer than the Hubble time (Mohr et al. 1999; Allen et al. 2001; Finoguenov et al. 2001; Xu et al. 2001; Arnaud et al. 2002). Although there is good indication that clusters are a population with rather regular properties, as suggested by these scaling laws, there is also direct evidence that it is not exactly self–similar. Deviations are most pronounced for the lower mass objects ($`T<2`$ keV), as perhaps expected since they have shorter cooling times and plausible energy injection mechanisms more readily compete with their gravitational energy (Ponman et al. 1999; Lloyd–Davies et al. 2000). More generally, one may write

$$Y(M,z)=Y_{15}(z)M_{15}^{5/3+\alpha }(1+z)^\gamma $$
(6)
The factor $`Y_{15}(z)`$ incorporates all the $`z`$ dependence of Eqs. (4,5) and is defined such that the self–similar model corresponds to $`\alpha =\gamma =0`$.

Although we have discussed Eq. (6) via scaling relations for $`T`$ and $`f_{\mathrm{gas}}`$, it is more pertinent to consider it as a single relation between observable SZ flux and cluster mass and redshift, one that expresses the proportionality between the observable flux and the gas thermal energy. Hydrodynamical simulations incorporating cooling and pre–heating (da Silva et al. 2003) indicate $`\alpha `$ between 0.1 and 0.2 and $`\gamma 0`$, with a very tight scatter about the relation, much tighter than corresponding relations for X–ray quantities. This confirms our expectations that the SZ flux is a more robust quantity than its X–ray counterparts. Even though very low mass clusters are included in the fit, only mild deviations from self–similarity in this relation are seen in the simulations. These deviations are most pronounced in the low mass systems, while the higher mass systems appear compatible with the self–similar scaling laws for $`T`$$`M`$ and $`f_{\mathrm{gas}}`$. Caution is still warranted, however, due to the unstable nature of cooling, issues of numerical resolution, and the fact that these codes model the gas as a single phase medium; the simulation volume also contains mostly low mass systems, with only a few clusters with $`M_{15}>0.2`$. Simulation results are nonetheless roughly consistent with X–ray observations (Borgani et al. 2002; Muanwong et al. 2002).

At our fiducial Planck flux limit of $`Y=10^4`$ arcmin2, we have checked that the counts and value of the angular function at 30 arcmins are little affected by clusters below $`T=2`$ keV, at least at $`\mathrm{\Omega }_M<0.6`$. We therefore concentrate on the self–similar case with $`\alpha =\gamma =0`$. On the other hand, we consider the normalization of Eq. (6) as a parameter free to vary within certain limits suggested by X–ray observations of $`f_{\mathrm{gas}}`$ and $`T`$. Moscardini et al. (2002) have argued that this freedom makes it difficult to use only SZ observations to constrain $`\sigma _8`$ and $`\mathrm{\Omega }_\mathrm{M}`$. To overcome this modeling uncertainty, we combine SZ observations of both the number counts and the angular correlation function with constraints arising from the local abundance of X–ray clusters. We shall find that the three kinds of observations are complementary and lead to constraints on the cosmological parameters.
To give a feel for the order of magnitude, we note that

$`Y_{15}(z)`$ $`=(7.4\times 10^5h^{7/6}\mathrm{arcmin}^2)\left({\displaystyle \frac{T_{}}{\mathrm{keV}}}\right)\left({\displaystyle \frac{f_{\mathrm{gas}}}{0.07h^{3/2}}}\right)\times `$ (7)
$`\left({\displaystyle \frac{\mathrm{\Delta }(z)E(z)^2}{178}}\right)^{1/3}\left[12{\displaystyle \frac{\mathrm{\Omega }_\mathrm{\Lambda }(z)}{\mathrm{\Delta }(z)}}\right]{\displaystyle \frac{1}{d_{\mathrm{ang}}^2(z)}}`$
$`Y_{}\left({\displaystyle \frac{\mathrm{\Delta }(z)E(z)^2}{178}}\right)^{1/3}\left[12{\displaystyle \frac{\mathrm{\Omega }_\mathrm{\Lambda }(z)}{\mathrm{\Delta }(z)}}\right]{\displaystyle \frac{1}{d_{\mathrm{ang}}^2(z)}}`$
For reference, $`T_{}=1.2`$ according to the simulations of Evrard et al. (1996) and $`f_{\mathrm{gas}}=0.07h^{1.5}`$ from Mohr et al. (1999). At an observation frequency of 2.1 mm, the maximum decrement of the thermal spectrum, a $`Y=7.4\times 10^5`$ arcmin2 corresponds to a flux density of $`7`$ mJy. In all the analysis a minimum cluster mass of $`10^{14}h^1`$M is imposed.

## 4 Results

### 4.1 The angular correlation function

Calculated angular correlation functions are shown in Figure 1 for three different limiting flux values ($`Y_{\mathrm{lim}}`$) in two different flat cosmologies ($`\mathrm{\Omega }_\mathrm{M}=0.3`$ and $`\mathrm{\Omega }_\mathrm{\Lambda }=0.7`$, with $`\sigma _8=0.9`$, and $`\mathrm{\Omega }_\mathrm{M}=1`$ with $`\sigma _8=0.7`$; e.g., Blanchard et al. 2000); the $`Y`$$`M`$ normalization is $`Y_{}=7.6\times 10^5h^{7/6}`$ arcmin2. Our angular correlation function is consistent with the results of Diaferio et al. (2003) at a separation of 30 arcmins. As shown by the latter authors, small scales are affected by the specific choice of clustering bias function $`b(M,z)`$, but the model differences drop to $`10`$ % at 30 arcmins. This is comparable to the statistical measurement errors expected in the case of the Planck survey, as discussed below. We therefore use the angular function at 30 arcmin separation in our analysis of cosmological parameters.
Figure 2 shows the cluster selection function $`\mathrm{\Phi }(z)`$, defined in Eq. (2), at $`Y_{\mathrm{lim}}=10^4`$ arcmin2 in the two cosmologies and for various normalizations $`\sigma _8`$. The broad break in the angular function in Figure 1 corresponds to the break in $`\xi _{\mathrm{dm}}`$ just beyond $`10h^1`$ Mpc. For example, the break occurs $`1\mathrm{deg}`$ for $`Y_{\mathrm{lim}}=10^4`$ arcmin2 in the critical model. According to Figure 2, the selection function peaks around $`z0.2`$ (for the chosen value of $`\sigma _8=0.7`$), projecting the break to an angular separation of $`\theta 10h^1/(3000h^1\times 0.2)1\mathrm{deg}`$. Figures 1 and 2 visually illustrate how the angular function encodes information on the catalog radial distribution.

Although the two models have comparable angular correlations at the bright end (the upper curves in Fig. 1 corresponding to $`Y_{\mathrm{lim}}=10^3`$ arcmin2), their dependence on catalog depth clearly differs. At the bright end, we are mainly observing the local cluster population, which is essentially the same in both models since the present–day abundance is the same and the density perturbation power spectrum is fixed ($`\mathrm{\Gamma }=0.25`$). Note that the low–density model has a slightly steeper slope at small separation, due to its greater non–linear evolution. The angular function of the low–density model decreases and shifts to the left more rapidly with survey depth than in the critical model. The overall trend is easily understandable and due to the fact that the selection function broadens and peaks at higher redshift with survey depth, moving correlations to smaller angular scales and generally washing out the signal as more clusters are projected along the line–of–sight. In bright galaxy surveys, which sample the local universe where space is approximately Euclidean and galaxy evolution may be ignored, the dependence of the angular function on survey depth follows an important scaling law that is independent of the underlying cosmological model. No such universal scaling law obtains in the SZ case, because the radial distribution extends out to large redshifts. For the relevant flux limits, a SZ catalog therefore samples evolution in the cluster population. Since the cluster population evolves less rapidly with redshift in the low–density model, the angular correlation function therefore shifts down more rapidly with survey depth than in the case of the critical model. The angular correlation function of SZ clusters is therefore a cosmological probe.

### 4.2 Combining the SZ angular function and counts

As an illustration of this probe, we next consider the dependence of the angular function and of the SZ counts on the cosmological parameters $`(\mathrm{\Omega }_\mathrm{M},\sigma _8)`$. In Figure 3 we show the predicted counts and angular correlation function for a set of flat models over a grid in the $`(\sigma _8,\mathrm{\Omega }_\mathrm{M})`$–plane and for different limiting flux values, $`Y_{\mathrm{lim}}=10^3`$ arcmin2 (top left), $`Y_{\mathrm{lim}}=10^4`$ arcmin2 (top right), and $`Y_{\mathrm{lim}}=10^5`$ arcmin2 (bottom left). Contours of the angular correlation function are shown as continuous lines and refer to its value at a separation of $`\theta =30`$ arcmins. The counts contours at the corresponding flux limits are given as dotted lines and labeled in units of 1/sq.deg. This figure is constructed for a self–similar cluster population with a $`Y`$$`M`$ normalization, $`Y_{}=7.6\times 10^5h^{7/6}`$ arcmin2 (see Eq. 7); the effects of varying ICM properties will be discussed below.
We see from the figure that the counts and angular function carry more complementary information as the catalog increases in depth. At the high flux limit of $`Y_{\mathrm{lim}}=10^3`$ arcmin2, the catalog primarily probes the local universe; the two sets of contours are parallel and there is no information permitting us to constrain the two cosmological parameters. As we move deeper, the two sets of contours begin to cross, indicating the we are obtaining useful cosmological information capable of constraining the two parameters. The limit $`Y_{\mathrm{lim}}=10^4`$ arcmin2 is representative of the Planck mission (Aghanim et al. 1997; Bartelmann 2001; and references therein), while $`Y_{\mathrm{lim}}=10^5`$ arcmin2 would correspond to deeper ground–based surveys (probably attaining the confusion limit).

### 4.3 Influence of the ICM

The single relation Eq. (6) incorporates all of the ICM physics, and from the nature of the SZ flux, we expect it to be rather tight with little scatter. There is virtually no direct observational information on this relation, while numerical simulations (da Silva et al. 2003) do indeed confirm a small scatter. They also indicate that there is only a slight deviation from the self–similarity with $`\alpha =0.10.2`$ and $`\gamma =0`$ down to masses well below $`10^{14}`$M. X–ray observations, on the other hand, demonstrate deviations from self–similarity. This is most clear for low–mass systems with temperatures below $`T2`$ keV, but even in the richer systems the observed $`L`$$`T`$ deviates from the self–similar expectation. These results are beginning to be reproduced by numerical simulations including cooling/heating (e.g., Borgani et al. 2002), including those that indicate only slight deviation from self–similarity for the SZ relation Eq. (6) (Muanwong et al. 2002). That X–ray observables show greater deviation from self–similarity than the SZ flux is not surprising, given that the former are more sensitive to spatial/temperature structure in the gas than the latter, as emphasized previously. This is precisely the advantage of SZ surveys.

As mentioned, low mass systems with the greatest deviation from self–similarity in X–rays contribute little to counts and the angular function at a flux limit of $`Y_{\mathrm{lim}}=10^4`$ arcmin2. For example, for $`\mathrm{\Omega }_\mathrm{M}<0.6`$ and values of $`\sigma _8`$ in the range 0.6–1, the uncertainty on $`\sigma _8`$ due to the presence of low mass systems with $`T<2`$ keV is around $`5\%`$, well within our estimated errors on the angular correlation function (Eq. 8 below). We therefore focus on the effects of changing the normalization, Eq. (7), of the self–similar SZ flux–Mass relation ($`\alpha =\gamma =0`$). Physically, this represents an uncertainty in the average thermal gas energy of galaxy clusters. Figure 4 shows the effect of changing $`Y_{}`$ at a flux limit of $`Y_{\mathrm{lim}}=10^4`$ arcmin2. From the figure, we see that changing this normalization stretches the counts and angular function contours vertically by roughly the same amount, their point of intersection (if any) remaining at roughly the same value of $`\mathrm{\Omega }_\mathrm{M}`$. For example, a factor 2 change in normalization moves the intersection by $`1020`$ % in $`\sigma _8`$ (e.g., at $`\mathrm{\Omega }_\mathrm{M}0.3`$). This represents an inherent systematic caused by modeling uncertainty associated with the ICM. Uncertainties related to ICM modeling have been extensively discussed by Moscardini et al. (2002) for the full 3D spatial correlation function of SZ clusters.
Including the uncertain SZ flux normalization, $`Y_{}`$, we now have three parameters to determine, and so additional information is needed; for example, an observational program to determine the $`Y`$$`M`$ relation on a representative sample of clusters. Another tactic is to use the constraint arising from the local cluster abundance as measured by the X–ray temperature function.

#### 4.3.1 Adding the local cluster abundance constraint

For a fixed value of $`T_{}`$, the local X–ray temperature function constrains the cosmological parameters $`\sigma _8`$ and $`\mathrm{\Omega }_\mathrm{M}`$ to a well–defined curve in the plane: $`\sigma _8\mathrm{\Omega }_\mathrm{M}^{0.6}0.6T_{}^{0.8}`$ (Pierpaoli et al. 2002). The exact relation depends somewhat on the chosen mass function. Consider such a constraint, shown as the middle dashed line in Figure 5. Supposing that we measure now both the SZ counts and angular function, finding 7.7 $`sq.deg`$ and $`w(\theta =30\mathrm{arcmins})=0.024`$, respectively, at $`Y=10^4`$ arcmin2 and given as the solid contours in the figure. Apart from the cosmological parameters $`\sigma _8`$ and $`\mathrm{\Omega }_\mathrm{M}`$, the SZ contours depend only on $`Y_{}`$, whose value is uncertain; by changing $`Y_{}`$, we move the SZ contours around the plane, as previously described. However, there is a unique value for $`Y_{}`$ that reduces the intersection of the three lines to a single point in the plane. By adjusting $`Y_{}`$ to this value, we simultaneously constrain the cosmological parameters and the normalization $`Y_{}`$. In the figure, this intersection lies on top of the true underlying model. The determined value of $`Y_{}`$ together with the known $`T_{}`$ imply a constraint on $`f_{\mathrm{gas}}`$555We are implicitly assuming that mean electron temperature, relevant to the SZ effect, and the observable X–ray temperature, parametrized by $`T_{}`$, are the same; this need not be the case.. By adding the information on the local cluster abundance, we have pinned down the relevant ICM physics and constrained the cosmological parameters, eliminating the primary uncertainty in using SZ clustering observations (Moscardini et al. 2002). Unfortunately, $`T_{}`$ is at present only poorly known, with simulations and observations indicating values generally in the range $`1.22`$ (Pierpaoli et al. 2002; Muanwong et al. 2002; Huterer $`\&`$ White 2002). We expect that it will be much better determined by the time large SZ catalogs become available, making this kind of analysis possible.

This example illustrates the utility of the angular function: with only the local cluster abundance and the SZ counts, we cannot determine the cosmological parameters $`\sigma _8`$ and $`\mathrm{\Omega }_\mathrm{M}`$, due to the uncertainty on $`Y_{}`$. Adding the angular function breaks the degeneracy. We are thus able to constrain the cosmology with the 2D SZ catalog, without recourse to redshift determinations. Constraints obtained in this way are particularly useful for their complementarity to constraints from CMB anisotropy and distance measurements with SNIa. For example, with SZ clusters we measure $`\sigma _8`$ directly on the relevant scales and constrain $`\mathrm{\Omega }_\mathrm{M}`$, as opposed to the physical density $`\mathrm{\Omega }_\mathrm{M}h^2`$ in the case of the CMB.

As mentioned above, clusters below $`2`$ keV contribute little to the contours at our fiducial flux limit of $`Y=10^4`$ arcmin2. At lower flux limits, more representative of future ground–based surveys, we could expect that the possible effects of deviations from self–similar scaling become more important, although simulations at present indicate only mild effects on the SZ flux.

### 4.4 Discussion of statistical errors

Before concluding, we briefly examine the statistical errors expected on a measurement of the SZ angular correlation function. Although a thorough analysis of this issue requires detailed simulations of a given survey, we may nonetheless make some general arguments to gain insight into what may eventually be achieved. We take the Planck survey as our example. The resolution of the Planck SZ catalog will be on the order of 5 arcmins (at best), so we can only expect to measure the correlations on larger scales (as the $`30\mathrm{}`$ we have chosen in our analysis), and the fiducial sensitivity expected is $`Y10^4`$ arcmin2 (Aghanim et al. 1997; Bartelmann 2001). Since the angular correlations are small in this context ($`0.1`$), we may estimate the (statistical) error on a measurement of $`w(\theta )`$ as the Poisson variance in the number of pairs, $`n_{\mathrm{pair}}`$, at this separation (e.g., Peebles 1980; Landy & Szalay 1993). This quantity is determined by the counts as follows: Suppose that we measure $`w`$ in a annular bin of width $`\mathrm{\Delta }\theta `$ at angular distance $`\theta `$ from a cluster. The mean number of clusters in this ring is $`n=2\pi \theta \mathrm{\Delta }\theta \mathrm{\Sigma }(Y_{\mathrm{lim}})`$, from which we deduce that the total number of pairs at this separation in a catalog of $`N`$ clusters is about $`n_{\mathrm{pair}}(1/2)N\times n=N\pi \theta \mathrm{\Delta }\theta \mathrm{\Sigma }`$. In other words, we estimate the error to be

$`\mathrm{\Delta }w(\theta ,Y_{\mathrm{lim}})`$ $`{\displaystyle \frac{1}{2\pi }}\left[\mathrm{\Sigma }(Y_{\mathrm{lim}})\right]^1{\displaystyle \frac{1}{\theta }}\left({\displaystyle \frac{\mathrm{\Delta }\theta }{\theta }}\right)^{1/2}`$ (8)
$`2.78\times 10^3\left({\displaystyle \frac{\mathrm{deg}}{\theta }}\right)\left({\displaystyle \frac{\mathrm{\Delta }\theta }{\theta }}\right)^{1/2}\left({\displaystyle \frac{\mathrm{\Sigma }}{\mathrm{deg}^2}}\right)^1`$
Notice that this statistical error depends essentially on the number counts $`\mathrm{\Sigma }`$. This leads to the gray shaded area in Figure 5

## 5 Discussion and Conclusion

We have calculated the angular correlation function of SZ–detected clusters in order to evaluate its usefulness for extracting cosmological information directly from a 2D SZ cluster survey, before 3D follow–up. The angular correlation function of SZ detected clusters differs from the angular power spectrum of SZ induced CMB anisotropies, which is dominated by the Poisson term. The different scaling of angular correlations with survey depth visually demonstrates the cosmological sensitivity of the angular function (See Figs 12). As illustration, we considered the parameter pair $`(\mathrm{\Omega }_\mathrm{M},\sigma _8)`$ in the context of flat $`\mathrm{\Lambda }`$–CDM models. We found that at sufficient depth (e.g., $`Y_{\mathrm{lim}}10^4`$ arcmin2, comparable to the Planck mission), the counts and angular function combined can constrain these parameters. Modeling uncertainty associated with the ICM may be reduced by measuring the $`T`$$`M`$ relation and adding the corresponding constraint from the local abundance of X–ray clusters; in this way, the two cosmological parameters and the SZ flux normalization $`Y_{}`$ may be found. Deeper ground–based surveys (e.g., $`Y_{\mathrm{lim}}10^5`$ arcmin2), will pick up a larger number of low–mass objects whose ICM properties require more careful modeling of deviations from self–similarity.

The accuracy and precision with which one will be able to measure the angular function (and the counts) is clearly an important issue. We only briefly touched on this point with a simple estimate of the expected statistical errors expected on $`w(\theta )`$ in the case of the Planck mission. A more detailed examination incorporating simulations of the SZ survey characteristics, including instrument noise and foreground contamination, as well as of the catalog extraction algorithms, is needed. The survey selection function will in reality depend on these details (Bartlett 2001; Melin et al. 2003; White 2003), and so will the errors on the counts and measured angular function.

It is clear that the most powerful constraints from an SZ survey will come from its measured redshift distribution. The interesting point is, however, that the method proposed here increases the immediate scientific return from an SZ survey by offering a way to obtain pertinent cosmological constraints using only a 2D SZ catalog, without recourse to follow–up observations.

###### Acknowledgements.

We are grateful to Lauro Moscardini, Antonaldo Diaferio and Antonio da Silva, and the anonymous referee for their very useful comments that helped the improvement of this paper. We would also like to thank Antonella De Luca and Milan Maksimovic, who first opened the contacts between us. S. Mei acknowledges support from the European Space Agency External Fellowship programme.


text:
# Turbulent magnetic Prandtl number and magnetic diffusivity quenching from simulations

## 1 Introduction

The concept of turbulent diffusion is often invoked when modeling large scale flows and magnetic fields in a turbulent medium. Turbulent magnetic diffusion is similar to turbulent thermal diffusion which characterizes the turbulent exchange of patches of warm and cold gas. This concept is also applied to turbulent magnetic diffusion which describes the turbulent exchange of patches of magnetic field with different strengths and direction. Reconnection of magnetic field lines is not explicitly required, but in the long run unavoidable if the magnetic power spectrum is to decrease toward small scales. The idea of Prandtl is that only the energy carrying eddies contribute to the mixing of large scale distributions of velocity and magnetic field structures. This leads to a turbulent magnetic diffusion coefficient $`\eta _\mathrm{t}\frac{1}{3}U\mathrm{}`$, where $`U`$ is the typical velocity and $`\mathrm{}`$ the scale of the energy carrying eddies. For the kinematic turbulent viscosity one expects similar values. Analytic theory based on the quasilinear approximation also produces similar (but not identical) values of $`\eta _\mathrm{t}`$ and $`\nu _\mathrm{t}`$ (e.g. Kitchatinov et al. 1994).

It is usually assumed that the values of $`\eta _\mathrm{t}`$ and $`\nu _\mathrm{t}`$ are independent of the molecular (microscopic) viscosity and magnetic diffusivity, $`\nu `$ and $`\eta `$. However, in the context of the geodynamo or in laboratory liquid metals the microscopic magnetic Prandtl number, $`P_\mathrm{m}=\nu /\eta `$ is very small ($`10^5`$). This has dramatic consequences for the magnetorotational instability (Balbus & Hawley 1991). This instability is generally accepted as the main mechanism producing turbulence in accretion discs (Balbus & Hawley 1998). For sufficiently small values of $`P_\mathrm{m}`$, however, this instability is suppressed (Rüdiger & Shalybkov 2002). On the other hand, the Reynolds number of the flow is quite large ($`10^5\mathrm{}10^6`$) and the flow therefore most certainly turbulent. This led Noguchi et al. (2002) to invoke a turbulent kinematic viscosity, $`\nu _\mathrm{t}`$, but to retain the microscopic value of $`\eta `$. The resulting effective magnetic Prandtl number they used was $`10^2`$ – big enough for the magnetorotational instability to develop. On may wonder, of course, why one should not instead use turbulent values for both coefficients, i.e. $`\nu _\mathrm{t}/\eta _\mathrm{t}1`$. This would lead to even more favorable conditions for the magnetorotational instability (Rüdiger et al. 2002).

Similar constraints have also been reported for the convection-driven geodynamo: Christensen et al. (1999) found that there is a minimum value of $`P_\mathrm{m}`$ of about 0.25 below which dynamo action does not occur at all. Similar results have also been reported by Cattaneo (2003). These results are disturbing, because both for the sun and for the earth, $`P_\mathrm{m}1`$. For $`P_\mathrm{m}`$ of order unity, on the other hand, earth-like magnetic configurations can more easily be reproduced (see Kutzner & Christensen 2002).

Because of these restrictions, one wonders whether the effective magnetic Prandtl number to be used is not $`P_\mathrm{m}`$, but rather the turbulent value, $`P_{\mathrm{m},\mathrm{t}}=\nu _\mathrm{t}/\eta _\mathrm{t}`$. This raises the important questions whether $`P_{\mathrm{m},\mathrm{t}}`$ is actually of order unity and whether it is independent of the microscopic value, $`P_\mathrm{m}`$. The aim of this paper is to estimate the value of $`P_{\mathrm{m},\mathrm{t}}`$ using turbulence simulations.

The knowledge of the value of $`P_{\mathrm{m},\mathrm{t}}`$ is also important for the solar dynamo. The qualitative properties of the dynamo depend on the relative importance of the large scale flows and hence on the magnitude of $`\eta _\mathrm{t}`$. If $`\eta _\mathrm{t}`$ is too large, the influence of a meridional flow of say 10 m/s is small so that only little modification can be expected for the basic $`\alpha \mathrm{\Omega }`$-dynamo (Roberts & Stix 1972). In this case, however, we know that conventional dynamo models of the solar activity cycle have difficulty to explain Spörer’s law of equatorward sunspot migration. The alternative that the resulting poleward migration can be overcompensated by an internal equatorward flow requires a sufficiently small value of $`\eta _\mathrm{t}`$, which implies that $`P_{\mathrm{m},\mathrm{t}}>1`$ (Choudhuri et al. 1995; Dikpati & Charbonneau 1999; Bonanno et al. 2002).

Given the importance of the value of the turbulent magnetic Prandtl number it is useful to assess the problem using three-dimensional simulations of turbulent flows. We determine $`\nu _\mathrm{t}`$ and $`\eta _\mathrm{t}`$ by measuring the decay rate of a large scale (mean) velocity and magnetic field, $`\overline{𝒖}`$ and $`\overline{B}`$, respectively. We emphasize that we are not addressing the question whether $`\nu _\mathrm{t}`$ and $`\eta _\mathrm{t}`$ can really be used in studies of the dynamo or the magnetorotational instability, for example.

We consider weakly compressible nonhelically forced turbulence and use a model similar to that of Brandenburg (2001), but with kinetic helicity fluctuating about zero. Dynamo action for such a model has recently been considered by Haugen et al. (2003), but it sets in only at magnetic Reynolds numbers above $`30`$, which is not the case in the present simulations. We begin however by first reviewing the basic results for the values of $`\nu _\mathrm{t}`$ and $`\eta _\mathrm{t}`$ within the framework of the quasilinear (Roberts & Soward 1975, Rüdiger 1989) and other approximations.

## 2 Results from quasilinear approximation

For steady homogeneous isotropic turbulence the correlation tensor is independent of $`x`$ and $`t`$, i.e.

$$u_i^{}(x,t)u_j^{}(x+\xi ,t+\tau )=Q_{ij}(\xi ,\tau ),$$
(1)
where angular brackets denote an ensemble average and primes fluctuations about the average. In the quasilinear approximation the transport coefficients are conveniently expressed in terms of the Fourier transformed correlation tensor, $`\widehat{Q}_{ij}(k,\omega )`$, which is normalized such that

$$Q_{ij}(\xi ,\tau )=\widehat{Q}_{ij}(k,\omega )e^{i(k\xi \omega \tau )}𝑑k𝑑\omega .$$
(2)
For the turbulent viscosity and the turbulent magnetic diffusivity one finds respectively (Rüdiger 1989)

$$\nu _\mathrm{t}=\frac{4}{15}\frac{\nu ^3k^6\widehat{Q}_{ll}(k,\omega )}{(\omega ^2+\nu ^2k^4)^2}𝑑k𝑑\omega ,$$
(3)
$$\eta _\mathrm{t}=\frac{1}{3}\frac{\eta k^2\widehat{Q}_{ll}(k,\omega )}{\omega ^2+\eta ^2k^4}𝑑k𝑑\omega .$$
(4)
Obviously, both quantities are of the same order of magnitude, but they are not identical. In the limits $`\nu ,\eta 0`$ the expressions are drastically simplified, i.e.

$$\nu _\mathrm{t}=\frac{1}{15}\underset{\mathrm{}}{\overset{\mathrm{}}{}}u^{}(x,t)u^{}(x,t+\tau )𝑑\tau $$
(5)
and

$$\eta _\mathrm{t}=\frac{1}{6}\underset{\mathrm{}}{\overset{\mathrm{}}{}}u^{}(x,t)u^{}(x,t+\tau )𝑑\tau ,$$
(6)
so that for the turbulent magnetic Prandtl number is

$$P_{\mathrm{m},\mathrm{t}}=\frac{\nu _\mathrm{t}}{\eta _\mathrm{t}}=\frac{2}{5}=0.4.$$
(7)
This results is similar to that of Nakano et al. (1979) for the thermal Prandtl number.

Rüdiger (1989) lists a number of other approaches for calculating turbulent transport coefficients, which all yield Prandtl numbers around or below unity. One particular approach is the renormalization group analysis which was applied to turbulence by Forster et al. (1977) for the case of a passive scalar, and later by Fournier et al. (1982) to the case with magnetic fields. These results are valid in the long-time large-scale limit, and the value of $`P_{\mathrm{m},\mathrm{t}}`$ turned out to be close to 0.7; see Eq. (23) of Fournier et al. (1982).

Kitchatinov et al. (1994) use a mixing length approximation where terms of the form $`\mathrm{d}/\mathrm{d}t\nu ^2`$ are replaced by $`\tau _{\mathrm{corr}}^1`$, where $`\tau _{\mathrm{corr}}`$ is the correlation time of the turbulence. They find $`\nu _\mathrm{t}=(4/15)\tau _{\mathrm{corr}}u_{\mathrm{rms}}^2`$ and $`\eta _\mathrm{t}=(1/3)\tau _{\mathrm{corr}}u_{\mathrm{rms}}^2`$, so $`P_{\mathrm{m},\mathrm{t}}=4/5=0.8`$. Yet another approach is the $`\tau `$-approximation where triple correlations are replaced by a damping term that is proportional to the quadratic moments (e.g. Kleeorin et al. 1996, Blackman & Field 2002). Here no Fourier transformation in time is used. This gives, as before, $`\eta _\mathrm{t}=(1/3)\tau u_{\mathrm{rms}}^2`$ (where $`\tau `$ is now interpreted as a relaxation time), but $`\nu _\mathrm{t}=(2/15)\tau u_{\mathrm{rms}}^2`$, so $`P_{\mathrm{m},\mathrm{t}}=2/5=0.4`$. This is half the value obtained from the mixing length approximation, but the same as in Eq. (7).

The fact that in all these cases $`P_{\mathrm{m},\mathrm{t}}`$ is less than unity can be traced back to the presence of the pressure term in the momentum equation. If this term is ignored (as in pressureless Burgers turbulence or ‘burgulence’) one always gets $`P_{\mathrm{m},\mathrm{t}}=1`$.

It is tempting to speculate that the discrepancy between the different analytic approaches is related to the validity of some idealizing assumptions made in order to apply the quasilinear and other approximations. Clearly, additional approaches are needed to get a more complete picture regarding the correct value of $`P_{\mathrm{m},\mathrm{t}}`$. It is nevertheless encouraging that $`P_{\mathrm{m},\mathrm{t}}`$ does not strongly deviate from unity.

In the remainder of this paper we estimate $`\nu _\mathrm{t}`$ and $`\eta _\mathrm{t}`$ numerically by considering the decay of an initial large scale velocity or magnetic field, respectively, in the presence of small scale turbulence.

## 3 The model

The equations describing compressible isothermal hydromagnetic flows with constant sound speed, $`c_\mathrm{s}`$, are

$$\frac{\mathrm{D}u}{\mathrm{D}t}=c_\mathrm{s}^2\mathrm{ln}\rho +\frac{J\times B}{\rho }+F_{\mathrm{visc}}+𝒇,$$
(8)
$$\frac{\mathrm{D}\mathrm{ln}\rho }{\mathrm{D}t}=u,$$
(9)
$$\frac{B}{t}=\times (u\times B)+\eta ^2B,$$
(10)
where $`u`$ is the velocity, $`\rho `$ the density, $`B`$ is the magnetic field, and $`J=\times B/\mu _0`$ is the current density with $`\mu _0`$ being the vacuum permeability. The viscous force is

$$F_{\mathrm{visc}}=\nu \left(^2u+\frac{1}{3}u+2𝗦\mathrm{ln}\rho \right),$$
(11)
where $`𝖲_{ij}=\frac{1}{2}(u_{i,j}+u_{j,i})\frac{1}{3}\delta _{ij}u`$ is the traceless rate of strain tensor.

We solve the equations using the Pencil Code111 http://www.nordita.dk/data/brandenb/pencil-code, which is a memory-efficient sixth-order finite difference code using the $`2N`$-RK3 scheme of Williamson (1980). For most of the simulations a resolution of $`128^3`$ meshpoints is used, but in Sect. 5 a higher resolution of up to $`512^3`$ meshpoints was necessary.

We focus on the case where the forcing, $`𝒇`$, occurs at a wavenumber around $`k_\mathrm{f}=10`$. The forcing is such that the turbulence is subsonic and nonhelical. We consider two different periodic initial conditions,

$$B=(\mathrm{cos}k_1z,0,0)B_0\text{(nonhelical)}$$
(12)
and

$$B=(\mathrm{cos}k_1z,\mathrm{sin}k_1z,0)B_0\text{(helical)},$$
(13)
where $`B_0`$ is the amplitude of the initial field. In the fully helical case one may expect a different decay time because the magnetic helicity is a conserved quantity in the limit of small magnetic diffusivity. For the velocity field we use similar initial conditions, but we do not expect this to be sensitive to helicity, because kinetic helicity is not conserved in the limit $`\nu 0`$, and would only be conserved in the unphysical case $`\nu =0`$.

A detailed discussion of the initial conditions may at first glance appear somewhat surprising, because for forced turbulent flows the initial conditions are normally forgotten after about one turnover time. This is indeed the case for hydrodynamic turbulence, but not for hydromagnetic turbulence if the magnetic field has net magnetic helicity. The reason is that, regardless of the level of turbulence, the net magnetic helicity can only change on the resistive time scale. Our results below confirm this and they are indeed in agreement with earlier model predictions (cf. Blackman & Brandenburg 2002). The situation would be different if the initial field was bi-helical, i.e. with oppositely helical contributions at different scales. This case has been studied elsewhere (Yousef & Brandenburg 2003).

In Fig. 1 we show kinetic and magnetic energy spectra of the run with $`\text{Re}=150`$ and $`R_\mathrm{m}=15`$ at three different times using a resolution of $`512^3`$ meshpoints. The kinetic energy shows indications of a short inertial range in $`15<k<40`$. Below the forcing scale, in $`2<k<9`$, velocity and magnetic fields are random and $`\delta `$-correlated in space, giving rise to a $`k^2`$ spectrum. The magnetic energy is substantially weaker than the kinetic energy. This is because here the magnetic Prandtl number is small, $`P_\mathrm{m}=0.1`$, and the magnetic Reynolds number is subcritical for dynamo action. With our definition of $`R_\mathrm{m}`$ the critical value lies around 25 (Haugen et al. 2003). The small scale magnetic energy is therefore maintained by constantly stirring the slowly decaying large scale field.

Given that the initial large scale field depends only on $`z`$, it makes sense to define a mean field by averaging over the $`x`$ and $`y`$ directions. Alternatively, one might define an average by Fourier filtering, but this has the disadvantage that not all the Reynolds rules are satisfied. For example, the average of a product of a mean and a fluctuating quantity would not vanish. However, for all practical purposes our horizontal average is nearly equivalent to a projection onto the $`k=k_1`$ Fourier mode. Indeed, the main reason for forcing at a large wavenumber, $`k_\mathrm{f}=10`$, is that we need some degree of scale separation. Without scale separation, there would be no way of distinguishing between mean and fluctuating fields. Since the velocity fluctuations are constantly driven via the forcing term, it would be impossible to measure any decay of the mean velocity. Nevertheless, even with scale separation there will always be a certain level of noise in the mean field whose energy is $`(k_1/k_\mathrm{f})^2`$ times smaller than energy of the fluctuations. This means that we can measure an exponential decay of the mean field only in a certain window where nonlinear effects are already weak, but were the noise level is not yet reached.

## 4 Results

### 4.1 Decay of $`\overline{u}`$ and $`\overline{B}`$

We begin by considering the decay of a helical large scale magnetic field and compare it with the decay of a large scale helical velocity field in a purely hydrodynamic simulation; see Fig. 2. Here, large scale velocity and magnetic fields are defined as horizontal averages over $`x`$ and $`y`$; the result is denoted by $`\overline{𝒖}`$ and $`\overline{B}`$, respectively. During the time interval when mean velocity and magnetic field decay exponentially, the corresponding decay rates are determined as

$$\lambda _u(\overline{𝒖})=\frac{\mathrm{d}\mathrm{ln}\overline{𝒖}{}_{}{}^{2}^{1/2}}{\mathrm{d}t},\lambda _B(\overline{B})=\frac{\mathrm{d}\mathrm{ln}\overline{B}^2^{1/2}}{\mathrm{d}t}.$$
(14)
In the graphs of $`\lambda _u(\overline{u})`$ and $`\lambda _B(\overline{B})`$ an exponential decay shows up as a plateau. The magnetic field decay is initially slow, so $`\lambda _B(\overline{B})`$ is initially not constant, but then it speeds up and $`\lambda _B(\overline{B})`$ reaches a plateau. The decay of the velocity field is immediately fast and $`\lambda _u(\overline{u})`$ lies immediately on a plateau. This suggests that the turbulent magnetic diffusivity is affected by the strong initial field that in turn gives rise to a quenching of the turbulent magnetic diffusivity. Strong means that the magnetic field strength is comparable with the equipartition field strength, $`B_{\mathrm{eq}}=\mu _0\rho u^2^{1/2}`$. The initially strong large scale flow and the associated vorticity, on the other hand, do not and are also not expected to affect the turbulent viscosity and the associated decay of this large scale flow. For $`|\overline{B}|B_{\mathrm{eq}}`$, however, both $`\overline{𝒖}`$ and $`\overline{B}`$ decay at the same rates, $`\lambda _u`$ and $`\lambda _B`$, respectively. This allows us to calculate

$$\nu _\mathrm{t}=\lambda _u/k_1^2,\eta _\mathrm{t}=\lambda _B/k_1^2,$$
(15)
where $`k_1`$ is the wavenumber of the initial large scale velocity and magnetic fields. From the present simulations, where $`k_\mathrm{f}/k_1=10`$, we find

$$\nu _\mathrm{t}\eta _\mathrm{t}=(0.8\mathrm{}0.9)\times u_{\mathrm{rms}}/k_\mathrm{f}\text{(for }\overline{B}^2B_{\mathrm{eq}}^2\text{)}.$$
(16)
Once $`|\overline{𝒖}|`$ has decreased below a certain level ($`<0.1u_{\mathrm{rms}}`$), it cannot decay further and continues to fluctuate around $`0.08u_{\mathrm{rms}}`$, corresponding to the level of the rms velocity of the (forced!) turbulence at $`k=k_1`$ (see the dashed line in Fig. 2).

The quenching of the magnetic diffusivity, $`\eta _\mathrm{t}=\eta _\mathrm{t}(\overline{B})`$, can be obtained from one and the same run by simply determining the decay rate, $`\lambda _B(\overline{B})`$, at different times, corresponding to different values of $`\overline{B}=|\overline{B}|`$; see Fig. 3. To describe departures from purely exponential decay we adopt a $`\overline{B}`$-dependent $`\eta _\mathrm{t}`$ expression of the form

$$\eta _\mathrm{t}(\overline{B})=\eta _{\mathrm{t0}}/(1+a\overline{B}^2/B_{\mathrm{eq}}^2),$$
(17)
where $`\eta _{\mathrm{t0}}`$ is the unquenched (kinematic) value of $`\eta _\mathrm{t}`$, described approximately by Eq. (16), and $`a`$ is a fit parameter. According to Cattaneo & Vainshtein (1991) the parameter $`a`$ is expected to be of the order of the magnetic Reynolds number based on the microscopic magnetic diffusivity,

$$R_\mathrm{m}=u_{\mathrm{rms}}k_\mathrm{f}/\eta .$$
(18)
Figure 3 suggests that $`a0.4R_\mathrm{m}`$.

Before we discuss the effective quenching behavior of $`\eta _\mathrm{t}`$ in more detail we should note that Eq. (17), and in particular the value of $`a`$, do not apply universally and depend on the field geometry. This is easily demonstrated by considering a nonhelical initial field. In that case the decay becomes unquenched already for $`\overline{B}^2/B_{\mathrm{eq}}^21`$. Equation (17) can still be used as a reasonable fit formula, but now $`a=1`$ produces a good fit (independent of $`R_m`$); see Fig. 4.

In the nonhelical case there is an initial phase where the field increases due to the wind-up of the large scale field. Since we measure $`\eta _\mathrm{t}`$ from the decay rate of the large scale field, this would formally imply negative values of $`\eta _\mathrm{t}`$. Traces of this effect can still be seen in Fig. 4 near $`\overline{B}^2/B_{\mathrm{eq}}^2=1`$. For this reason our method can only give reliable results if $`|\overline{B}|\begin{array}{c}<\\ \end{array}0.8B_{\mathrm{eq}}`$. In the case of a helical initial field, on the other hand, we have $`\overline{J}\times \overline{B}=0`$, i.e. the large scale field is force-free and interacts only weakly with the turbulence. In particular, there is no significant amplification from the initial wind-up of the large scale field.

### 4.2 Comparison with the dynamical quenching model

In the case of a helical field and for $`\overline{B}^2/B_{\mathrm{eq}}^2\begin{array}{c}>\\ \end{array}R_\mathrm{m}^1`$ the slow decay of $`\overline{B}`$ is related to the conservation of magnetic helicity. As discussed already by Blackman & Brandenburg (2002), this behavior is related to the phenomenon of selective decay (e.g. Montgomery et al. 1978) and can be described by the dynamical quenching model. This model goes back to an early paper by Kleeorin & Ruzmaikin (1982, see also Kleeorin et al. 1995), but it applies even to the case where the turbulence is nonhelical and where there is no $`\alpha `$ effect in the usual sense. However, the magnetic contribution to $`\alpha `$ is still non-vanishing because it is driven by the helicity of the large scale field.

To demonstrate this quantitatively we solve, in the one mode approximation ($`k=k_1`$) with $`\overline{B}=\widehat{B}\mathrm{exp}(\mathrm{i}k_1z)`$, the mean-field induction equation

$$\frac{\mathrm{d}\widehat{B}}{\mathrm{d}t}=\mathrm{i}k_1\times \widehat{𝓔}\eta k_1^2\widehat{B}$$
(19)
together with the dynamical $`\alpha `$-quenching formula \[Eq. (13) of Blackman & Brandenburg (2002)\]

$$\frac{\mathrm{d}\alpha }{\mathrm{d}t}=2\eta k_\mathrm{f}^2\left(\alpha +\stackrel{~}{R}_\mathrm{m}\frac{\mathrm{Re}(\widehat{𝓔}^{}\widehat{B})}{B_{\mathrm{eq}}^2}\right),$$
(20)
where

$$\widehat{𝓔}=\alpha \widehat{B}\eta _\mathrm{t}\mathrm{i}k_1\times \widehat{B}$$
(21)
is the electromotive force, and $`\stackrel{~}{R}_\mathrm{m}`$ is defined as the ratio $`\eta _{\mathrm{t0}}/\eta `$, which is expected to be close to the value of $`R_\mathrm{m}`$ as defined by Eq. (18).

In Fig. 5 we show the evolution of $`\overline{B}/B_{\mathrm{eq}}`$ for helical and nonhelical initial conditions, $`\widehat{B}(1,\mathrm{i},0)`$ and $`\widehat{B}(1,0,0)`$, respectively. In the case of a nonhelical field, the decay rate is not quenched at all, but in the helical case quenching sets in for $`\overline{B}^2/B_{\mathrm{eq}}^2\begin{array}{c}>\\ \end{array}R_\mathrm{m}^1`$.

In the helical case, the onset of quenching at $`\overline{B}^2/B_{\mathrm{eq}}^2R_\mathrm{m}^1`$ is well reproduced by the simulations. In the nonhelical case, however, some weaker form of quenching sets in when $`\overline{B}^2/B_{\mathrm{eq}}^21`$ (Fig. 4). We refer to this as standard quenching (e.g. Kitchatinov et al. 1994) which is known to be always present. In Blackman & Brandenburg (2002) this was modeled by allowing in Eq. (21) $`\eta _\mathrm{t}`$ to be $`\overline{B}`$-dependent. They adopted the formula

$$\eta _\mathrm{t}=\eta _{\mathrm{t0}}/(1+\stackrel{~}{g}|\overline{B}|/B_{\mathrm{eq}})$$
(22)
and found that, for a range of different values of $`R_\mathrm{m}`$, $`\stackrel{~}{g}=3`$ resulted in a good description of the simulations of cyclic $`\alpha \mathrm{\Omega }`$-type dynamos (Brandenburg et al. 2002). We emphasize that this $`\eta _\mathrm{t}`$ is not used in a diagnostic way as in Eq. (17), but rather in the numerical solution of Eqs (19) and (20). The resulting decay law, shown as a dotted line in Fig. 5, agrees now with the decay law seen in the turbulence simulations (Fig. 2). The helical case with $`\stackrel{~}{g}=3`$ is still compatible with the simulations.

## 5 Independence of microscopic viscosity

Finally we need to show that the turbulent magnetic Prandtl number is indeed independent of the microscopic magnetic Prandtl number. In Fig. 6 we plot the decay rates, obtained by differentiating $`\mathrm{ln}\overline{B}(t)`$, for three different values of the microscopic viscosity, keeping $`\eta `$ fixed. The resulting values of the flow Reynolds number, $`\text{Re}=u_{\mathrm{rms}}k_\mathrm{f}/\nu `$, vary between 20 and 150, giving $`P_\mathrm{m}`$ in the range between 0.1 and 1. Within plot accuracy the three values of $`\lambda _B`$ turn out to be identical in the interval where the decay is exponential.

The duration of this interval is $`u_{\mathrm{rms}}k_\mathrm{f}\mathrm{\Delta }t200`$, which is comparable to the time interval in Fig. 2 during which the decay is exponential. In one of the three cases ($`R_\mathrm{m}=20`$) the initial field was rather strong, so that it takes a long time before the magnetic helicity constraint becomes unimportant so that the field can decay exponentially ($`u_{\mathrm{rms}}k_\mathrm{f}t800`$).

The numerical resolution used in most of the models is $`128^3`$ mesh points. However, as Re is increased, higher resolution is required. For Re=80 we used $`256^3`$ mesh points and for Re=150 we used $`512^3`$ mesh points. This implies mesh a Reynolds number, $`u_{\mathrm{rms}}\mathrm{\Delta }x/\nu `$, based on the mesh spacing $`\mathrm{\Delta }x`$, of about 18. Empirically we know that larger values are not generally possible.

## 6 Conclusions

The turbulence simulations presented here have shown that the turbulent magnetic Prandtl number is always of order unity, regardless of the values of the microscopic magnetic Prandtl number. Under the assumption of incompressibility, both the quasilinear approximation and the renormalization group approach give turbulent magnetic Prandtl numbers somewhat below unity, which is related to the pressure term in the momentum equation. Here we find instead $`P_{\mathrm{m},\mathrm{t}}1`$. There are several plausible reasons for this discrepancy: (i) our simulations are actually weakly compressible, (ii) they are non-steady and, (iii) the idealizing assumptions made in the analytic approaches may not be justified.

Our results have also shown that, for nonhelical magnetic fields, the turbulent magnetic diffusivity is quenched when the magnetic energy becomes comparable to the kinetic energy. For helical magnetic fields, however, an apparent suppression of the decay rate is observed which agrees with predictions from a dynamical quenching model. If this suppression is described by an algebraic expression, quenching would set in for magnetic energies much below the kinetic energy.

The present work demonstrates that the dynamical quenching approach is not restricted to dynamos, but it can also deal with decay problems, as was already mentioned in Blackman & Brandenburg (2002). The dynamical quenching model is usually formulated in terms of $`\alpha `$, but for helical mean fields $`\overline{J}`$ and $`\overline{B}`$ are parallel and the separation into contributions from $`\alpha \overline{B}`$ and $`\eta _\mathrm{t}\overline{J}`$ becomes less meaningful. It is for this reasons that an $`\alpha `$ term appears in the description of the decay of helical fields, rather than a dynamical contribution to $`\eta _\mathrm{t}`$-quenching.

The remaining quenching of $`\eta _\mathrm{t}`$ that affects both helical and nonhelical fields is consistent with an algebraic quenching formula that is non-catastrophic, i.e. independent of the microscopic magnetic diffusivity.

Although our results suggest that the turbulent magnetic Prandtl number is of order unity, we cannot claim that it is safe to use turbulent viscosity and magnetic diffusivity in a simulation of the dynamo or the magnetorotational instability, for example, as a replacement of a fully resolved simulation. First of all, the functional form of the turbulent transport coefficients is for realistic turbulent flows more complicated and involves in practice tensorial rather than scalar coefficients. Numerical evidence for this has been presented elsewhere in the context of shear flow turbulence (Brandenburg & Sokoloff 2002). Furthermore, there will be additional terms such as the $`\alpha `$-effect (see Sect. 4.2) and the AKA-effect (Frisch et al. 1987; see also Brandenburg & Rekowski 2001). Most importantly, turbulent transport may be nonlocal, as is well known in meteorology when modeling atmospheric flows (Stull 1984, Ebert et al. 1989), where the turbulent transport is described by so-called transilient matrices (see also Miesch et al. 2000 for examples of astrophysical convection). Nonlocal transport means that the transport coefficients have to be replaced by integral kernels. In Fourier space, the convolution with an integral kernel corresponds to a multiplication with a wavenumber dependent factor. There is indeed some evidence that the main contribution comes only from the smallest wavenumbers (Brandenburg & Sokoloff 2002). This is primarily a consequence of a lack of scale separation in naturally forced turbulence, such as shear flows or convection. In the present context, however, this is not an issue because we have deliberately considered the case where the scale of the turbulent eddies is much smaller than the scale of the large scale field ($`k_\mathrm{f}/k_1=10`$).

Finally, we wish to point out that studies of instabilities (e.g. the magnetorotational or the dynamo instability) using turbulent transport coefficients can sometimes lead to paradoxical situations. In the case of solar convection, for example, one expects from mixing length theory that turbulent viscosity and thermal diffusivity are on the order of a few times $`10^{12}\mathrm{cm}^2\mathrm{s}^1`$. However, using such values in a global model of the sun leads to an instability (Rüdiger 1989, Rüdiger & Spahn 1992), which is in fact nothing but a repetition of the original convection instability that leads to turbulence in the first place (Tuominen et al. 1994). It is therefore plausible that the actual values of the turbulent transport coefficients should rather be close to the those for marginal stability. This would lead to a global constraint similar to the magnetic helicity constraint that governs the nonlinear behavior of the $`\alpha `$-effect in helical hydromagnetic turbulence. At present, however, there is no theoretical framework that allows self-consistent modeling of convection using mean-field theory.

###### Acknowledgements.

We thank an anonymous referee for making useful suggestions and drawing our attention to the paper by Fournier et al. (1982). Use of the supercomputers in Odense (Horseshoe), Trondheim (Gridur), and Leicester (Ukaff) is acknowledged.
```
$Header: /home/brandenb/CVS/tex/mhd/turb_prandtl/paper.tex,v 1.51 2004/02/02 17:49:53 tarek Exp $

```