text:
# Normal parameters for an analytic description of the CMB cosmological parameter likelihood

## 1 Introduction

The challenge of turning a CMB dataset into cosmological parameter constraints is one that has been solved by a series of data compressions. The information in the time–ordered data is compressed into a map (Wright et al., 1996; Tegmark, 1997a). The information in the map is then compressed into a power spectrum (Tegmark, 1997b; Bond et al., 1998, 2000; Wandelt et al., 2001; Bartlett et al., 2000; Hivon et al., 2002). Finally, the information in the power spectrum is then compressed into cosmological parameters (Lineweaver et al., 1997; Benoit, 2002).

This last step, however, is more of a data “explosion” than a data compression. Although the number of parameters is indeed small ($`10`$) the non–Gaussian distribution of their errors may be characterized by $`10^{10}`$ numbers for grid–based likelihood evaluation (Tegmark & Zaldarriaga, 2000) or perhaps as little as $`30,000`$ for the Monte-Carlo Markov Chain method (Christensen et al., 2001).

This explosion has undesirable consequences. The full constraining power of the data is not used by the cosmological community. Typically authors report not the full (cumbersome) likelihood, but projections or marginalizations of it down to one or at most two–dimensional spaces. Given the near–degeneracies that exist in the probability distribution, and its non-Gaussian nature, these final steps lose significant information.

Our work here provides the final step of compression in the data analysis pipeline. We compress the probability distribution of the cosmological parameters down to the $`50`$ parameters of an analytic form for the likelihood (a mean and covariance matrix). This step makes the full information in the parameter likelihood function easily useable.

Our procedure is analogous to the ‘radical compression’(Bond et al., 2000) used for compressing the uncertainty in the power spectrum estimates, whose distribtuion is also non–Gaussian. In both cases a non–linear variable transformation is used, with the property that the transformed variables are well–approximated by a normal distribution.

We were inspired to search for a set of normal parameters by a proposal that the $`C_l`$ can have a linear dependence on a set of well-chosen parameters that span the space of possible models (Kosowsky et al., 2002). Clearly such a parameter set would be tremendously valuable for parameter estimation since the $`C_l`$ for any given model could be calculated practically instantaneously simply by summing the terms in the first order Taylor expansion. Such a set would also, given sufficiently Gaussian errors on $`C_l`$, have errors with a normal distribution.

Although inspired by Kosowsky et al. (2002), our goal here is entirely different. It is to find a set of parameters (the normal parameters) for which the likelihood is Gaussian. Although the approximate linearity demonstrated by Kosowsky et al. (2002) gave us confidence to try to find a set of normal parameters, it by no means guaranteed success. The linear-response approximation has not yet been demonstrated to be sufficiently accurate for purposes of parameter estimation given MAP data (and certainly not given current data) and the errors on $`C_l`$ are not actually Gaussian.

While Kosowsky et al. (2002) are developing a tool for parameter estimation from CMB data, our work is aimed at being able to easily take full advantage of such parameter estimation. The most straightforward application of CMB data compressed in this form is for simultaneous analysis of CMB constraints with those from other cosmological probes. In the past to do this analysis one would have to entirely reproduce the reduction from bandpowers (plus their window functions, offset log–normal parameters, Fisher matrices, and calibration uncertainties) to cosmological parameter likelihoods. This procedure requires the organization of a lot of data as well as hundreds of thousands of angular power spectrum calculations. In contrast, with the compression to normal parameters in hand one need only perform a variable transformation and evaluate a $`10`$–dimensional Gaussian. We use a similar set of parameters as KMJ, and demonstrate explicitly the validity of the normal approximation by numerical computation of the likelihood given current CMB data.

One can also simplify a probability distribution by choosing linear combinations of the parameters that diagonalize the covariance matrix, the so–called parameter eigenmodes (Efstathiou & Bond, 1999). However, the resulting parameters are still dependent and non–Gaussion if the distribution of the original parameters is non-Gaussian. We do find diagonalization of the covariance matrix of the normal parameters to be useful, both for our fitting procedure and for gaining an understanding of the nature of the constraints the data places on the parameter space.

We test our procedure on current data and provide the mean and covariance matrix to the normal parameters that best fit the likelihood given current data. As we go to press, this fit is now out of date due to the data from the Wilkinson Microwave Anisotropy Probe (WMAP) 111http://map.gsfc.nasa.gov. However, we expect these data to further improve the validity of the normal approximation to the likelihood and thus our method will be of lasting value.

In section II we describe the normal parameters. In section III we describe our calculation of the exact likelihood via the Monte Carlo Markov Chain method (Christensen et al., 2001). In section IV we describe our procedure for finding the mean and covariance matrix that provide the best fit to our MCMC–calculated likelihood. In section V we present our results. In section VI we discuss applications and conclude.

## 2 The normal parameters

The normal parameters are a non-linear combination of the six cosmological parameters that we consider here – $`\omega _b`$, $`\omega _d`$, $`\mathrm{\Omega }_\mathrm{\Lambda }`$, $`z`$, $`A(2.73\times 10^6)^2\mathrm{Mpc}^3P_\mathrm{\Phi }(0.05\mathrm{Mpc}^1)`$ and $`n_s`$ (baryon density, dark matter density, dark energy density, reionization redshift, amplitude and scalar tilt of primordial power spectrum respectively). In particular, we replace $`\mathrm{\Omega }_\mathrm{\Lambda }`$, $`A`$ and $`n_s`$ by

$`\mathrm{\Theta }_s`$ $``$ $`\mathrm{\Theta }_s^Er_s/d_A`$ (1)
$`A^{}`$ $``$ $`{\displaystyle \frac{A}{76,000}}\left({\displaystyle \frac{0.05\mathrm{Mpc}^1}{k_{\mathrm{pivot}}}}\right)^{1n_s}e^{2\tau }`$ (2)
$`t`$ $``$ $`{\displaystyle \frac{1}{\sqrt{w_b}}}\mathrm{\hspace{0.17em}2}^{n_s1}`$ (3)
where $`r_s`$ is the sound horizon at recombination, $`d_A`$ is the angular diameter distance to the recombination surface and we take $`k_{\mathrm{pivot}}=0.067\mathrm{Mpc}^1`$. The distinction between $`\mathrm{\Theta }_s`$ and $`\mathrm{\Theta }_s^E`$ is explained below. The numerical factors are chosen so that $`A^{}`$ is of order unity. We use $`\tau `$ to denote the optical depth to Thomson scattering from here back to some time before reionization and after recombination. The resulting variable set, $`\{\omega _b,\omega _d,\mathrm{\Theta }_s,A^{},t,z\}`$, has a probability distribution, given CMB data, that is well–approximated by a Gaussian as we will demonstrate below.

We chose the normal parameters by considering what combinations of cosmological parameters affect the features in the angular power spectrum over the range in which the data are highly constraining ($`100l500`$). The parameter $`A^{}`$ describes the overall amplitude, $`t`$ is chosen to correlate with the ratio of the amplitude of the second peak to the first peak (at fixed $`\omega _m`$), and $`\mathrm{\Theta }_s`$ scales the features horizontally.

Hu et al. (2001) have investigated this phenomenology as well. For $`\mathrm{\Theta }_S`$ we use their approximation to the angular size of the sound horizon, which has the advantage of being an algebraic expression222We set $`\mathrm{\Theta }_s\pi /l_A`$ where $`l_A`$ is given by their equations A3-A5., rather than the exact sound horizon, $`\mathrm{\Theta }_s^E`$. This distinction is important because they can differ by more than the uncertainty in $`\mathrm{\Theta }_S`$. Hu et al. (2001) also define a parameter $`H_2`$ analogous to our $`t`$ but with a different parameter dependence.

The importance of $`\mathrm{\Theta }_s`$ for understanding the $`C_l`$ is widely recognized (Efstathiou & Bond, 1999; Tegmark & Zaldarriaga, 2000; Hu et al., 2001; Kaplinghat et al., 2002; Kosowsky et al., 2002). It is one of the best–determined cosmological quantities: $`\mathrm{\Theta }_s=0.^{}59\pm 0.^{}01`$ (Knox et al., 2001).

The parameter $`A^{}`$ sets the overall amplitude of the spectrum at $`l100`$. At these $`\mathrm{}`$ values, Thomson scattering depresses the amplitude by $`\mathrm{exp}(2\tau )`$. The value of $`k_{\mathrm{pivot}}`$ was chosen to decorrelate $`A^{}`$ and $`n_S`$. In an idealized case with uniform relative error bars on $`l(l+1)C_l/(2\pi )`$ from some $`l_{\mathrm{min}}`$ to $`l_{\mathrm{max}}`$ we would expect $`k_{\mathrm{pivot}}=\sqrt{l_{\mathrm{min}}l_{\mathrm{max}}}/\eta _0`$ where $`\eta _0=(14\pm 0.6)`$ Gpc is approximately the coordinate distance to the last–scattering surface. Taking $`l_{\mathrm{min}}=10`$ and $`l_{\mathrm{max}}=1000`$ we expect $`k_{\mathrm{pivot}}0.007\mathrm{Mpc}^1`$. At fixed $`\tau `$ (or $`z`$) we find this value works well. Allowing $`\tau `$ to vary leads to a correlation with $`n_S`$. A result of this correlation is that decorrelation between $`A^{}`$ and $`n_S`$ is best done with a much larger $`k_{\mathrm{pivot}}`$ of about $`0.067\mathrm{Mpc}^1`$. The $`\tau `$-$`n_S`$ correlation arises because both rely considerably on $`C_l`$ measurements at $`l<100`$.

The least “normal” of the normal parameters is $`z`$, which is not very well–constrained by CMB data. The situation may improve quite soon with large–angle polarization data from MAP (Kaplinghat et al., 2002). Or, it is already not a problem if one interprets observation of a Gunn–Peterson trough in a $`z=6.3`$ quasar as indicating $`z_{\mathrm{RI}}6.3`$ (Becker et al., 2001; Fan et al., 2002). We caution that though the data indicate the amount of neutral Hydrogen in the inter–galactic medium is rising from zero, how much the fraction of free electrons is increasing with $`z`$ is highly unconstrained (Kaplinghat et al., 2002).

We assume a step function transition of the ionization fraction from 0 to 1 at redshift $`z`$. For the flat models we consider this results in an optical depth of

$$\tau =0.038\omega _bh/\omega _m\left[\sqrt{\mathrm{\Omega }_\mathrm{\Lambda }+\mathrm{\Omega }_m(1+z)^3}1\right]$$
(4)
(Hu & White, 1997; Kaplinghat et al., 2002). We settled on $`z`$ as a normal parameter after trial and error. We have found $`z`$ to be more normal than either $`\tau `$ or $`\mathrm{exp}2\tau `$.

We have not included tensor perturbations, curvature or dark energy models with $`wp/\rho 1`$ in our analysis and so our results strictly only apply with these assumptions. Including all these variations (at fixed $`\mathrm{\Theta }_s`$) would only alter the $`C_l`$ at $`l60`$. Their inclusion would therefore affect our constraints on $`z`$ and $`n_S`$ (and therefore $`t`$), but not $`A^{}`$, $`\mathrm{\Theta }_s`$, $`\omega _b`$ or $`\omega _d`$.

## 3 Likelihood Calculation

What we want to know is, given the data and any other assumptions we make about the world, what is the probability distribution of the parameters? This posterior probability distribution can be calculated by use of Bayes’ theorem which states:

$$P(\stackrel{}{\theta }|d)P(d|\stackrel{}{\theta })P_{\mathrm{prior}}(\stackrel{}{\theta }).$$
(5)
where $`d`$ refers to data and the proportionality constant is chosen to ensure $`P(\stackrel{}{\theta }|d)𝑑\stackrel{}{\theta }=1`$. With a uniform prior this simply reduces to $`P(\stackrel{}{\theta }|d)P(d|\stackrel{}{\theta })`$. This probability of the data given the parameters is, when thought of as a function of the parameters, called the likelihood, $`(\stackrel{}{\theta })`$.

Often we are interested in the posterior probability distribution for one or two parameters alone. This marginalized posterior is given by integrating over the other parameters. For example:

$$P(\theta _1,\theta _2|d)=\mathrm{\Pi }_{i=3}^n𝑑\theta _i(\stackrel{}{\theta })P_{\mathrm{prior}}(\stackrel{}{\theta })$$
(6)
where $`n`$ is the number of parameters. We use the prior to incorporate non–CMB information such as that the redshift of reionization must be greater than 6.3 (Becker et al., 2001).

In the following subsections we discuss first how we evaluate the likelihood function at a single point and then how we evaluate it over a large parameter space and produce marginalized posterior distributions.

### 3.1 Likelihood evaluation

Here we take the data to be the measured averages of the angular power spectrum (called bandpowers), $`D_i^d`$. The expected signal contribution to $`D_i^d`$ is given by an average over the power spectrum:

$$D_i^s=\underset{l}{}u_{\alpha (i)}^2f_{il}𝒞_l(\stackrel{}{\theta })$$
(7)
where $`u_{\alpha (i)}`$ is the calibration parameter for dataset $`\alpha `$ and $`𝒞_ll(l+1)C_l/(2\pi )`$ is the angular power spectrum.

The uncertainty in the $`D_i`$ is non–Gaussian but well–approximated by the offset log–normal form of Bond et al. (2000). Specifically, $`=\mathrm{exp}(\chi ^2/2)`$ where

$`\chi ^2`$ $`=`$ $`{\displaystyle \underset{i,j}{}}\left(Z_i^\mathrm{s}Z_i^\mathrm{d}\right)M_{ij}^Z\left(Z_j^\mathrm{s}Z_j^\mathrm{d}\right)+\chi _{\mathrm{expt}}^2;`$ (8)
$`\chi _{\mathrm{expt}}^2`$ $``$ $`{\displaystyle \underset{\alpha }{}}{\displaystyle \frac{(u_\alpha 1)^2}{\sigma _\alpha ^2}}+(b\overline{b})^2/\sigma _b^2;`$ (9)
$`Z_i^\mathrm{d}`$ $``$ $`\mathrm{ln}(D_i^d+x_i);`$ (10)
$`Z_i^\mathrm{s}`$ $``$ $`\mathrm{ln}\left(D_i^s+x_i\right);`$ (11)
$`M_{ij}^Z`$ $``$ $`M_{ij}\left(D_i^d+x_i\right)\left(D_j^d+x_j\right)\text{no sum};`$ (12)
where $`M_{ij}`$ is the weight matrix for the band power data $`D_i^d`$. The expt label is for experimental parameters. These include the calibration parameters and a beam–width parameter, $`b`$, for the one experiment with significant and quantified uncertainty in their beam width. For simplicity, we take the prior probability distribution for the experimental parameters to be normally distributed. Since the datasets have already been calibrated, the mean of the calibration parameters is at $`\overline{u}_\alpha =1`$. The calibration parameter index, $`\alpha `$, is a function of $`i`$ since different power spectrum constraints from the same dataset all share the same calibration uncertainty.

We include bandpower data from Boomerang, the Degree Angular Scale Interferometer (DASI; Halverson et al., 2002), Maxima (Lee et al., 2001)), the COsmic Background Explorer(COBE; Bennett et al., 1996), VSA and CBI. The weight matrices, band powers and window functions are publicly available for both VSA and DASI. For COBE we approximate the window functions as tophat bands; all other information is available in Bond et al. (2000) and in electronic form with the DASh package. For Boomerang, CBI and Maxima we approximate the window functions as top-hat bands, the weight matrices as diagonal and the log-normal offsets, $`x`$, as zero. The Boomerang team report the uncertainty in their beam full-width at half-maximum (fwhm) as $`12.9\pm 1.4`$ minutes of arc. We follow them in modelling the departure from the nominal (non-Gaussian) beam shape as a Gaussian. For the Boomerang $`D_i^s`$, we replace $`𝒞_l`$ in Eq. 7 with $`𝒞_l\mathrm{exp}(l^2b^2)`$. We do not allow the fwhm to go below 11’.5 or above 14’.3; i.e. we give such fluctuations zero probability. To reduce our sensitivity to beam errors, for Boomerang and Maxima we only use bands with maximum $`l`$-values less than 1000. For CBI we use their broader mosaic’ed fields with bandpowers that extend to $`l=1900`$.

### 3.2 Exploring the Parameter Space

Our first step in exploring the high–dimensional parameter space is the creation of an array of parameter values called a chain, where each element of the array, $`\stackrel{}{\theta }`$, is a location in the $`n`$-dimensional parameter space. The chain has the useful property once it has converged that $`P(\stackrel{}{\theta }R)=N(\stackrel{}{\theta }R)/N`$ where the left–hand side is the posterior probability that $`\stackrel{}{\theta }`$ is in the region $`R`$, $`N`$ is the total number of chain elements and $`N(\stackrel{}{\theta }R)`$ is the number of chain elements with $`\stackrel{}{\theta }`$ in the region $`R`$. Once the chain is generated one can then rapidly explore one–dimensional or two–dimensional marginalizations in either the original parameters, or in derived parameters, such as $`t_0`$. Calculating the marginalized posterior distributions is simply a matter of histogramming the chain.

The chain we generate is a Monte Carlo Markov Chain (MCMC) produced via the Metropolis–Hastings algorithm described in Christensen et al. (2001). The candidate–generating function for an initial run was a normal distribution for each parameter. Subsequent runs used a multivariate–normal distribution with cross–correlations between cosmological parameters equal to those of the posterior as calculated from the initial run.

The covariance matrix of the candidate–generating function, $`C_G`$ is actually a scaled version of the posterior covariance matrix, $`C_P`$: $`C_G=aC_P`$. Proper choice of $`a`$ is important if the movement through parameter space is to be efficient. If $`a`$ is too low, the acceptance rate will be high, but the typical step size will be small and a full sampling of the parameter space will take a long time. If $`a`$ is too high then most steps will be to regions of greatly reduced probability, the acceptance rate will therefore be low and once again a full sampling of the parameter space will take a long time. As a useful measure of the speed at which we sample the parameter space, we define the step distance to be:

$$\mathrm{\Delta }\sigma _{n,n1}=\sqrt{\underset{i}{}\mathrm{\Delta }\theta _iC_{P,ij}^1\mathrm{\Delta }\theta _j};$$
(13)
where $`\mathrm{\Delta }\theta _i`$=$`(\theta _n)_i(\theta _{n1})_i`$, $`n`$ is the iteration number and $`i`$ enumerates the parameters. In other words we are using the inverse covariance matrix as the distance metric for the parameter space. We automatically adjust $`a`$ in the first 20,000 samples of a run to maximize this average step distance. We find it has a fairly broad plateau between 0.3 and 0.7. In practice we calculate the square of the distance with the slightly simpler expression: $`_i(\mathrm{\Delta }\theta _i)^2/C_{P,ii}`$.

All of our results are based on an MCMC run consisting of 754,000 iterations. For the “burn-in” the initial 20,000 samples were discarded, and the remaining set was thinned by accepting every 25th iteration, resulting in nearly 30,000 samples. The acceptance rate was 43% with $`a=.37`$ and an average step distance of 0.25 indicating highly efficient movement through the parameter space. We used the CODA software (Best et al., 1995) to confirm that the chain passed the Referty-Lewis convergence diagnostics and the Heidelberger-Welch stationarity test.

## 4 Fitting to the Likelihood

We do not fit directly to the ten–dimensional likelihood or even to its marginalization down to the six cosmological parameters. The reason is that in the high–dimensional spaces, the likelihood is very sparsely sampled. For a fairly coarse grid with 8 steps in each direction, there would be $`6^8>250,000`$ grid points and yet we have only 30,000 samples. Only with the marginalizations down to even lower–dimensions does the likelihood become densely sampled.

Instead we fit to the two–dimensional marginalized likelihoods and from these fits reconstruct the full result. We typically use 50 bins in each parameter which results in about 100 samples in the most-likely of the 2500 bins.

We do this fitting in the normal parameter space. That is, at each point in the chain we calculate the value of the two normal parameters of interest and histogram that point accordingly. For every such pair of normal parameters we find $`\overline{\theta }_i`$ and $`F_{ij}^{2d}`$ that give the best–fit two–dimensional Gaussian likelihood where

$$_{2d}^N\mathrm{exp}\left[\frac{1}{2}\left(\theta _i\overline{\theta }_i\right)F_{ij}^{2d}\left(\theta _j\overline{\theta }_j\right)\right].$$
(14)
From these $`n\times (n1)/2`$ ($`=15`$ for $`n=6`$) fits we reconstruct the six–dimensional Gaussian. We exploit a property of Gaussian distributions that the inverse of $`F_{ij}^{2d}`$ gives the elements of the covariance matrix in the full higher–dimensional space: $`C_{ij}=(F^{2d})_{ij}^1`$. Thus we uniquely determine the 15 independent off-diagonal elements of the 6 by 6 covariance matrix. The diagonal elements are over–determined with five estimates of $`C_{ii}`$ for each $`i`$; we average these together for the final estimate. The full covariance matrix is then inverted, giving the 6–dimensional Fisher-matrix of the likelihood function. The $`\overline{\theta }_i`$ are also averaged to produce an estimate of the likelihood maximum.

Note that each of the off-diagonal elements is obtained from a 2D gaussian fit to different variable pairs. Therefore there is no guarantee that the covariance matrix built from these 2x2 blocks will necessarily be positive-definite (which a true covariance matrix should be). We employ an iterative procedure to solve the above problem. We first increase the diagonal terms of $`C`$ until $`C`$ is positive-definite. We then invert $`C`$ to get the Fisher matrix $`F=C^1`$ and find its eigenmodes. Now we identify these eigenmodes as the “new normal parameters” and perform the 2D fit procedure again (this step is not computationally intensive). The resulting covariance matrix (from the 2D fits) will be close to diagonal and hence positive definite. We repeat this procedure typically about 15 times by which time the changes in the covariance matrix are less than 3%. We have also verified that the final result is insensitive to the arbitrary adjustments to the covariance matrix (to make it positive-definite) in the first step.

The end result is an approximate likelihood, $`^N`$, with most likely values $`\overline{\theta }_i`$ and covariance matrix $`C_{ij}`$:

$$2\mathrm{ln}^N=\left(\theta _i\overline{\theta }_i\right)C_{ij}^1\left(\theta _j\overline{\theta }_j\right).$$
(15)
If we denote the cosmological parameters with $`\stackrel{}{c}`$ and the normal parameters with $`\stackrel{}{\theta }`$, the likelihood of the cosmological parameters is given by

$$^c(\stackrel{}{c})=^N(\stackrel{}{\theta }(\stackrel{}{c}))|\frac{\stackrel{}{c}}{\stackrel{}{\theta }}|$$
(16)
where

$$|\frac{\stackrel{}{c}}{\stackrel{}{\theta }}|=\frac{\mathrm{\Theta }_s}{\mathrm{\Omega }_\mathrm{\Lambda }}\frac{A^{}}{A}\left(\frac{0.05\mathrm{Mpc}^1}{k_{\mathrm{pivot}}}\right)^{n_s1}b\mathrm{ln}2$$
(17)
is the Jacobian of the variable transformation.

## 5 Results

In Table 1 we show the parameters of $`^N`$ that make it a good approximation to our MCMC–calculated likelihood. For ease of interpreting the magnitude of the off–diagonal elements we have shown the correlation matrix which is given by $`C_{ij}/\sqrt{C_{ii}C_{jj}}`$ instead of $`C_{ij}`$ itself.

We compare the fit to the exact likelihood in Fig. 1 which shows 6 of the 15 possible marginalizations down to two dimensions. The marginalization of the normal fit is done by a Monte Carlo process. We use the normal fit to rapidly generate a chain whose elements are samples from the normal distribution. We can then manipulate this chain (e.g., to marginalize down to 2 dimensions) just like we manipulate the Markov chains.

All six marginalizations down to one dimension are shown in Fig. 2. One can see that the fit provides an excellent match to the exact one–dimensional likelihoods.

Note that although the fit is Gaussian in the six–dimensional space, the priors implicit in the marginalization down to one dimension result in slightly non–Gaussian one–dimensional likelihoods. This is because the marginalization is not done over the full parameter space but has the constraints $`z>6`$, $`0<\mathrm{\Omega }_\mathrm{\Lambda }<0.9`$, $`0.45<h<0.95`$ and $`0.02<\mathrm{\Omega }_m<1.1`$.

Our normal fit to the likelihood can also be used to obtain marginalized constraints on the cosmological parameters. Marginalizations down to two dimensions and one dimension are shown in Figures 3 and 4 respectively. Note that the highly non–Gaussian distributions of the cosmological parameters are described well by the fit that is Gaussian in the normal parameters.

The eigenmodes and eigenvalues $`\lambda `$ are shown in Table 2. These are the eigenmodes of the fractional Fisher matrix, $`F_{ij}^f=\overline{\theta }_iF_{ij}^a\overline{\theta }_j`$ where $`\overline{\theta }_i`$ are the best–fit values of the normal parameters. We look at the fractional Fisher matrix so its eigenvalues quantify the fractional error and not the absolute error. The eigenvectors are ordered from highest eigenvalue to lowest eigenvalue. We see that current CMB data constrain four parameters at the 10% level or better. As expected, the best–constrained eigenvector is predominantly $`\mathrm{\Theta }_s`$ and the worst–constrained is almost entirely $`z`$. The second–from–worst is $`\omega _d`$, though even it is farily well–constrained with $`\sqrt{\lambda }=0.16`$. Eigenvectors 2, 3 and 4 are mostly $`t`$, $`A^{}`$ and $`\omega _b`$ respectively, though each with significant amounts of other parameters mixed in.

We have also explored the likelihood distribution of the parameter set proposed by Kosowsky et al. (2002). In 2-D projections, most pairs of parameters do look Gaussian. An exception is any pair involving their variable $`𝒵\mathrm{exp}(2\tau )`$, which is significantly less Gaussian than the variable we have used to parameterize reionization: the redshift of reionization, $`z`$.

## 6 Discussion and Conclusions

We have defined a set of parameters that, given current data, have a likelihood that is normally distributed. Of course, our simple fit to the likelihood surface will very soon be out-of-date due to the expected MAP data. Still, we expect the technique to continue to be useful. Generally, as data become more constraining, the Gaussian approximation becomes better.

The challenge in the future will be to provide good parameterizations for those parameters that remain poorly determined. Some of these only have their effects at low $`\mathrm{}`$, such as the parameters governing reionization. For these it may be useful to define parameters for which $`\mathrm{ln}C_l`$ is linear, since the errors in this quantity are very nearly normally distributed at low $`\mathrm{}`$ Bond et al. (2000). Other physical effects, such as gravitational lensing, are only important at high $`\mathrm{}`$ where, if instrument noise is the dominant contribution to error in $`C_l`$, then $`C_l`$ will be nearly normally distributed. For these one would want to use parameterizations for which $`C_l`$ is linear.

An analytic form for the likelihood makes the actual application of CMB parameter constraints to cosmological questions much easier. One application is to quickly calculate constraints on various combinations of the cosmological parameters such as $`\sigma _8`$ (Holder, 2002) or the age, $`t_0`$ (Knox et al., 2001). Another is to combine the CMB results with those from other probes to improve constraints or search for inconsistencies (Wang et al., 2002). Or one can forecast expected constraints from a combination of the CMB and other probes such as supernovae (Frieman et al., 2002) or galaxy cluster number counts (Levine et al., 2002). Although these applications are possible without an analytic form for the CMB likelihood, such a form greatly simplifies the analysis.

Use of normal parameters will also improve the efficiency of Monte Carlo evaluations of the likelihood, by providing an analytic generating function that is well–matched to the actual likelihood. Although the statistical properties of the post–convergence chain do not depend on the generating function used, the time required for convergence depends critically on the generating function. In particular the chain rapidly converges if the generating function approximates the likelihood well. Clearly, this issue becomes more important for larger parameter spaces. A practical way to implement the use of normal parameters in generating Monte Carlo Markov chains would be to create a small chain, derive the normal parameters and the gaussian approximation to the likelihood, and then use this gaussian likelihood as the generating function.

In this paper we have explicitly given an analytic form for the likelihood (given a certain data set) of cosmological parameters and shown it to be an excellent approximation to the MCMC–calculated likelihood. We view this as the final step in a long data-reduction chain that begins with the time–ordered data.

We thank N. Christensen, A. Kosowsky and R. Meyer for useful conversations. This work was supported in part by NASA grant NAG5-11098.


format: text

text:
# Spectroscopy of Globular Clusters in the Fornax dwarf galaxy

## 1 Introduction

The Fornax dwarf spheroidal (dSph) galaxy has one of the most complex star formation histories of any of the Local Group dSphs (Stetson et al., 1998). While a few of these galaxies seem to have only old stars (e.g., Draco, Sculptor), such dSphs as Leo I and Carina show star formation up to $`2`$ Gyr ago (Dolphin, 2002). In contrast, several studies (Stetson et al. 1998; Buonanno et al. 1999, hereafter B99) have found an age range of $``$ 0.5–12 Gyr in the Fornax field, dominated by intermediate-age stars. This complex history could be due in part to the galaxy’s high luminosity (and mass); at M$`{}_{V}{}^{}13.5`$ it is second only to the Sagittarius dSph.

Fornax is unique among the Local Group dSphs in that it has a significant population of globular clusters—five—which results in a specific frequency $`S_N20`$. This is the highest $`S_N`$ known for any galaxy. However, increasing evidence suggests that other dSphs may have once had their own GCs as well. Sagittarius, for example, is now known to have at least four GCs (M54, Ter 4, Ter 7, Arp 2). Dinescu et al. (2000) and Martinez-Delgado et al. (2002) have argued that the tidally disrupting GC Pal 12 may be associated with Sagittarius as well. More speculatively, Yoon & Lee (2002) have proposed that a subset of the Oosterhoof II Milky Way GCs (including M92 and M68) may have been tidally stripped from the Large Magellanic Cloud. This suggests that GC populations in dSph/dIrr galaxies may have been more common is the past, but that Fornax was “dynamically lucky” in that its distance from the Galactic Center ($`140`$ kpc) may have prevented similar stripping of its GCs. Thus a close study of the Fornax GC system (GCS) could shed light on what may have been a common phenomenon among dwarf galaxies. While the discovery of other such systems would not explain the large difference in $`S_N`$ between Fornax and, for example, regular dwarf spheroidals, it would help address why Fornax currently appears to be a singular object.

For its small size, the Fornax GCS is surprisingly complex. For example, Rodgers & Roberts (1994) found that while the radial profiles of clusters 1 and 2 (hereafter we refer to the GCs as H1–H5, according to Hodge 1961) are well-fit by King profiles with tidal truncation (though the H1 sample was poor), they could not find any evidence of such truncation in H3–H5. While it is difficult to distinguish between a large tidal radius and no truncation at all, even a large difference in tidal radii between the two groups is not easily explained. There appear to be no simple correlations between these properties and galactrocentric distance, metallicity, or age. The radial profile of H5 is even more complex: Mackey & Gilmore (2003) find that the core of the cluster is best fit by a power law (rather than the flat, constant density core of a King model) and suggest that H5 may be a post core-collapse GC.

Oh et al. (2000) show that the expected timescale for mass segregation in the galactic potential is $`1`$ Gyr and that the GCs should have sunk to the center of the galaxy. That this is not the case is evidence that some other mechanism is kinematically exciting the clusters into their present orbit. They suggest that the most likely scenario is one in which tidal effects from the Galaxy torque the cluster orbits. Though dynamical friction within Fornax will counteract this force, the present widely spread GC locations would be possible only if the galaxy were undergoing significant mass loss (at least 30% of the original total). This might indicate that a Hubble time ago Fornax was much more massive, possibly explaining its high GC specific frequency. Were this scenario true, we would expect to see tidally stripped stars with high proper motions (and chemical similarities to Fornax field stars) in the area surrounding the galaxy. No such studies have yet been published. A recent proper motion study of Fornax (Piatek et al., 2002) suggests that it may not even be bound to the Galaxy, although a bounded low-eccentricity or even circular orbit is consistent with the data if the galaxy is near perigalacticon. The Keplerian period for such an orbit would be 4–8 Gyr; the hypothetical mass loss could have peaked during the last perigalactic passage.

In this paper, we present a spectroscopic survey of the Fornax GCs, deriving ages and metallicities, with a view to constraining the formation and evolution of the galaxy.

## 2 Observations and Data Reduction

The observations were made using the Low Resolution Imaging Spectrometer (Oke et al., 1995) in longslit mode on the Keck I telescope. All five clusters were observed during a two night run in December 1996. A single 300 s exposure was taken for each of the clusters. These observations were supplemental to our main program. The $`V`$ band magnitudes of the five GCs are H1 (15.6), H2 (13.5), H3 (12.6), H4 (13.4), and H5 (13.4), all from Webbink (1985).

All observations were made with 600 line mm-1 gratings, blazed at 5000 Å. This gave a spectral resolution of $``$ 6 Å and a dispersion of 1.28 Å pixel-1. The total wavelength range was generally 3800–6200 Å, but the usable range was often shorter by several hundred angstroms, due both to low signal-to-noise (S/N) ratios for several objects as well as the poor instrumental response of the LRIS “red arm” at the blue end of the spectrum (these spectra were taken before installation of the “blue arm”).

The standard data reduction was performed using IRAF111IRAF is distributed by the National Optical Astronomy Observatories, which are operated by the Association of Universities for Research in Astronomy, Inc., under cooperative agreement with the National Science Foundation.. Raw images were debiased and then flatfielded, using a normalized composite sky flat. After cosmic ray removal using the IRAF task *cosmicrays*, the spectra were wavelength calibrated with HgKr comparison arc lamp spectra. These spectra were flux-calibrated using the flux standard Feige 24. The spectra were corrected to zero redshift before the measurement of any line indices. For reference, Figure 1 shows the spectra of the five GCs.

## 3 Lick/IDS Indices

Elemental abundances were measured using Lick/IDS absorption-line indices (Trager et al. 1998 and references therein), from spectra convolved with wavelength-dependent Gaussian kernel to more closely match the Lick/IDS resolution. See Larsen & Brodie (2002) and Strader et al. (2003) for more extensive discussions of this issue.

We first measured a set of Lick/IDS indices (listed in Table 1) for comparison to stellar evolutionary models. Though we did measure indices for the faint cluster H1, the errors are too large for their inclusion in the tables, and we exclude H1 from subsequent analysis. To check that the indices were self-consistent, we plotted the major metallicity indicators (Mg2, Mgb, Fe5335) against each other. In each case, within the error bars, the indices fell along a 1:1 relation.

Our measurements of Lick/IDS Mg and Ca features have large index errors, translating into 0.3–0.4 dex errors in \[$`\alpha `$/Fe\], making it impossible to determine whether a given cluster is truly $`\alpha `$-enhanced. An additional complicating factor is that $`\alpha `$-enhanced SSP models (e.g., Milone et al. 2000; Thomas et al. 2002) converge at very low metallicities. For these reasons we were unable to make any accurate estimates of \[$`\alpha `$/Fe\] for the Fornax GCs.

## 4 Metallicities

Brodie & Huchra (1990, hereafter BH90) defined a procedure for estimating the metallicity of globular clusters with low resolution spectra by taking a weighted mean of various elemental absorption-line indices sensitive to \[Fe/H\]. Indices defined over a wavelength range of interest are calculated with respect to a pseudocontinuum, defined using regions to either side of the feature bandpass.

No resolution correction was made for measuring BH90 metallicity indices, since the BH90 calibration primarily used data near our resolution or better. Then, as described in §2, the spectra were flux calibrated before any abundance measurements were made. The minor role of flux calibration in measuring indices has been extensively documented (Faber et al., 1985; Kissler-Patig et al., 1998; Larsen & Brodie, 2002). To verify this, we measured several strong features on the spectra both before and after flux calibration. The differences between the fluxed and unfluxed line indices were generally $``$ 3–5% of the error bars, and thus insignificant.

We estimated \[Fe/H\] for clusters H2–H5 according to the formulae in BH90, with the exception of a slight modification in calculating the variance, as described in Larsen & Brodie (2002). The feature indices used in the spectroscopic measurements are listed in Table 2. The resulting individual index \[Fe/H\] estimates, as well as the composite values, are listed in Table 3.

## 5 Ages

To estimate the ages of the GCs in our sample, we plotted various combinations of the age-sensitive indices (e.g H$`\beta `$ and H$`\gamma `$) against the metallicity features (e.g., Mg2, Mgb, Fe5270, Fe5335). Figure 2 shows H$`\beta `$ plotted against a combined Fe5270-Fe5335 index (denoted $`<`$Fe$`>`$) for clusters H2–H5. Superimposed on the plots are models obtained from Thomas et al. (2002). The GCs H2–H4 are coeval within the errors, while H5 appears to be younger by a couple of Gyr at the 2$`\sigma `$ level. An H$`\beta `$ vs \[MgFe\] diagram (Figure 3) supports the slightly younger age of H5. One advantage to using the composite index \[MgFe\] is that is it insensitive to \[$`\alpha `$/Fe\] variations. The lack of *absolute* agreement between the index measurements and Simple Stellar Population (SSP) models may be due to overestimation of the rise in H$`\beta `$ absorption at very low metallicities (R. Schiavon, private communication). These diagrams provide support at the $`2\sigma `$ level for an age spread in the Fornax GCs, with H5 several Gyr younger than the other clusters, despite its relatively low metallicity.

Our Fornax data are compared to the Lick-corrected Large Magellanic Cloud (LMC) sample of globular clusters from Beasley, Hoyle & Sharples (2002) on the Thomas et al. (2002) models in Figures 2 & 3. In the H$`\beta `$-$`<`$Fe$`>`$ planes the Fornax clusters span similar ranges in metallicity, although in the mean the LMC sample is slightly more metal-rich. Two of the LMC clusters show strong H$`\beta `$ absorption on the SSP grids; Beasley et al. attribute this to blue horizontal branch (HB) effects rather than younger ages. The Fornax cluster H5 also exhibits a slightly elevated H$`\beta `$ line-strength, with a very similar metallicity (and presumably age) to the LMC cluster NGC 1916. The Buonanno et al. (1998) CMD of H5 clearly shows a blue HB. In principle, the stronger H$`\beta `$ index, and hence younger inferred age of the cluster, could be due to the presence of the blue HB (e.g., Lee et al. 2000). However, we note that H2 and H3 have very similar blue HB morphologies and yet have relatively weak H$`\beta `$ lines. Our Fornax data is not well-constrained by the SSP models in Figure 3, since three of the four clusters lie to the left/below the grids. However, our general conclusion is that these GCs are old, with cluster H4 somewhat more metal-rich than the rest.

To probe these results, we undertook a similar exercise using the age-sensitive indices H$`\gamma `$ and H$`\delta `$, with inconclusive results. These indices are can be less dependent on metallicity, if defined properly (for example, see the H$`\gamma `$ definition in Vazdekis & Arimoto 1999). However, this reduced metal-sensitivity is only possible because the index range is very small (less than one Angstrom). This, in turn, means that the spectra must be very high S/N to obtain an accurate index value. Another complicating factor is that there is a significant age-metallicity degeneracy in current SSP models for these indices at very low \[Fe/H\] values, and thus the measured H$`\gamma `$ and H$`\delta `$ indices for our clusters are compatible with isochrones from 7–14 Gyr. The indices H$`\gamma _A`$ and H$`\delta _A`$, as defined in Worthey & Ottaviani (1997), are listed in Table 4.

For GCs H2-H5 we have also derived a metallicity estimate using the H$`\beta `$ vs $`<`$Fe$`>`$ model grid. In Table 5 we list this estimate (called the “Lick Indices” method) along with that derived from the BH90 method (see §4).

## 6 Comparison with Previous Works

One immediately noticeable feature in Table 3 is how little spread there is in our metallicity estimates, amounting to only $`0.3`$ dex among clusters H2–H5. This is in contrast to the larger spread ($``$ 0.6–0.9 dex) found by some early investigators (e.g., Dubath et al. 1992), which was primarily driven by the high metallicity of H4 with respect to the other clusters. Compared to the Dubath et al. \[Fe/H\] estimates, we find that H3 and H5 are slightly more metal-rich, with \[Fe/H\] = $``$1.8 and $``$1.7 respectively, while H4 is somewhat more metal-poor at \[Fe/H\] = $``$1.5. The sum of these small shifts substantially reduces the overall dispersion. Were we able to include H1, it would likely not affect this result, since previous investigations uniformly have found that within the error bars H1 and H5 have the same metallicity (Dubath et al., 1992; Buonanno et al., 1985, 1998). For comparison purposes, we have listed our \[Fe/H\] estimates together with past measurements using a variety of methods in Table 5. Most previous researchers find H4 to be slightly more metal-rich that the other four GCs. The case of H4 is discussed further in §6.1.

Several of the less common methods in Table 5 merit explanation: Dubath et al. (1992) cross-correlated the cluster spectrum with that of a template; calculating the area under the peak of the cross-correlation function is a measure of the effective temperature of the RGB of the cluster and hence of the metallicity. Demers et al. (1990) converted the $`BV`$ color of H1 to a metallicity using the Zinn & West (1984) relation. Finally, the “Zinn Q” method of Zinn & Persson (1981) uses a transformation between the cluster metallicity and the quantity $`Q_{39}`$, which is a composite photometric measure of the curvature and line blanketing in the cluster spectrum.

In terms of GC ages, we found the GCs to be uniformly old with H5 perhaps a few Gyr younger. The primary alternative sources of Fornax GC age measurements are the color-magnitude diagrams (CMDs), the most accurate of which are constructed from HST images (B99, Buonanno et al. 1998). To compare our age measurements with the CMD methods, we chose the cluster pair H3 and H5, as the error bars on H2 are somewhat larger and unknown parameters may be affecting the H4 CMD. Buonanno et al. (1998) found essentially no age difference between the two clusters using either the $`\mathrm{\Delta }V_{\mathrm{HB}}^{\mathrm{TO}}`$ method ($`0.1\pm 2.8`$ Gyr) or the $`\mathrm{\Delta }(VI)_{\mathrm{RGB}}^{\mathrm{TO}}`$ method ($`0.1\pm 0.8`$ Gyr). However, the error bars on the measurements include *only* errors in the measurement of $`V_{\mathrm{HB}}`$ and $`V_{\mathrm{TO}}`$. If other sources of error are included, such as the uncertainty in the HB MV-\[Fe/H\] relation, the possible age spread of 2–3 Gyr found by this work agrees with the CMD findings within the error bars.

### 6.1 The Strange Case of H4

Early spectroscopic and photometric studies of H4 indicated it was relatively metal-rich at \[Fe/H\] $`1.4`$, compared to the other Fornax GCs, which seemed to be uniformly metal-poor (Harris & Canterna, 1977; Zinn & Persson, 1981; Dubath et al., 1992). Interestingly, this value is close the mean of the Fornax field (\[Fe/H\] = $``$1.36; Buonanno et al. 1999). Detailed studies of the CMD of the cluster were difficult to make due to its high density and location near the center of the galaxy. Ground-based attempts failed to achieve a CMD accurate enough to substantially constrain the metallicity, although CMD metallicity estimates were lower than the spectroscopic values (Buonanno et al., 1985; Beauchamp et al., 1995). However, B99 were able to construct a relatively clean CMD using HST Wide Field and Planetary Camera 2 (WFPC2) images. By comparing the slope of the cluster’s RGB to that of Galactic GCs, they found the metallicity of H4 was \[Fe/H\] = $`2.01\pm 0.20`$, in sharp disagreement with earlier measurements and our work.

B99 considered several explanations for the discrepancy. The first was a significantly larger reddening for H4 than for the other clusters, but by comparing the H4 CMD to that of H2, they argued that the amount of reddening needed to account for the apparent metallicity difference in the two methods was improbably high. Field star contamination was also considered as a possible factor, since the surrounding field stars are of a higher metallicity than the cluster stars (as ascertained from the HST CMD). This might especially be an issue for H4 as it lies close to the center for the galaxy.

Their third explanation concerns the nature of the relationship between \[Fe/H\] and the RGB. Noting that the correlations between broadband colors or elemental line indices and the metallicity of the cluster was dependent primarily on the nature of the RGB, and that the calibration of these relations is based upon Galactic GCs, they argue that were H4 somehow different from the bulk of Galactic GCs, such relationships might not hold. Thus, the spectroscopic \[Fe/H\] estimates might be invalid. To support this view, they show that the shape of the H4 CMD, with an extremely red HB and a relatively bright SGB, differs substantially from that of the other Fornax GCs as well as that of a “classic” metal-poor Milky Way GC, M68, in that the CMDs cannot be shifted to fit each other with reasonable metallicity and age differences. B99 suggest that an anomalous HB and SGB are the source of the discrepancy in that they cause the integrated cluster spectrum to appear to be more metal-rich. While a parameter which *only* affected the HB and SGB and not the slope of the RGB could plausibly exist, B99 offer no suggestions.

B99 describe their fourth explanation of the difference as “speculative.” They note that with the appropriate reddening corrections, the CMD of H4 is a close fit to that of Ruprecht 106 (hereafter R106), an outer halo Milky Way GC believed to be 3-4 Gyr younger than the typical old, metal-poor Galactic GC. Indeed, there has been a similar discrepancy between the spectroscopic and RGB slope \[Fe/H\] measurements for this cluster, with the former method yielding a much higher metallicity than the latter (Brown et al., 1997; Buonanno et al., 1993)

Importantly, Brown et al. (1997) also found that R106 had anomalous values of \[O/Fe\] and \[$`\alpha `$/Fe\]. The former measurement yielded a value of \[O/Fe\] = $`0.0\pm 0.13`$. In a typical metal-poor cluster \[O/Fe\] $``$ 0.3–0.4, although there is a great deal of variation both from cluster to cluster and within some individual GCs (Carney, 2001). A similar difference is seen in the $`\alpha `$-element abundance—Brown et al. (1997) found \[$`\alpha `$/Fe\] = $`0.05\pm 0.05`$, compared to a typical value of \[$`\alpha `$/Fe\] $`0.3`$ (Carney, 1996). Brown et al. (1997) also found an $`\alpha `$-element underabundance in the young outer-halo Galactic GC Pal 12, the only other Milky Way cluster for which a lack of $`\alpha `$-element enhancement is seen (Carney, 2001). No complete, deep $`VI`$ CMD of Pal 12 currently exists (the closest is Rosenberg et al. 1998).

These results suggest a detailed examination of \[$`\alpha `$/Fe\] and \[O/Fe\] values in H4 is in order. The two common methods for measuring \[O/Fe\] in red giants and in horizontal branch stars (forbidden \[OI\] lines at $`6300`$ and 6344 Å and the IR triplet) both require high-resolution spectroscopy. Furthermore, the two methods are at odds with each other, differing by $`0.6`$ dex at low metallicities (Carney, 2001). High-resolution spectra and a solution to this discrepancy will be needed before accurate calculations and comparisons of \[O/Fe\] can be made.

A qualitative argument for a higher metallicity for H4 can be seen in Figure 1. The metal lines in H4 are clearly stronger than those in H5, e.g., the Mg2 and Fe5270 index bands. Note in particular the increased strength of the G Band, which measures CH, and is closely correlated with \[Fe/H\] (Brodie & Huchra, 1990). Based *only* on our data, the metallicity difference between H4 and the other clusters is significant at the 2$`\sigma `$ level.

In sum, we believe that H4 is $`0.30.4`$ dex more metal-rich than the other Fornax GCs, and that the discrepancy between our spectral \[Fe/H\] estimate (as well as those of previous spectroscopic and photometric studies; see Table 5) and that based on the RGB slope could be caused by a solar-like \[$`\alpha `$/Fe\] ratio in the cluster. High resolution spectra are needed to confirm this suggestion.

We offer a final, alternative explanation of the nature of H4 that could resolve the controversy. We find that H4 is old like clusters H1–H3, and yet slightly more metal-rich. It is also known to be the most compact cluster and is centrally located. We speculate that it is actually the nucleus of the Fornax dSph galaxy (this same possibility has been suggested by Hardy 2002). As such, it would be comparable to the Galactic GC M54, which has been suggested to be the remnant nucleus of the Sagittarius (Sgr) dSph galaxy. Like H4, M54 (Layden & Sarajedini 2000) is as old as Terzan 8 and Arp 2 (two other Sgr GCs) and slightly more metal-rich (though this result is not statistically significant). Interestingly, Layden & Sarajedini find evidence for multiple Sgr field star formation episodes with ages of 11 and 5 Gyr, with the most recent burst occurring 0.5 Gyr ago. As we shall discuss in §7, Fornax has a similarly bursty star formation history.

## 7 The Fornax Star Formation History and Discussion

Buonanno et al. (1999) found that star formation in Fornax began $`12`$ Gyr ago, with bursts occurring $``$7, 4 and 2.5 Gyr ago and with the most recent episode in the last $``$0.5 Gyr. The wide field imaging study of 40,000 stars by Saviane et al. (2000) indicates that the oldest population represents a relatively small component of the galaxy, with the dominant intermediate population having a mean age of 5.4 $`\pm `$ 1.7 Gyr. Saviane et al. also found that the young burst continued until as recently as 200 Myr ago. This latest star formation episode may be the last, as Fornax appears to have little or no ISM gas left (Young 1999).

Recently, Gallert et al. (2002) obtained spectra of the Ca II triplet feature for individual RGB stars using VLT/FORS. The resulting age-metallicity relation shows a rapid increase in metallicity from \[Fe/H\] $``$ –2 to –1 at 12 Gyr, followed by a near constant \[Fe/H\] $``$ –1 until a few Gyr ago.

In terms of spatial distribution, the oldest population is the most extended, indicating that the initial burst occurred over the whole galaxy. Subsequent bursts were progressively more centrally concentrated, perhaps forming as gas was funneled toward the galaxy center. This indicates an outside-in formation process extending over several Gyr. It is possible that interactions between Fornax and the Milky Way provided the mechanism for such a process, but defining the past orbits of these galaxies is difficult (see Piatek et al. 2002 for a recent HST proper motion study of the galaxy).

The GCs H1–H3 and H5 are extremely metal-poor (using previous results for H1) and all can likely be associated with the oldest field star population, despite the tentative relative youth of H5. A similar situation is seen with old GCs in other Local Group galaxies (e.g., the Large Magellanic Cloud). H4 has been discussed earlier—it is centrally located and may be slightly more metal-rich than the other GCs, so it is tempting to associate it with the second burst of star formation. However, our analysis does not provide evidence that it is younger than the other GCs. If we associate the four metal-poor GCs or all five of the GCs with the oldest field population, which represents $``$ 30% of the galaxy light (Saviane et al. 2000), then the specific frequency (SN) of this population rises from $`20`$ to 50–65.

There are several possible interpretations of this result. Assuming that dwarf galaxies are the building blocks of more massive galaxies (as argued in hierarchical scenarios), then it suggests that early protogalactic fragments may have formed GCs with efficiencies much higher than today. Alternatively, Harris (2001) has stressed the interpretation of Fornax and other high SN dwarf galaxies as being *field star poor* rather than GC rich. In one possible scenario, winds from massive stars and supernovae associated with the first generation of Fornax stars blow a substantial portion of the remaining gas out of the galaxy’s shallow potential well, which temporarily halts star formation. The subsequent cooling and infall of part of this gas (which may be dynamically connected to interactions with the Milky Way) then forms the basis for later bursts of star formation.

The present study seems to raise more questions than it answers. Are $`\alpha `$-element anomalies in H4 responsible for the discrepancy between spectroscopic and CMD \[Fe/H\] estimates? If so, what is the cause of these anomalies? A classic scenario would be that H4 formed at the tail end of a starburst, after enrichment of the interstellar medium by both $`\alpha `$-element rich Type II supernovae (SNe) and the longer-lived, Fe-bearing Type Ia SNe. This might account for its higher metallicity, but begs the question as to why it is *older* than the presumably $`\alpha `$-enhanced GC H5.

We could plausibly conjecture that the Fornax dSph was constructed out of several Searle-Zinn fragments (Searle & Zinn, 1978), each of which may have produced one or more GCs and had its own chemical enrichment period before coalescing into a larger object (see Stetson et al. 1996 for a recent discussion of the varied Searle-Zinn models). Such a scenario might then separate out the age and abundance issues discussed above.

Future high resolution spectroscopy of Fornax GC red giants, field stars, and possible candidates for stripped stars will be necessary to untangle the complex interplay of dynamical and chemical evolution in the galaxy.

## 8 Conclusions

We have analyzed spectra of the five globular clusters of the Fornax dSph. We find that clusters H2, H3, and H5 are very metal-poor with \[Fe/H\] $`1.8`$. We discuss the controversial case of cluster H4 and conclude that it is slightly more metal-rich than the other four GCs with \[Fe/H\] $``$ –1.5, thus supporting previous spectral and photometric results but contrary to recent CMD-based claims. We also argue for $`\alpha `$-element anomalies in H4, similar to the outer-halo Milky Way GCs Ruprecht 106 and Pal 12. We speculate that H4 may be the nucleus of Fornax, comparable to M54 in the Sagittarius dSph. There is evidence for an age spread among the five clusters, with H5 possibly younger than the other clusters by a few Gyr.


text:
# Thermal Pions at Finite Isospin Chemical Potential

## Abstract

The density corrections, in terms of the isospin chemical potential $`\mu _I`$, to the mass of the pions are studied in the framework of the $`SU(2)`$ low energy effective chiral lagrangian. The pion decay constant $`f_\pi (T,\mu _I)`$ is also analized. As a function of temperature for $`\mu _I=0`$, the mass remains quite stable, starting to grow for very high values of $`T`$, confirming previous results. However, there are interesting corrections to the mass when both effects (temperature and chemical potential) are simultaneously present. At zero temperature the $`\pi ^\pm `$ should condensate when $`\mu _I=m_\pi `$. This is not longer valid anymore at finite $`T`$. The mass of the $`\pi _0`$ acquires also a non trivial dependence on $`\mu _I`$ due to the finite temperature.

Pions play a special role in the dynamics of hot hadronic matter since they are the lightest hadrons. Therefore, it is quite important to understand not only the temperature dependence of the pion’s Green functions but also their behavior as function of density, through the chemical potential. The dependence of the pion mass and decay constant on temperature $`m_\pi (T)`$, $`f_\pi (T)`$ has been studied in a variety of frameworks, such as thermal QCD-Sum Rules DFL , Chiral Perturbation Theory (low temperature expansion) GL , the Linear Sigma Model LCL , the Mean Field Approximation Bar1 , the Virial Expansion Schenk , etc. In fact, the pion propagation at finite temperature has been calculated at two loops in the frame of chiral perturbation theory Schenk2 ; T . There seems to be a reasonable agreement that $`m_\pi (T)`$ is essentially independent of $`T`$, except possibly near the critical temperature $`T_c`$ where $`m_\pi (T)`$ increases with $`T`$ and that $`f_\pi (T)`$ vanishes for the critical temperature.

The introduction of in-medium processes via isospin chemical potential has been studied at zero temperature son ; toublan ; splittorff  in both phases ($`|\mu _I|m_\pi `$) at tree level. The problem with both, temperature and density, has been worked out for barionic chemical potential with Chiral Perturbation Theory AE-N . It is also possible to find certain region of the stable pion gas in which the pion number is locally conserved ayala .

Usually, there are two procedures to extract the information of $`m_\pi `$ and $`f_\pi `$ in the frame of chiral perturbation theory. The first one is to compute the Axial-Axial correlator which provides us with the decay constant and the mass corrections. GL ; GL2 ; Schenk2

$$d^4xe^{ipx}0|A_\mu ^a(x)A_\nu ^b(0)|0=\delta ^{ab}\frac{p_\mu p_\nu f_\pi ^2}{p^2m_\pi ^2}$$
(1)
In the second method, radiative corrections to the propagators are considered together with the realization of PCAC, $`0|A_\mu ^a|\pi ^b(p)=ip_\mu \delta ^{ab}f_\pi `$ making then use of appropriate counterterms. The use of counterterms is not necessary in the Axial-Axial correlator method. We have checked that both methods leave the same answers. Let us proceed in the frame of the $`SU(2)`$ chiral perturbation theory. The most general chiral invariant expression for a QCD-extended lagrangian, GL2 ; pich  under the presence of external hermitian-matrix auxiliary fields, has the form

$`_{QCD}(s,p,v_\mu ,a_\mu )`$ $`=`$ $`_{QCD}^0`$ (2)
$`+`$ $`\overline{q}\gamma ^\mu (v_\mu +\gamma _5a_\mu )q`$
$``$ $`\overline{q}(si\gamma _5p)q,`$
where $`v_\mu `$, $`a_\mu `$, $`s`$ and $`p`$ are vector, axial, scalar and pseudoscalar fields. The vector current is given by

$$J_\mu ^a=\frac{1}{2}\overline{q}\gamma _\mu \tau ^aq.$$
(3)
When $`v,a,p=0`$ and $`s=M`$, being $`M=diag(m_u,m_d)`$ the mass matrix, we obtain the usual $`QCD`$ lagrangian. This procedure is formal, in the sence that we reproduce the usual QCD lagrangian with current masses. However, we would like to notice that a scalar field in chiral lagrangian models the spontaneus break of chiral symmetry through a non vanishing vacuum expectation value. In this sense if we take for $`s=M`$, these masses should be actually constituent quark masses, while in the QCD lagrangian we have current masses. Nevertheless this is a formal step which tries only to motivate what follows in the context of effective pion lagrangian.

The effective action with finite isospin chemical potential is given by

$`_{QCD}^I`$ $`=`$ $`_{QCD}(M,0,0,0)+\mu ^au^\mu J_\mu ^a`$ (4)
$`=`$ $`_{QCD}(M,0,\mu u_\mu ,0)`$
where $`\mu ^a=(0,0,\mu _I)`$ is the third isospin component, $`\mu =\mu ^a\tau ^a/2`$ and $`u_\mu `$ is the 4-velocity between the observer and the thermal heat bath. This is required in order to describe in a covariant way this system, where the Lorentz invariance is broken since the thermal heath bath represents a privileged frame of reference.

Proceeding in the same way, now in the low-energy description where only pion degrees of freedom are relevant, let us consider the most general chiral invariant lagrangian ordered in a series of powers of the external momentum. We will start with the $`𝒪(p^2)`$ chiral lagrangian

$`_2`$ $`=`$ $`{\displaystyle \frac{f^2}{4}}Tr\left[(D_\mu U)^{}D^\mu U+U^{}\chi +\chi ^{}U\right]`$ (5)
with

$`D_\mu U`$ $`=`$ $`_\mu Ui[v_\mu ,U]i\{a_\mu ,U\}`$
$`\chi `$ $`=`$ $`2B(s+ip)`$
$`U`$ $`=`$ $`\overline{U}^{\frac{1}{2}}(e^{i\pi ^a\tau ^a/f})\overline{U}^{\frac{1}{2}}`$ (6)
$`\overline{U}`$ is the vacuum expectation value of the field $`U`$ and $`B`$ in the previous equation is an arbitrary constant which will be fixed when the mass is identified setting $`(m_u+m_d)B=m^2`$. The most general $`𝒪(p^4)`$ chiral lagrangian has the form

$`_4`$ $`=`$ $`\alpha _1\left(Tr\left[(D_\mu U)^{}D^\mu U\right]\right)^2`$ (7)
$`+`$ $`\alpha _2Tr\left[(D_\mu U)^{}D_\nu U\right]Tr\left[(D^\mu U)^{}D^\nu U\right]`$
$`+`$ $`\alpha _3\left(Tr\left[\chi U^{}+U\chi ^{}\right]\right)^2`$
$`+`$ $`\alpha _4Tr\left[(D_\mu U)^{}D^\mu U\right]Tr\left[\chi U^{}+U\chi ^{}\right]`$
$`+`$ $`\alpha _5\left[L_{\mu \nu }UR^{\mu \nu }U^{}\right]`$
$`+`$ $`i\alpha _6Tr\left[L_{\mu \nu }D^\mu U(D^\nu U)^{}+R_{\mu \nu }(D^\mu U)^{}D^\nu U\right]`$
$`+`$ $`\alpha _7\left(Tr\left[\chi U^{}U\chi ^{}\right]\right)^2`$
$`+`$ $`\alpha _8Tr\left[\chi U^{}\chi U^{}+U\chi ^{}U\chi ^{}\right]`$
$`+`$ $`\alpha _9Tr\left[L_{\mu \nu }L^{\mu \nu }+R_{\mu \nu }R^{\mu \nu }\right]`$
$`+`$ $`\alpha _{10}Tr[\chi ^{}\chi ]`$
with

$`L_{\mu \nu }=_\mu l_\nu _\nu l_\mu +i[l_\mu ,l_\nu ],`$ $`l_\mu =v_\mu a_\mu `$
$`R_{\mu \nu }=_\mu r_\nu _\nu r_\mu +i[r_\mu ,r_\nu ],`$ $`r_\mu =v_\mu +a_\mu `$ (8)
The different coupling constants $`\alpha _i`$in the previous expression are related to the couplings introduced by GL1 . Here we use the prescription of Scherer .

The effective action with finite chemical potential in terms of pion degrees of freedom has the same form as eq.4, where the different external fields are defined in eq.6. In this paper we will consider one loop corrections, up to the fourth order in the fields, to the lagrangian $`_2`$ and the free part, i.e the tree level part of $`_4`$ with renormalized fields. This procedure is standard, GL2 ; holstein . We will concentrate on the phase where $`\mu _I<m_\pi `$, where the vacuum expectation value $`\overline{U}=1`$. The interacting part $`_4`$ involves higher powers in the momentum of the pion fields. The constants $`\alpha _i`$ present in $`_4`$ are known from decay and scattering measurements. Therefore, we have the following lagrangians $`_{i,j}`$

$`_{2,2}`$ $`=`$ $`{\displaystyle \frac{1}{2}}\left[\left(\pi _0\right)^2m^2\pi _0^2\right]+\left|_I\pi \right|^2m^2\left|\pi \right|^2`$
$`_{2,4}`$ $`=`$ $`{\displaystyle \frac{1}{4!}}{\displaystyle \frac{m^2}{f^2}}\pi _0^4+{\displaystyle \frac{1}{6f^2}}[4|_I\pi |^2|\pi |^2`$
$`+(|\pi |^2)^2+m^2\left(\right|\pi |^2)^2]`$
$`+{\displaystyle \frac{1}{6f^2}}[2|_I\pi |^2\pi _0^22(\pi _0)^2|\pi |^2`$
$`+\pi _0^2|\pi |^2+m^2\pi _0^2|\pi |^2]`$
$`_{4,2}`$ $`=`$ $`2{\displaystyle \frac{m^2}{f^2}}\left[l_4\right|_I\pi |^2m^2(l_3+l_4)|\pi |^2`$ (9)
$`+{\displaystyle \frac{1}{2}}l_4(\pi _0)^2{\displaystyle \frac{1}{2}}m^2(l_3+l_4ϵ_{ud}^2l_7)\pi _0^2]`$
with $`l_i`$ the original parameters of Gasser & Leutwyler $`SU(2)`$ lagrangian

$`l_4`$ $`=`$ $`8\alpha _4`$
$`l_3`$ $`=`$ $`16\alpha _3+8\alpha _88\alpha _4`$
$`l_7`$ $`=`$ $`16\alpha _78\alpha _8`$ (10)
where the subindexes ($`i,j`$) in the lagrangian denote the order in powers of momentum and fields, respectively, and

$$_{I\pm }\pm i\mu _Iu.$$

This definition of the covariant derivative is natural, since we know weldon ; actor  that the chemical potential is introduced as the zero component of an external “gauge” field. In the previous expression,

$$|\pi |^2\pi ^+\pi ^{},|_I\pi |^2=(_I\pi )_+(_I\pi )_{}.$$

We will neglect $`ϵ_{ud}^2=(m_um_d)^2/(m_u+m_d)^2`$ because it only shifts in a small quantity the neutral pion mass and we are interested in the thermal and density evolution of the masses.

For renormalizing with counterterms we introduce the following decomposition

$`_{eff}`$ $`=`$ $`_{2,2}+_{2,4}^r+_{4,2}^r`$
$`_{2,2}`$ $`=`$ $`_{2,2}^r+\delta ,`$ (11)
where the $`r`$ index denote the lagrangian with renormalized fields.

Setting $`\pi _0=\sqrt{Z_0}\pi _0^r`$ and $`\pi _\pm =\sqrt{Z_\pm }\pi _\pm ^r`$ in $`_{2,2}`$, we have

$`\delta `$ $`=`$ $`{\displaystyle \frac{1}{2}}\delta _{Z_o}\left[(\pi _0^r)^2m^2(\pi _0^r)^2\right]`$ (12)
$`+\delta _{Z_\pm }\left[\left|_I\pi ^r\right|^2m^2\left|\pi ^r\right|^2\right]`$
with $`\delta _{Z_i}=Z_i1`$.

First, let us consider the temperature and density corrections to the pion propagator. Since our calculation will be at the one loop level, we do not need the full formalism of thermo field dynamics, including thermal ghosts and matrix propagators. The propagator

$$D_\pm (x)=D(x;\pm \mu _I)+D_\beta (x;\pm \mu _I)$$
(13)
for charged pions at the tree level will be given by an extension, for a non-vanishing chemical potential, of the well known Dolan-Jackiw propagators for scalar fields weldon . Note that since there is no chemical potential associated to the neutral pion, the thermal propagator $`D_0`$ will be the usual one

$`D_0(x)=D(x;0)+D_\beta (x;0)`$ (14)
where, in momentum space

$`\underset{¯}{\mathrm{D}}(k;\pm \mu _I)`$ $`=`$ $`{\displaystyle \frac{i}{k_\pm ^2m^2+iϵ}}`$
$`\underset{¯}{\mathrm{D}}_\beta (k;\pm \mu _I)`$ $`=`$ $`2\pi n_B(|ku|)\delta (k_\pm ^2m^2)`$ (15)
with

$$k_\pm k\mu _Iu,n_B(x)=\frac{1}{e^{\beta x}1}$$

are the shifted momentum and the Bose-Einstein factor.

We will use the $`\overline{MS}`$-scheme, and we renormalize as usual at $`T=0`$, since the thermal corrections are finite. The self energy for charged and neutral pions including the counterterms has the form

$`\mathrm{\Sigma }_\pm (p)`$ $`=`$ $`[A_\pm \delta _{Z_\pm }]p_\pm ^2[A_\pm ^{}\delta _{Z_\pm }]m^2+A_\pm ^{\prime \prime }up_\pm `$
$`\mathrm{\Sigma }_0(p)`$ $`=`$ $`[A_0\delta _{Z_0}]p^2[A_0^{}\delta _{Z_0}]m^2`$ (16)
with

$`A_\pm `$ $`=`$ $`{\displaystyle \frac{1}{3f^2}}\left[D_\pm (0)+D_0(0)\right]2{\displaystyle \frac{m^2}{f^2}}l_4`$
$`A_\pm ^{}`$ $`=`$ $`{\displaystyle \frac{1}{6f^2}}\left[2D_\pm (0)D_0(0)\right]2{\displaystyle \frac{m^2}{f^2}}(l_3+l_4)`$
$`A_\pm ^{\prime \prime }`$ $`=`$ $`{\displaystyle \frac{2}{f^2}}uD_\pm (0)`$
$`A_0`$ $`=`$ $`{\displaystyle \frac{2}{3f^2}}D_\pm (0){\displaystyle \frac{m^2}{f^2}}l_4`$
$`A_0^{}`$ $`=`$ $`{\displaystyle \frac{1}{6f^2}}\left[3D_0(0)2D_\pm (0)\right]{\displaystyle \frac{m^2}{f^2}}(l_3+l_4)`$ (17)
Our prescription to fix the counterterm $`\delta _{Z_\pm }`$ is to impose that $`\mathrm{\Sigma }`$ does not depend on $`p^2`$, so, $`\delta _{Z_i}=A_i`$. In this way, the renormalized propagators will take the form

$`i\underset{¯}{\mathrm{D}}_\pm ^r(p)^1`$ $`=`$ $`p_\pm ^2A_\pm ^{\prime \prime }up_\pm m^2[1A_\pm ^{}+A_\pm ]`$
$`i\underset{¯}{\mathrm{D}}_0^r(p)^1`$ $`=`$ $`p^2m^2[1A_0^{}+A_0]`$ (18)
where $`\alpha _i`$ terms absorbs the divergences

$`ł_i`$ $`=`$ $`{\displaystyle \frac{\gamma _i}{32\pi ^2}}\left[\overline{ł}_i+\mathrm{ln}{\displaystyle \frac{4\pi m^2}{\mathrm{\Lambda }^2}}{\displaystyle \frac{2}{d4}}\gamma +1\right]`$ (19)
in which the $`\gamma _i`$ terms are tabulated GL2 ; holstein , being $`\mathrm{\Lambda }`$ a scale factor.

We identify $`m_{\pi ^+}`$ and $`m_\pi ^{}`$ from the solution of $`\underset{¯}{\mathrm{D}}_\pm ^r(p)^1|_{𝐩=0}=0`$ in the frame where the heath bath is at rest ($`u=(1,\mathrm{𝟎})`$). We get the well known result for $`T=\mu _I=0`$

$$m_\pi =m\left(1\alpha _\pi \overline{l}_3/4\right)$$
(20)
is identified with the physical mass. $`\alpha _\pi =(m_\pi /4\pi f_\pi )^2`$ is the perturbative term that fixes the scale of energies in the theory (for energies below $`4\pi f_\pi `$) so we neglect the $`𝒪(g^2)`$ factors. This allows us to set $`mm_\pi `$ in all radiative corrections (and also $`ff_\pi `$). The procedure is the same for $`m_{\pi ^0}`$

It is important to remark that radiative corrections will leave a dependence on the chemical potential for the pion mass only for finite values of temperature. In a strict sense, this procedure does not allow us to say nothing new for an eventual chemical potential dependence of the masses at $`T=0`$ (cold matter) which is already included in $`_2`$. In this case, $`T=0`$, we have to follow the usual procedure, son ; toublan , of computing the minimum of the effective potential in $`_2`$ when the chemical potential is taken into account, without considering radiative corrections. This enables to identify a phase structure where a non trivial vacuum appears for higher values of $`\mu _I`$, $`|\mu _I|>m_\pi `$characterized by the appearance of a condensate $`\pi ^{}`$. (The opposite occurs for negative values of the chemical potential, where the vacuum state is a condensate $`\pi ^+`$). At $`T=0`$ when $`\mu _I=m_\pi `$, the mass of $`\pi ^{}`$ vanishes.

For finite $`T`$ and $`\mu _I`$, we find the following expression for the masses

$`m_{\pi ^\pm }(T,\mu _I)`$ $`=`$ $`m_\pi \left[1+\alpha _\pi I_0\pm \left(\mu _I/m_\pi 4\alpha _\pi J\right)\right]`$
$`m_{\pi ^0}(T,\mu _I)`$ $`=`$ $`m_\pi \left[1+\alpha _\pi \left(2II_0\right)\right]`$ (21)
with

$`I`$ $`=`$ $`{\displaystyle _1^{\mathrm{}}}𝑑x\sqrt{x^21}[n_B(m_\pi x\mu _I)+n_B(m_\pi x+\mu _I)]`$
$`J`$ $`=`$ $`{\displaystyle _1^{\mathrm{}}}𝑑xx\sqrt{x^21}[n_B(m_\pi x\mu _I)n_B(m_\pi x+\mu _I)]`$
$`I_0`$ $`=`$ $`I(\mu _I=0)`$ (22)
Note that our convention for the chemical potential sign is contrary to the one adopted in the paper by Kogut and Toublan toublan , who extended previous results by Son and Stephanov son .

If the chemical potential of the charged pions vanishes, i.e for symmetric matter, at finite T we get the well known result for $`m_\pi (T)`$ due to chiral perturbation theory GL , see also LCL . However, due to radiative corrections to the neutral pion propagator, its mass will acquire a non trivial chemical potential dependence for finite values of temperature. In the approach where the minimum of the effective potential is calculated (for finite $`\mu _I`$ and $`T=0`$), the mass of the neutral pion remains constant.

We show in Fig.1 a tridimensional picture for the behavior of the mass of the neutral pion. Note that when $`\mu _I=0`$, $`m_{\pi _0}(T)=m_{\pi ^\pm }(T)`$.

From Fig.2 we see that at zero temperature, we agree with the usual prediction, $`m_\pi ^+=m_\pi +\mu _I`$. In fact, at zero temperature the $`\pi ^+`$ should condensate when $`\mu _I=m_\pi `$ (the inverse situation occurs for $`\pi ^{}`$). Now, this situation changes if temperature starts to grow. The condensation point disappear at $`\mu _I=m_\pi `$; in $`\mu _I=m_\pi `$ the mass start to decrease. For small $`T`$ (for example inside an neutron star), this effect is neglegible.

In connection with the behavior of $`f_\pi (T,\mu _I)`$ when $`\mu _I<m_\pi `$, we have make used of PCAC, which provides us with a relation between the renormalized propagator and the pion decay constant.

The Axial current is obtained as the functional derivative of the action with respect to $`a_\mu ^a`$, with $`a_\mu =a_\mu ^a\tau ^a/2`$

$$A_\mu ^a=\frac{\delta 𝒮}{\delta a_\mu ^a}(M,0,\mu u_\mu ,0)$$
(23)
The axial current is

$`A_{(1,1)\mu }^\pm `$ $`=`$ $`f(_\mu ^I\pi )^\pm `$
$`A_{(1,3)\mu }^\pm `$ $`=`$ $`{\displaystyle \frac{2}{3f}}\{\pi ^0[\pi ^0(_\mu ^I\pi )^\pm \pi ^\pm _\mu \pi ^0]`$
$`+\pi ^\pm [\pi ^{}(_\mu ^I\pi )^\pm \pi ^\pm (_\mu ^I\pi )^{}]\}`$
$`A_{(3,1)\mu }^\pm `$ $`=`$ $`{\displaystyle \frac{m^2}{f}}2l_4(_\mu ^I\pi )^\pm `$
$`A_{(1,1)\mu }^0`$ $`=`$ $`f_\mu \pi ^0`$
$`A_{(1,3)\mu }^0`$ $`=`$ $`{\displaystyle \frac{2}{3f}}\left\{2|\pi |^2_\mu \pi ^0\pi ^0_\mu |\pi |^2\right\}`$
$`A_{(3,1)\mu }^0`$ $`=`$ $`{\displaystyle \frac{m^2}{f}}2l_4_\mu \pi ^0`$ (24)
Now, the effective axial current at $`𝒪(p^3)`$ will be

$`A_\mu ^i`$ $``$ $`A_{(1,1)\mu }^i+A_{(1,3)\mu }^{ir}+A_{(3,1)\mu }^{ir}`$ (25)
$`=`$ $`\sqrt{Z_i}A_{(1,1)\mu }^{ir}+A_{(1,3)\mu }^{ir}+A_{(3,1)\mu }^{ir}`$
with $`i=\{\pm ,0\}`$. We will take

$$\sqrt{Z_i}=\sqrt{1+\delta _{Z_i}}1+\frac{1}{2}\delta _{Z_i}+𝒪(\delta _{Z_i}^2).$$

The value of the $`\delta _{Z_i}`$ are the same as those obtained in the mass renormalization.

After taking into account the different tadpole diagrams which correct the coupling of the current to one pion states, we find

$`0|A_\mu ^\pm |\pi ^{}(p)`$ $`=`$ $`ip_{\pm \mu }\left[f_\pi f2\alpha _\pi (I+I_0)\right]`$
$`\pm iu_\mu f\alpha _\pi 8J`$
$`0|A_\mu ^0|\pi ^0(p)`$ $`=`$ $`ip_\mu \left[f_\pi f4\alpha _\pi I\right]`$ (26)
with

$`f_\pi `$ $`=`$ $`f\left(1+\alpha _\pi \overline{l}_4\right)`$ (27)
Now, we can set $`ff_\pi `$, $`mm_\pi `$ in all $`𝒪(\alpha _\pi )`$ terms, since any correction will be of order $`\alpha _\pi ^2`$ (including $`\alpha _\pi `$), then we define the effective decay constant as the part proportional to $`p_\mu `$, so

$`f_{\pi ^\pm }(T,\mu _I)`$ $``$ $`f_\pi [12\alpha _\pi (I+I_0)]`$
$`f_{\pi ^0}(T,\mu _I)`$ $``$ $`f_\pi [14\alpha _\pi I]`$ (28)
For $`\mu _I=0`$ we agree with the well known results of Gasser and Leutwyler GL . For an increasing finite chemical potential, the $`f_\pi (T)`$ couplings decrease faster. This effect is enhanced for $`f_{\pi _0}(T)`$ and is related to the fact that $`f_{\pi _0}(T)`$only receives radiative corrections from charged pion tadpoles.

In heavy ion collisions, a finite value of $`\mu _I`$ means that, at least locally, we would expect more pions with definite charge than in the symmetric case. According to this picture, the production rate of dileptons from pion annhilation should be supressed. Probably, the detection of such kind of effects will demand a higher center of mass energy.

In order to explore the region where $`|\mu _I|>m_\pi `$, associated to a new phase where the condensates occur, we need to redefine our fields as fluctuations around the configuration corresponding to a minima of the effective potential in $`_2`$. At present we are working on it, but it is possible to extrapolate, for $`Tm_\pi `$ and $`\mu _Im_\pi `$ the condensation point in such a way that we actually remain in the first phase. However the curve in the $`\mu _IT`$ plane that separates both phases is only reliable in the parameters region mentioned before where in the thermal factors in eq.(22), we have taken the approximation $`n_B(m_\pi x\pm \mu _I)exp[\beta (m_\pi x\pm \mu _I)]`$. A complete analysis of the phase can be found in splittorff2 . The phase diagram is shown in Fig. 5 in accordance with son . However, for higher values of $`\mu _I`$ changes abruptly and our approximation is no longer valid.

Acknowledgements: The work of M.L. has been supported by Fondecyt (Chile) under grant No.1010976. C.V. acknowledges support from a Conicyt Ph.D fellowship (Beca Apoyo Tesis Doctoral).


text:
# Superfluid transition in quasi2D Fermi gases

## Abstract

We show that atomic Fermi gases in quasi2D geometries are promising for achieving superfluidity. In the regime of BCS pairing for weak attraction, we calculate the critical temperature $`T_c`$ and analyze possibilities of increasing the ratio of $`T_c`$ to the Fermi energy. In the opposite limit, where a strong coupling leads to the formation of weakly bound quasi2D dimers, we find that their Bose-Einstein condensate will be stable on a long time scale.

Recent progress in trapping and cooling of Fermi isotopes of K  and Li  has shown the ability to go far below the temperature of quantum degeneracy and to manipulate independently the trapping geometry, density, temperature and interparticle interaction. The Duke experiment  presents intriguing results on the possibility of achieving a superfluid phase transition in the two-component Fermi gas of 6Li.

Two-dimensional Fermi gases have striking features not encountered in 3D. In the superfluid state, thermal fluctuations of the phase of the order parameter strongly modify the phase coherence properties. The interaction strength depends logarithmically on the relative energy of the colliding atoms. For degenerate Fermi gases this energy is of the order of the Fermi energy $`\epsilon _F`$ which is proportional to the 2D density $`n`$. Accordingly, the exponential dependence of the BCS transition temperature on the interaction strength transforms into a power law dependence on the density: $`T_cn^{1/2}`$ . This suggests a unique possibility to cross the critical point by adiabatically expanding a degenerate Fermi gas. Since the ratio $`T/\epsilon _F`$ remains unchanged, the temperature scales as $`n`$ and decreases with density faster than $`T_c`$.

Experimentally it is possible to achieve the quasi2D regime by confining the atoms in one direction so tightly that the corresponding level spacing exceeds the Fermi energy. Under this condition the degenerate Fermi gas is kinematically two-dimensional. Thus far, this regime has been reached for Cs atoms  and for Bose-Einstein condensates of Na  and Rb .

In the quasi2D regime the mean-field interaction between particles exhibits a similar logarithmic dependence on the particle energy as in the purely 2D case . The amplitude of the $`s`$-wave scattering turns out to be sensitive to the strength of the tight confinement . This opens new handles on manipulations of the interparticle interaction and superfluid pairing.

In this Rapid Communication we show that atomic Fermi gases in quasi2D geometries can become strong competitors of 3D gases in achieving superfluidity. The ability to increase the interparticle interaction by tuning the trap frequencies gives an opportunity to realize a transition from the standard BCS pairing in the case of weak attraction to the limit of strong interactions and pairing in coordinate space. In the latter case one eventually gets a dilute system of weakly bound quasi2D dimers of fermionic atoms, which can undergo Bose-Einstein condensation. For the BCS case, we calculate the critical temperature $`T_c`$ to second order in perturbation theory and discuss possibilities of increasing the ratio $`T_c/\epsilon _F`$. In the other extreme, we find that the interaction between the quasi2D dimers is repulsive, and their collisional relaxation and decay are strongly suppressed. This allows us to conclude that BEC of these composite bosons will be stable on a long time scale.

We consider an ultracold two-component Fermi gas in the quasi2D regime and confine ourselves to the $`s`$-wave interaction and superfluid pairing between atoms of different components. We assume that the characteristic radius of the interaction potential is much smaller than the harmonic oscillator length in the tightly confined direction, $`l_0=(\mathrm{}/m\omega _0)^{1/2}`$, where $`m`$ is the atom mass, and $`\omega _0`$ is the confinement frequency. Then the interaction problem involves two length scales: $`l_0`$ and the 3D scattering length $`a`$. For $`a<0`$ and $`|a|l_0`$, there is a peculiar quasi2D weakly bound $`s`$-state of two particles, with the binding energy

$$\epsilon _0=0.915(\mathrm{}\omega _0/\pi )\mathrm{exp}(\sqrt{2\pi }l_0/|a|)\mathrm{}\omega _0.$$
(1)
In this case the coupling constant for the intercomponent interaction takes the form $`g=(4\pi \mathrm{}^2/m)\mathrm{ln}^1(\epsilon _0/\epsilon )`$, where the relative collision energy $`\epsilon `$ is assumed to be either much smaller or much larger than $`\epsilon _0`$ (see  and refs. therein). As in degenerate Fermi gases one has $`\epsilon \epsilon _F`$, the interaction is attractive ($`g<0`$) if the density is sufficiently high and one satisfies the inequality

$$\epsilon _0/\epsilon _F1.$$
(2)
Thus, the inequality (2) is the necessary condition for the BCS pairing. For finding the critical temperature $`T_c`$ below which the formation of Cooper pairs becomes favorable, we go beyond the simple BCS approach and proceed along the lines of the theory developed by Gor’kov and Melik-Barkhudarov for the 3D case .

The critical temperature $`T_c`$ is determined as the highest temperature for which the linearized equation for the order parameter (gap) $`\mathrm{\Delta }=g\widehat{\mathrm{\Psi }}\widehat{\mathrm{\Psi }}`$ has a nontrivial solution . Assuming that the quasi2D gas is uniform in two in-plane directions, the gap equation in the momentum space takes the 2D form

$`\mathrm{\Delta }(𝐪)`$ $``$ $`{\displaystyle }\{{\displaystyle \frac{\mathrm{}^2}{m}}f(z,𝐪,𝐪^{})[K(q^{})+{\displaystyle \frac{1}{z\xi (q^{})+i0}}]`$ (3)
$`+`$ $`\delta V(𝐪,𝐪^{})K(q^{})\}\mathrm{\Delta }(𝐪^{}){\displaystyle \frac{\mathrm{d}^2q^{}}{(2\pi )^2}},`$ (4)
where $`K(q)=(1/2\xi (q))\mathrm{tanh}(\xi (q)/2T)`$, $`\xi (q)=\mathrm{}^2q^2/2m\mu `$, and $`\mu \epsilon _F=\pi \mathrm{}^2n/m`$ is the chemical potential. The first term in the rhs of Eq.(3) results from the direct interaction between particles, and we renormalized the interaction potential in terms of the scattering amplitude $`f`$ (vertex function). The latter is a solution of the quasi2D scattering problem. The parameter $`z`$ has a meaning of the total energy of colliding particles in their center of mass reference frame. It is of the order of $`\epsilon _F`$ and drops out of the final answer. The term $`\delta V(𝐪,𝐪^{})`$ describes the modification of the interparticle interaction due to the presence of other particles (many-body effects). The leading contributions to this term are second order in the scattering amplitude. They are shown in Fig.1 and correspond to an indirect interaction between two particles when one of them interacts with a particle-hole pair virtually created from the ground state (filled Fermi sea) by the other particle. These second order contributions are important for the absolute value of the critical temperature (preexponential factor in the 3D case), whereas higher order terms involving more interaction events can be neglected.

The amplitude $`f`$ is independent of the momenta $`𝐪,𝐪^{}`$ and we will use the 2D relation (see )

$$f=4\pi \mathrm{ln}^1\{\epsilon _0/(z)\}.$$
(5)
The last term in the rhs of Eq.(3) is a small correction, since $`\delta Vf^2`$. Accordingly, the quantity $`|f|`$ represents a small parameter of the theory. In the quasi2D regime ($`z\mathrm{}\omega _0`$) the motion of particles in the tightly confined direction provides a correction to Eq.(5), which is $`(z/\mathrm{}\omega _0)f^2`$ . It is much smaller than $`\delta V`$ and will be omitted.

Equations (3) and (5) show that the momentum dependence of the order parameter appears only due to the second order term that contains many-body contributions to the interparticle interaction $`\delta V`$. The latter is a function of $`p=|𝐪+𝐪^{}|`$ and rapidly decays for $`p>2q_F`$, where $`q_F=\sqrt{2m\mu }/\mathrm{}`$ is the Fermi momentum. For $`p2q_F`$ the quantity $`\delta V`$ is almost constant. Therefore, one has $`\mathrm{\Delta }(𝐪^{})\mathrm{\Delta }(q_F)`$ in a wide momentum range near the Fermi surface. Then, for $`q=q_F`$ a direct integration of Eq.(3) yields

$`\mathrm{\Delta }(q_F)`$ $`=`$ $`{\displaystyle \frac{f(z)}{4\pi }}\mathrm{\Delta }(q_F)\mathrm{ln}\left[{\displaystyle \frac{(2\mu z)}{\pi ^2T^2}}\mathrm{exp}(2\gamma )\right]`$ (6)
$``$ $`\delta V(q_F,q_F){\displaystyle \frac{m}{2\pi \mathrm{}^2}}\mathrm{ln}\left(C{\displaystyle \frac{\mu }{T}}\right)\mathrm{\Delta }(q_F),`$ (7)
where $`\gamma 0.5772`$ is the Euler constant, and $`C`$ is a numerical factor determined by the momentum dependence of $`\mathrm{\Delta }`$ and $`\delta V`$. The calculation of $`\delta V(q_F,q_F)`$ is straightforward and gives $`\delta V(q_F,q_F)=(\mathrm{}^2/2\pi m)f^2(2\mu )`$. Then, using Eq.(5) we obtain the critical temperature

$`T_c(2\mu /\pi )\mathrm{exp}\left(\gamma 1|2\pi \mathrm{Re}f^1(2\mu )|\right).`$

The exponent in this equation should be large and the quantity $`\mathrm{Re}f^1(2\mu )`$ should be negative. As the chemical potential is $`\mu \epsilon _F`$, from Eq.(5) one sees that these requirements are reached under the condition (2). Using Eqs. (1) and (5) the critical temperature takes the form

$$T_c=\frac{\gamma \sqrt{2\epsilon _0\epsilon _F}}{\pi e}=0.16\sqrt{\epsilon _F\mathrm{}\omega _0}\mathrm{exp}\left(\sqrt{\frac{\pi }{2}}\frac{l_0}{|a|}\right)\epsilon _F.$$
(8)
The relative correction to this result is of the order of $`1/|\mathrm{ln}(\epsilon _0/\epsilon _F)|1`$.

Note that Eq. (8) predicts by a factor of $`e`$ smaller value for the critical temperature than a simple BCS calculation . This means that the attractive interaction between particles becomes weaker once we take into account the polarization of the medium.

The ratio $`T_c/\epsilon _F`$ is not necessarily very small. For example, using Feshbach resonances the scattering length is tunable over a wide interval of negative values . Keeping the exponential term equal to $`0.05`$ in Eq.(8), with $`\omega _0100`$ kHz we obtain $`T_c/\epsilon _F0.1`$ for 2D densities $`n10^9`$ cm-2 ($`T_c40`$ nK).

As in the purely 2D case , the transition temperature $`T_cn^{1/2}`$ and the ratio $`T_c/\epsilon _F`$ increases with decreasing density as $`n^{1/2}`$. This is a striking difference from the 3D case, where this ratio decreases exponentially with density. In the presence of the in-plane confinement, one can approach the BCS transition in a degenerate Fermi gas by adiabatically expanding the quasi2D trap in the in-plane direction(s). As the degeneracy parameter $`T/\epsilon _F`$ is conserved in the course of the adiabatic expansion, the ratio $`T/T_c`$ will decrease as $`n^{1/2}`$. Equations (1) and (8) also show that one can increase $`\epsilon _0`$ and $`T_c/\epsilon _0`$ by tuning $`|a|`$ to larger values or by making the tight confinement stronger and thus decreasing $`l_0`$.

What happens if $`\epsilon _0`$ and $`\epsilon _F`$ become comparable with each other, i.e. one reaches the quasi2D resonance for two-body collisions? Then Eq.(8) leads to $`T_c\epsilon _F`$ and is no longer valid. In fact, for $`\epsilon _0>\epsilon _F`$ the formation of bound quasi2D dimers of distinguishable fermions becomes energetically favorable and one encounters the problem of Bose-Einstein condensation of these bosonic molecules. Thus, an increase of the ratio $`\epsilon _0/\epsilon _F`$ from small to large values is expected to provide a transformation of the BCS pairing to molecular BEC. This type of crossover has been discussed in literature in the context of superconductivity  and in relation to superfluidity in 2D films of 3He . The idea of using a Feshbach resonance for achieving a superfluid transition in the BCS-BEC crossover regime in ultracold 3D Fermi gases has been proposed in refs. .

We will not consider the crossover regime and confine ourselves to the limiting case of molecular BEC ($`\epsilon _0\epsilon _F`$). A subtle question is related to the stability of the expected condensate, which depends on the interaction between the molecules. For the repulsive interaction one will have a stable molecular BEC, and the attractive interaction should cause a collapse.

The molecule-molecule scattering is a 4-body problem described by the Schrödinger equation

$`[{\displaystyle \frac{\mathrm{}^2}{m}}(_{𝐫_1}^2+_{𝐫_2}^2){\displaystyle \frac{\mathrm{}^2}{2m}}_𝐑^2+U(r_1)+U(r_2)`$ (9)
$`+{\displaystyle \underset{\pm }{}}U({\displaystyle \frac{𝐫_1+𝐫_2}{2}}\pm 𝐑)E]\mathrm{\Psi }(𝐫_1,𝐫_2,𝐑)=0.`$ (10)
Here $`𝐫_1`$ is the distance between two given distinguishable fermions, $`𝐫_2`$ is the distance between the other two, $`𝐑`$ is the distance between the centers of mass of these pairs, and $`U`$ is the interatomic potential. The total energy is $`E=2\epsilon _0+\epsilon `$, with $`\epsilon `$ being the relative molecule-molecule kinetic energy.

The interaction between molecules is present only at intermolecular distances of the order of or smaller than the size of a molecule $`d_{}=\mathrm{}/\sqrt{m\epsilon _0}`$. Therefore, at energies $`\epsilon \epsilon _0`$ the scattering between molecules is dominated by the $`s`$-wave channel and can be analyzed on the basis of the solution of Eq.(9) for $`\epsilon =0`$. For large $`R`$ the corresponding wavefunction is $`\mathrm{\Psi }(𝐫_1,𝐫_2,𝐑)K_0(r_1/d_{})K_0(r_2/d_{})\mathrm{ln}(\alpha R/d_{})`$, where the decaying Bessel function $`K_0(r_{1,2}/d_{})`$ represents the 2-body bound state. The parameter $`\alpha `$ is a universal constant which can be found by matching the quantity $`\mathrm{ln}(\alpha R/d_{})`$ with the solution of Eq.(9) at short distances. Finally, matching $`\mathrm{ln}(\alpha R/d_{})`$ with the wavefunction of free relative motion of two molecules at distances $`d_{}R\mathrm{\Lambda }_\epsilon `$, where $`\mathrm{\Lambda }_\epsilon =\mathrm{}/\sqrt{m\epsilon }`$ is their de Broglie wavelength, we obtain the coupling constant (scattering amplitude) for the interaction between molecules:

$$g_m=(2\pi \mathrm{}^2/m)\mathrm{ln}^1(2\alpha ^2e^{2\gamma }\epsilon _0/\epsilon )>0;\epsilon \epsilon _0.$$
(11)
A precise value of $`\alpha `$ is not important as it gives rise to higher order corrections in Eq.(11). However, in order to make sure that this constant is neither anomalously large nor anomalously small we have integrated Eq.(9) numerically. For this purpose, it is convenient to transform Eq.(9) into an integral equation for a function which depends only on three independent coordinates. This has been done by using the method of ref. . Our calculations lead to $`\alpha 1.6`$. They show the absence of 4-body weakly bound states and confirm an intuitive picture that the interaction between two molecules can be qualitatively represented by means of a purely repulsive potential with the range $`d_{}`$. For the interaction between Bose-condensed dimers, in Eq.(11) one has $`\epsilon =2n_mg_m\epsilon _0`$, where $`n_m`$ is the density of the dimers (see  and refs. therein). We thus conclude that a Bose condensate of these weakly bound dimers is stable with respect to collapse.

The 2D gas of bosons becomes Bose-condensed below the Kosterlitz-Thouless transition temperature $`T_{KT}`$  which depends on the interaction between particles. According to the recent quantum Monte Carlo simulations , for the 2D gas with the coupling constant (11) the Kosterlitz-Thouless temperature is given by

$$T_{KT}=(\pi \mathrm{}^2n_m/m)\mathrm{ln}^1\left[(\eta /4\pi )\mathrm{ln}\left(1/n_md_{}^2\right)\right],$$
(12)
where the numerical factor $`\eta 380`$. For $`\epsilon _F\epsilon _0`$, the density of dimers $`n_mn/2`$ and the parameter $`(1/n_md_{}^2)2\pi \epsilon _0/\epsilon _F`$. Then, for $`\epsilon _0/\epsilon _F=10`$, Eq.(12) gives $`T_{KT}/\epsilon _F0.1`$ and at densities $`10^8`$ cm-2 the transition temperature in the case of 6Li is $`T_{KT}30`$ nK.

The weakly bound dimers that we are considering are molecules in the highest rovibrational state and they can undergo collisional relaxation and decay. The relaxation process occurs in pair dimer-dimer or dimer-atom collisions. It produces diatomic molecules in deep bound states and is accompanied by a release of the kinetic energy. The size of these deeply bound molecules is of the order of the characteristic radius of the interatomic potential $`R_el_0`$, and their internal properties are not influenced by the tight confinement. Therefore, the relaxation can be treated as a 3D process and it requires the presence of at least three fermionic atoms at distances $`R_e`$ between them. Since at least two of them are identical, the relaxation probability acquires a small factor $`(kR_e)^2`$ compared to the case of bosons, where $`k1/d_{}=\sqrt{\epsilon _0/\mathrm{}\omega _0}/l_0`$ is a characteristic momentum of atoms. The 3D density of atoms in the quasi2D geometry is $`n/l_0`$. Thus, qualitatively, the inverse relaxation time can be written as $`\tau _{\mathrm{rel}}^1\alpha _{\mathrm{rel}}n(R_e/l_0)^2(\epsilon _0/\mathrm{}\omega _0)/l_0`$, where $`\alpha _{\mathrm{rel}}`$ is the relaxation rate constant for the highest rovibrational states of 3D molecules of two bosonic atoms. We estimate $`\tau _{\mathrm{rel}}`$ keeping in mind the recent measurements for Rb2 molecules  which give $`\alpha _{\mathrm{rel}}3\times 10^{11}`$ cm3/s. For $`l_0`$ in the interval from $`10^5`$ to $`10^4`$ cm, the suppression factor $`(R_e/l_0)^2(\epsilon _0/\mathrm{}\omega _0)`$ ranges from $`10^3`$ to $`10^5`$ and at 2D densities $`n10^8`$ cm-2 we find the relaxation time $`\tau _{\mathrm{rel}}`$ of the order of a second or larger.

Dimer-dimer pair collisions can lead to the formation of bound trimers, accompanied by a release of one of the atoms. The formation of deeply bound trimer states will be suppressed at least in the same way as the relaxation process discussed above. Therefore, it is important that there are no weakly bound trimers in (quasi)2D. Just as in 3D , this can be established by using the zero-range approach ($`R_e0`$). We have performed this analysis along the lines of the 3D work . Qualitatively, the symmetry of the 3-fermion system containing two identical fermions provides a strong centrifugal repulsion that does not allow the presence of 3-body bound states. This is in contrast to 2D bosons where one has two fully symmetric trimer bound states .

Thus, the life-time of quasi2D dimers of fermionic atoms is rather long and one easily estimates that it greatly exceeds the characteristic time of elastic collisions. One can even think of achieving BEC in the initially non-condensed gas of dimers produced out of a non-superfluid atomic Fermi gas under a decrease of $`n`$ or $`l_0`$.

In conclusion, we have found the temperature of superfluid phase transition in two-component quasi2D Fermi gases. Our results are promising for achieving this transition in both the regime of BCS pairing and the regime of BEC of weakly bound dimers.

We acknowledge fruitful discussions with C. Salomon, A. Mosk, C. Lobo, and J.T.M. Walraven. This work was financially supported by the the Dutch Foundations NWO and FOM, by Deutsche Forschungsgemeinschaft (DFG), by INTAS, and by the Russian Foundation for Basic Research. Le LKB est UMR 8552 du CNRS, de l’ENS et de l’Université P. et M. Curie.