format: text

text:
# The Sunyaev-Zel’dovich Effect Spectrum of Abell 2163

## 1 INTRODUCTION

The thermal Sunyaev-Zel’dovich effect (SZE) is a spectral distortion in the Cosmic Microwave Background (CMB) which results from the inverse Compton scattering of low-energy CMB photons by thermally distributed hot electrons in massive clusters of galaxies. The thermal SZE manifests itself as an increment in the Wien part of the CMB spectrum and a decrement in the Rayleigh-Jeans part of the spectrum (Sunyaev & Zel’dovich, 1972), with a crossover wavelength of $`1.38`$ mm. There is also a kinetic SZE, caused by the bulk radial motion of the cluster with respect to the CMB rest frame. This effect shifts the CMB Planck spectrum to a slightly lower (higher) temperature for positive (negative) peculiar velocities where positive is defined as motion away from the observer.

In this paper, we present a recent 1 cm interferometric detection of the SZE in the cluster Abell 2163, obtained with the BIMA and OVRO arrays. We combine this measurement with a previous result from the Sunyaev-Zel’dovich Infrared Experiment (SuZIE), a six element bolometer array (Holzapfel et al. 1997a, hereafter H1997a) on the 10.4 m Caltech Submillimeter Observatory. The SZE in Abell 2163 was also observed by Diabolo, a two element bolometer array mounted on the IRAM 30 m telescope (Desert et al., 1998), and their result is included in Table 1. The measurement has considerable uncertainty and carries very little weight in the spectral fit relative to the SuZIE data. Due to the complication of estimating the effects of systematics on the Diabolo measurement, we have chosen to exclude it from our spectral fit. Thus, combining the four measurements from SuZIE and OVRO/BIMA, we construct the SZE spectrum and use it to determine the central Compton $`y`$-parameter and constrain the radial peculiar velocity of the cluster. We also address the concern that Galactic dust along the line of sight to the cluster may be a significant contaminant in the SuZIE bands (Lamarre et al., 1998). Using data from the Infrared Space Observatory (ISO) and the Infrared Astronomical Satellite (IRAS) we determine the level of dust contamination in the SuZIE measurements.

In §2, we describe our observations, data reduction, and dust analysis. In §3 we describe the theoretical model for the SZE spectrum and present the results of fitting our data to this model; we also consider several systematic uncertainties affecting the fit. All uncertainties are quoted at the 68% confidence level unless otherwise specified.

## 2 OBSERVATIONS AND DATA REDUCTION

The 1 cm observations were carried out at the Berkeley Illinois Maryland Association (BIMA) millimeter array and the Owens Valley Radio Observatory (OVRO) millimeter array during the summers of 1996, 1997, and 1998. The telescopes were outfitted with centimeter-wavelength receivers specifically designed to detect the SZE in distant galaxy clusters (Carlstrom et al., 1996). The receivers use cooled high electron mobility transistor (HEMT) amplifiers (Pospieszalski et al., 1995) sensitive to radiation between 0.83 and 1.15 cm with typical receiver temperatures in the range 11–20 K. Typical system temperatures (scaled to above the atmosphere) for the Abell 2163 observations range between 45 and 55 K.

We observed Abell 2163 at OVRO in 1996 and 1998 for a total of 37 hours with six 10.4 m telescopes providing a $`4\mathrm{}.2`$ FWHM primary beam; we used two 1 GHz wide bands centered at 1.0 cm and 1.05 cm. This cluster was observed at BIMA in 1996, 1997 and 1998 for 11 hours with nine 6.1 m telescopes providing a $`6\mathrm{}.6`$ FWHM primary beam; we used an 800 MHz bandwidth centered at 1.05 cm. At both observatories, we interleaved cluster measurements with observations of the radio point source 1549+026 ($`3.60\pm 0.01`$ Jy) every 25 minutes to monitor system gains, and we used Mars for our absolute flux calibration (Grego et al., 2000). We reduced the data using the MIRIAD software package (Sault et al., 1995) at BIMA and the MMA package (Scoville et al., 1993) at OVRO. In all observations, we flagged data from baselines where one telescope was shadowed by another, data that were not bracketed in time by observation of a phase calibrator, and data showing evidence of poor atmospheric coherence or spurious correlations. Figure 1 shows a deconvolved SZE image of Abell 2163 constructed from BIMA data. SZE images are produced with the software package DIFMAP (Pearson et al., 1994). In Figure 1 we have applied a Gaussian taper with a half power radius of 1.5 k$`\lambda `$ to the $`u\text{-}v`$ data. A bright point source has been removed from the image in Figure 1; see §2.1 for more details.

An interferometer measures the Fourier transform of the sky brightness, hence we fit a model to the data in the Fourier plane where the noise characteristics and spatial filtering of the interferometer are well understood. The model is constructed in the image plane, multiplied by the primary beam of the telescope, then Fourier transformed to the conjugate plane, known as the $`u`$-$`v`$ plane, where it is compared with the data.

We model the cluster gas as a spherical isothermal-$`\beta `$ model (Cavaliere & Fusco-Femiano, 1978), where the SZE intensity change is given by

$$\mathrm{\Delta }I_\nu (\theta )=\mathrm{\Delta }I_{\nu \text{0}}\left(1+\frac{\theta ^2}{\theta _c^2}\right)^{(13\beta )/2}.$$
(1)
Here $`\mathrm{\Delta }I_{\nu \text{0}}`$ is the central SZE decrement/increment at frequency $`\nu `$, $`\theta _c`$ is the angular core radius of the cluster, and $`\beta `$ is a power law index. We adopt the model shape parameters $`\theta _c=1\mathrm{}.20\pm 0\mathrm{}.11`$ and $`\beta =0.616\pm 0.031`$ from an analysis of ROSAT X-ray data (Holzapfel et al. 1997b, hereafter H1997b). With these shape parameters we find $`\mathrm{\Delta }I_{\text{30}}`$$`=0.048\pm 0.006\text{ MJy sr}^1`$ (corresponding to a CMB temperature decrement $`\mathrm{\Delta }T_{\text{CMB}}=1780\pm 210`$ $`\mu `$K) where $`\mathrm{\Delta }I_{\text{30}}`$ is the SZE decrement at 30 GHz. The statistical uncertainty in the central decrement is determined by varying $`\mathrm{\Delta }I_{\text{30}}`$ while keeping the shape parameters fixed; all other parameters are free to assume their best-fit values. The uncertainties in $`\theta _c`$ and $`\beta `$ quoted above lead to an additional $`0.001\text{ MJy sr}^1`$ uncertainty on the 1 cm decrement.

The 1 cm decrement and combined statistical uncertainty are listed in Table 1, which also contains all spectral SZE data for this cluster available to date. The SuZIE $`\mathrm{\Delta }I_{\nu \text{0}}`$ values are from fits to the SuZIE data using the average of the measured spectral responses of the individual channels. The $`\mathrm{\Delta }I_{\nu \text{0}}`$ value from the Diabolo experiment is from Desert et al. (1998). As the table indicates, we have made small corrections in the SuZIE bands to account for dust contamination along the line of sight to the cluster. These corrections are described in more detail in §2.2.

### 2.1 Radio Point Sources

Point sources are identified in DIFMAP using a high-resolution map produced with data from baselines 20 meters and longer. We find one point source in the Abell 2163 data offset $`40\mathrm{}`$ in RA and $`16\mathrm{}`$ in declination from the SZE pointing center, a position consistent with that of a variable radio source with an inverted spectrum detected by the VLA (M. Birkinshaw 2001, private communication) at 2 cm and 6 cm. The point source has been removed from the image in Figure 1.

The flux and position of the detected point source is jointly fit with the SZE decrement, but we must consider the uncertainty in the decrement introduced by point sources in the field with flux densities below our detection threshold. Such a point source near the cluster SZE center would cause us to underestimate the magnitude of the decrement, whereas a point source in a beam sidelobe could lead to an overestimate. To determine the resulting uncertainty we have randomly added point sources to the data as follows. Progressing in steps of 1 $`\mu `$Jy, we choose a Poisson deviate of the expected number of point sources at each flux density. We add sources to the field accordingly out to a radius of $`5\mathrm{}`$, well past the half power point of the BIMA beam where faint point sources should have a negligible effect on the decrement. After reaching our $`3\sigma `$ flux density limit, we refit the SZE model without accounting for the added point sources.

To determine the number of point sources expected at each flux density, we first determine a point source number count versus flux density relationship of the form

$$dn/dSS^\alpha \text{ mJy}^1\text{ arcmin}^2$$
(2)
from our sample of 41 clusters observed at BIMA, where $`n`$ is the number density of point sources and $`S`$ is their flux density. We find that the distribution of sources beyond $`10\mathrm{}`$ from cluster centers agrees well with the expected distribution from a larger sample of point sources at 31 GHz (C. Pryke 2002, private communication). Inside $`10\mathrm{}`$, however, we detect nearly 30 times more point sources than expected due to the strong contribution of cluster cD galaxies. We therefore exclude the inner $`10\mathrm{}`$ of the map when selecting point sources to construct $`dn/dS`$ and address the issue of central point sources in more detail below. We choose $`240\mathrm{}`$ as an outer limit in selecting point sources as it is a good trade-off between encompassing a large number of detected sources and allowing a low beam-corrected flux density threshold. We use only point sources with beam-corrected flux densities greater than 2 mJy, since a 2 mJy point source would have a beam-attenuated flux density of 0.75 mJy at $`240\mathrm{}`$ and 0.75 mJy is a good fiducial $`3\sigma `$ threshold for all 41 maps. Lastly, we assume the power law index determined from the Pryke sample, $`\alpha =2.4`$, in Equation (2). Solving for the normalization, we find

$$dn/dS=2.0\times 10^2\left(\frac{S}{\mathrm{mJy}}\right)^{2.4}\text{mJy}^1\text{ arcmin}^2.$$
(3)
We then randomly distribute point sources according to this empirical relationship, starting from a flux density of 10 $`\mu `$Jy (tests of the Monte-Carlo showed that starting at a lower flux density has a negligible effect on the decrement) up to a flux density of 300 $`\mu `$Jy, corresponding to the $`3\sigma `$ limit of the more sensitive OVRO Abell 2163 data. The resulting distribution of best-fit decrements is nearly Gaussian and centered at the original decrement with a standard deviation of 0.002 MJy sr-1. We adopt this as the statistical uncertainty due to randomly distributed undetected point sources, and it has been combined in quadrature with the other statistical uncertainties in Table 1.

We attempt to estimate the uncertainty from undetected point sources a second way, using information on the distribution of point sources in the field from observations at longer wavelengths. Sources with flux densities greater than $`2`$ mJy at 21 cm appear in the NVSS catalog (Condon et al., 1998). The catalog contains four point sources within $`240\mathrm{}`$ of Abell 2163 that we do not detect in our 1 cm observations. Searching the 100 $`\mu `$Jy rms OVRO map at the positions of these sources, we find no evidence of emission from three of the four; their positions correspond to signals at the $`0.2\sigma `$, $`0.4\sigma `$, and $`1\sigma `$ levels. The fourth, however, displaced $`30\mathrm{}`$ in RA and $`127\mathrm{}`$ in declination from the map center, corresponds to a $`2\sigma `$ signal in the OVRO map. We do not marginalize over the fluxes of all four sources because of the loss of degrees of freedom in the fit, but marginalizing over the $`2\sigma `$ source we find no change in either the decrement or the uncertainty in the decrement, which is dominated by the noise in the map.

To account for the fact that on-center point sources are statistically more common we place a point source model at the estimated position of Abell 2163’s bright central galaxy (RA $`16^\text{h}15^\text{m}49^\text{s}.0`$, Dec $`6^{}8\mathrm{}42\mathrm{}.7`$ from Digitized Sky Survey) and marginalize over its flux. The new best-fit decrement is $`0.049`$ MJy sr-1, an increase in magnitude of just $`0.001\text{ MJy sr}^1`$ from the original decrement. The new decrement falls easily within the 68% uncertainty of the original; we revisit this result when we discuss systematic uncertainties (§3.1).

### 2.2 Analysis of Dust Emission

Lamarre et al. (1998) (hereafter L1998) determined that IR cirrus flux is a non-negligible contaminant in the three SuZIE bands by including submillimeter data from the PRONAOS balloon experiment (Serra, 1998) in their fit to the SuZIE and Diabolo data,. They accounted for dust by fitting to a combined SZE-spectral model and modified blackbody model, and found that removing the dust signal changed $`\mathrm{\Delta }I_{\nu \text{0}}`$ significantly. Specifically, they found that the 1.1 mm signal shrinks by $`0.051\text{ MJy sr}^1`$ to $`0.236\text{ MJy sr}^1`$, the 1.4 mm signal changes by $`0.018\text{ MJy sr}^1`$ to $`0.124\text{ MJy sr}^1`$, and the 2.1 mm band changes by $`0.004\text{ MJy sr}^1`$ to $`0.387\text{ MJy sr}^1`$. L1998 assumed in their analysis that the SuZIE scan pattern is identical to the PRONAOS chopping scheme. However, PRONAOS performed a $`6\mathrm{}.2`$ chop at constant elevation (see L1998 for details) while SuZIE uses a drift scan at constant declination (H1997b). Thus, the experimental filters are different and will not necessarily detect the same dust signal.

To better estimate the level of dust contamination in the high frequency SuZIE channels, we have applied the SuZIE observing scheme to extrapolated dust maps at the SuZIE wavelengths. The expected dust contamination is listed in Table 2 and the dust-corrected SuZIE measurements are given in Table 1. Here we describe the procedure.

For clarity, we summarize the map construction procedures as applied to the 1.1 mm band; the 1.4 mm and 2.1 mm procedures are identical. The 1.1 mm dust map is created by extrapolating from three existing dust maps centered on Abell 2163; ISO maps at 90 $`\mu `$m and 180 $`\mu `$m, and an IRAS map at 100 $`\mu `$m. We extract a $`15\mathrm{}\times 1\mathrm{}.75`$ central region from each dust map, consistent with the area of the SuZIE scans. We then smooth each map to a resolution consistent with the SuZIE beam size ($`1\mathrm{}.75`$), taking into account the intrinsic resolution of the IRAS and ISO experiments. We use a fiducial pixel to pixel uncertainty of $`0.2\text{ MJy sr}^1`$ (cf. Beichman et al. 1988) for the IRAS map; for the ISO maps, we adopt the values estimated in L1998 of $`0.14\text{ MJy sr}^1`$ and $`0.18\text{ MJy sr}^1`$ at 90 and 180 $`\mu `$m, respectively. For each ISO map, we then bin the $`0\mathrm{}.25`$ pixels into $`0\mathrm{}.75`$ pixels, consistent with binning performed in the SuZIE analysis. For each new $`0\mathrm{}.75`$ pixel, we take the average of the original three $`0\mathrm{}.25`$ pixels as the intensity. The IRAS pixels are roughly twice the size of the binned ISO pixels, so we linearly interpolate to estimate the intensity of interleaved pixels, which we artificially add to the IRAS map. We then fit a modified blackbody to these dust maps at each pixel position; since we have three dust maps, each fit is to three data points, or one “pixel triplet”. This allows us to extrapolate the intensities out to 1.1 mm. We fit a modified blackbody of the form

$$F_d\frac{\nu ^{3+\gamma }}{e^{h\nu /kT_d}1},$$
(4)
where $`\gamma `$ is the spectral index, $`T_d`$ is the dust temperature, and we normalize the model with the dust flux at 180 $`\mu `$m, $`F_d(180)`$. Following L1998, we fix the spectral index at $`\gamma =2`$. Upon fitting each dust pixel triplet to this model, we find $`F_d(180)`$ in the range $`25.1`$ to $`26.1`$ MJy sr-1 and $`T_d`$ in the range 21.3 to 21.9 K over all triplets. We then use each triplet’s best fit $`F_d(180)`$ and $`T_d`$ with Gaussian error propagation to extrapolate the intensity at that position and create our 1.1 mm dust map.

After repeating this process at 1.4 and 2.1 mm, we perform the simulated scan following the detailed discussion of SuZIE’s instrumental specifications and observing scheme in H1997b. In our simulation we scan across the $`0\mathrm{}.75`$ pixels of the dust maps in the direction of increasing right ascension (RA). For each pixel along the scan, we record the difference between its intensity and that of the pixel advanced $`4\mathrm{}.6`$ in RA; pixel-to-pixel uncertainties are combined in quadrature. This procedure is equivalent to SuZIE allowing the sky to drift overhead with the earth’s rotation and recording the difference in signal between two bolometer array elements separated by $`4\mathrm{}.6`$ on the sky. We do not treat SuZIE’s $`2\mathrm{}.3`$ “triple beam chop” here as it has little sensitivity to linear changes in brightness across the sky and therefore is less sensitive to dust; the results presented here serve as an upper limit to the dust contribution (see H1997b for more details). Our simulation thus creates “differenced” maps of the dust brightness.

We perform one-dimensional $`\chi ^2`$ fits in the parameter $`\mathrm{\Delta }I_{\nu \text{0}}`$ to the dust maps, fitting to a SuZIE-differenced $`\beta `$ model with the pixel-to-pixel dust noise as our uncertainty. We use a Monte Carlo procedure and successively fit to $`10^6`$ different realizations of the dust; in a single realization, each pixel has added to it some Gaussian-distributed amount of the pixel-to-pixel noise. The results of the Monte Carlo are highly Gaussian as expected; we thus take the most frequently returned $`\mathrm{\Delta }I_{\nu \text{0}}`$ as our best fit and those that enclose $`68.3\%`$ of the returned $`\mathrm{\Delta }I_{\nu \text{0}}`$’s as our $`1\sigma `$ uncertainties. The results are shown in Table 2 and reflected in the last column of Table 1. We find a dust contamination level roughly seven times smaller than that reported in L1998, and small compared to the statistical uncertainties of the SZE measurements.

## 3 SZE SPECTRAL MODEL AND RESULTS

We have fit the dust-corrected SZE spectral data in Table 1 to a model which consists of both the thermal and kinetic SZE components. The thermal component, $`\mathrm{\Delta }I_\text{T}`$, can be written as

$$\mathrm{\Delta }I_\text{T}=I_\text{0}yf(x)(1+\delta _\text{T})$$
(5)
where $`x=h\nu /kT_{\text{CMB}}`$, $`T_{\text{CMB}}=2.728`$ K is the CMB temperature (Fixsen et al., 1996), $`\delta _\text{T}`$ is a relativistic correction to the thermal effect good to fifth order in $`kT_e/m_ec^2`$ (Itoh et al., 1998), $`I_\text{0}2(kT_{\text{CMB}})^3/(hc)^2`$, and

$$f(x)=\frac{x^4e^x}{(e^x1)^2}\left[\frac{x(e^x+1)}{e^x1}4\right].$$
(6)
The quantities $`h`$, $`k`$, and $`c`$ are the Planck constant, Boltzmann constant, and speed of light, respectively. The quantity $`y`$ is the Compton $`y`$-parameter, and is given by

$$y=\left(\frac{kT_e}{m_ec^2}\right)n_e\sigma _T𝑑l,$$
(7)
where $`T_e`$ is the electron temperature of the intracluster medium (ICM), $`m_e`$ is the electron mass, $`n_e`$ is the electron density, and $`\sigma _T`$ is the Thomson scattering cross section. The Compton $`y`$-parameter is proportional to the pressure integrated along the line of sight ($`dl`$) through the cluster; we parameterize it in terms of a central value, which we denote $`y_\text{0}`$.

The kinetic component, $`\mathrm{\Delta }I_\text{K}`$, depends on both $`y`$ and the radial peculiar velocity of the cluster, $`v_\text{p}`$, through the equation

$$\mathrm{\Delta }I_\text{K}=(m_ec^2/kT_e)I_\text{0}y\frac{v_\text{p}}{c}g(x)(1+\delta _\text{K})$$
(8)
where $`g(x)=x^4e^x/(e^x1)^2`$ and $`\delta _\text{K}`$ is a relativistic correction to the kinetic effect good to third order in $`kT_e/m_ec^2`$ (Nozawa et al., 1998). Corrections of order $`v_\text{p}`$$`{}_{}{}^{2}/c^2`$ to the thermal and kinetic SZE are neglected since $`v_\text{p}`$$`c`$. The total intensity is then the sum of the thermal and kinetic components,

$$\mathrm{\Delta }I_\nu =\mathrm{\Delta }I_\text{T}+\mathrm{\Delta }I_\text{K}.$$
(9)
We perform a least squares fit to this model in the parameters $`y_\text{0}`$ and $`v_\text{p}`$ assuming an ICM temperature $`T_e`$ of $`12.4_{1.9}^{+2.8}`$ keV from H1997b.

The resulting $`v_\text{p}`$ and $`y_\text{0}`$ from the least squares fit to the four data points are $`y_\text{0}=3.56\pm 0.41\times 10^4`$ and $`v_\text{p}=410_{850}^{+1030}\text{ km s}^1`$ with $`\chi ^2=0.8`$ for two degrees of freedom. We compare this with the results of our fit to the SuZIE data only, $`y_\text{0}=3.72\pm 0.51\times 10^4`$ and $`v_\text{p}=290_{820}^{+1020}\text{ km s}^1`$. The addition of the 1 cm point improves the $`y_\text{0}`$ constraints by 18%. The constraints on $`v_\text{p}`$ hardly change when the 1 cm point is included since constraints on peculiar velocity are dominated by the SuZIE data points which straddle the thermal SZE null.

Figure 2 shows the SZE spectrum generated from the fit to all four frequency bands. There is excellent consistency amongst all of the measurements. The data were obtained from independent experiments, each with its own telescope, detector, and observing strategy, operating at four different wavelengths. The consistency amongst the data attests to the reliability and accuracy with which the SZE is now being measured.

### 3.1 Systematic Uncertainties

In this section we address several additional sources of uncertainty in our spectral fits; the results of this analysis are summarized in Table 3.

We first consider the systematic uncertainty associated with the dust emission calculated in §2.2. When we alter the dust intensity at the 68% level, the spectral fit results change by $`{}_{0.01}{}^{+0.00}\times 10^4`$ for $`y_\text{0}`$ and $`{}_{30}{}^{+0}\text{ km s}_{}^{1}`$ for $`v_\text{p}`$. The dust uncertainties are sufficiently small as to have very little impact on the spectral fit results. We find a more significant change when we repeat the dust analysis with a smaller spectral index, $`\gamma `$. Finkbeiner et al. (1999) recently suggested based on fits to FIRAS data, that while $`\gamma =2`$ is appropriate for wavelengths $`600`$ $`\mu `$m, the spectral index at longer wavelengths is closer to 1.5. Setting $`\gamma =1.5`$ for the entire dust spectrum we find that $`y_\text{0}`$ increases by $`0.02\times 10^4`$ and $`v_\text{p}`$ decreases by $`70\text{ km s}^1`$. We combine these in quadrature with the changes that resulted from scaling the dust at the 68% level and adopt the total as the systematic uncertainty due to dust contamination.

Temperature fluctuations in the CMB are a potential source of confusion for SZE measurements. Measurements of the kinetic SZE are particularly susceptible since the spectral dependence of the kinetic SZE and CMB fluctuations are identical, absent relativistic corrections to the SZE at the few percent level. Recent BIMA observations at 28.5 GHz have provided good constraints on CMB anisotropy at the angular scales of our measurement (Dawson et al., 2001; Holzapfel et al., 2000); they find $`\mathrm{\Delta }I_\nu <4.8\times 10^4\text{ MJy sr}^1`$ at 95% confidence, which is just 1% of $`\mathrm{\Delta }I_{\text{30}}`$. This uncertainty is negligible given the precision of our measurement. The level of CMB contamination may be larger in the SuZIE bands, however, particularly near the null in the thermal SZE.

To test the level of CMB contamination in the SuZIE bands we first use the software package CMBfast (Zaldarriaga & Seljak, 2000) to generate a CMB power spectrum up to a multipole $`\mathrm{}=2500`$. We assume a universe consisting of cold dark matter with a cosmological constant ($`\mathrm{\Lambda }`$CDM) and use conventional parameter values $`(\mathrm{\Omega }_b,\mathrm{\Omega }_{cdm},\mathrm{\Omega }_\mathrm{\Lambda },\tau _c,n_s,h)=(0.04,0.26,0.70,0,0.95,0.72)`$ where $`\mathrm{\Omega }_b`$ is the fractional baryon density of the universe, $`\mathrm{\Omega }_{cdm}`$ is the fractional cold dark matter density, $`\mathrm{\Omega }_\mathrm{\Lambda }`$ is the fractional density of vacuum energy, $`\tau _s`$ is the optical depth to reionization, $`n_s`$ is the spectral index of the power spectrum, and $`H_0=100h`$ km s-1 Mpc-1 is the Hubble constant. The values chosen are good estimates given the range of values currently constrained by CMB anisotropy measurements (Pryke et al., 2001; Netterfield et al., 2001; Stompor et al., 2001); $`h`$ comes from the Hubble Space Telescope Hubble Constant key project (Freedman et al., 2001). To account for secondary anisotropies we conservatively assume a flat power spectrum for $`\mathrm{}>2500`$, corresponding to arcminute angular scales. The flat band power we use corresponds to $`\mathrm{\Delta }T_{\text{CMB}}/T_{\text{CMB}}5\times 10^6`$, which is consistent with expectations of secondary anisotropies (Hu & Dodelson, 2001) and the BIMA limit of Dawson et al. (2001). We then use this power spectrum to produce mock CMB fields at SuZIE’s resolution. After simulating the SuZIE scan on these maps, we fit them to an isothermal $`\beta `$-model and use a Monte Carlo procedure similar to that used with the dust to determine the contribution of CMB to the total signal in each SuZIE band. Upon fitting to $`10^4`$ different maps for each band, we find an rms of $`\pm 0.020\text{ MJy sr}^1`$ at 1.1 mm, $`\pm 0.022\text{ MJy sr}^1`$ at 1.4 mm, and $`\pm 0.017\text{ MJy sr}^1`$ at 2.1 mm. Further tests using only primary anisotropy $`(\mathrm{}2500)`$ or only secondary anisotropy $`(\mathrm{}2500)`$ indicate that primary anisotropy is responsible for roughly 80% of the total CMB contribution in all three bands. Including the rms total anisotropy contributions in the SuZIE signals listed in Table 1, we redo the spectral fit and find changes of $`\pm 0.05`$ in $`y_\text{0}`$ and $`{}_{400}{}^{+420}\text{ km s}_{}^{1}`$ in $`v_\text{p}`$. The peculiar velocity suffers more from confusion with CMB anisotropy than does the Compton $`y`$-parameter, as anticipated. We adopt these results as the systematic uncertainties due to CMB contamination.

We also consider confusion with extended emission from the radio halo surrounding Abell 2163, which is the brightest such halo yet discovered (Feretti et al., 2001). To test its effect on the 1 cm decrement, we combine the OVRO data with the 1.4 cm NVSS map of the cluster region, which shows a large region of extended emission surrounding the cluster. We first scale the NVSS map to account for the unknown spectrum of the radio halo; we make this scale factor a free parameter and redo the fit described in §2. The decrement is increased by $`0.003\text{ MJy sr}^1`$ when we include this extra parameter, and the scaling is consistent with a spectral index $`\alpha 1.2`$ where the intensity has a $`\nu ^\alpha `$ dependence. We note, however, that our fit is consistent with $`\alpha =\mathrm{}`$ (i.e., no radio halo) within 68% confidence. Furthermore, the spectral index may not be constant across the halo, so the best-fit $`\alpha `$ must be taken as a crude estimate. H1997b also address radio halo confusion, and estimate a 2% effect at 2.1 mm. This assumes a spectral index of $`1.5`$, however; using our best-fit value $`\alpha =1.2`$, the effect increases to $`5\%`$. When we alter the 1 cm and 2.1 mm points accordingly and redo the spectral fit, we find that $`y_\text{0}`$ increases by $`0.17\times 10^4`$ to $`3.73\times 10^4`$, and $`v_\text{p}`$ increases by $`50\text{ km s}^1`$ to $`460\text{ km s}^1`$. We adopt these as estimates of the systematic uncertainties in $`y_\text{0}`$ and $`v_\text{p}`$ due to radio halo confusion.

We use a similar treatment to determine the systematic uncertainty due to point source contamination. In §2 we accounted for the possibility of an undetected on-center point source in the field at 1 cm, and found that it could increase the magnitude of the measured decrement by $`0.001\text{ MJy sr}^1`$. Such a faint source would not affect the SuZIE data unless it had a strongly inverted spectrum; we assume that it does not. The brighter off-center point source detected at 1 cm could affect the SuZIE data if its spectrum were inverted at millimeter wavelengths. Although its flux rises from 1 to 3 mJy between 6 cm and 2 cm (M. Birkinshaw 2001, private communication), we found at the time of our observations that the flux of the source at 1 cm varied between 1.1 and 1.4 mJy. This corresponds to a spectral index $`\alpha `$ between 1.1 and 1.5. Extrapolating with such a spectral index to 2.1 mm, we find a flux of less than 0.2 mJy which would have a negligible effect on the decrement. However, because of its variability and unusual spectrum, we conservatively assume a flat spectral index ($`\alpha =0`$) for $`\lambda <1`$ cm. Changing the 1 cm decrement to account for a faint central point source and altering the SuZIE 2.1 mm decrement to account for a 1.5 mJy point source, we find an increase of $`0.07\times 10^4`$ in $`y_\text{0}`$ and $`50\text{ km s}^1`$ in $`v_\text{p}`$. We include these in Table 3 as our systematic uncertainty due to point source contamination.

We test the importance of our assumed 12.4 keV electron temperature to $`y_\text{0}`$ and $`v_\text{p}`$ by performing the spectral fit over a large range in $`T_e`$. An increase in $`T_e`$ trades off with a decrease in the cluster optical depth $`\tau =n_e\sigma _T𝑑l`$, as can be seen in equation (7); thus so long as $`\tau `$ remains a free parameter, Compton $`y`$ is only weakly sensitive to temperature changes. In the $`1\sigma `$ range of $`T_e`$ ($`10.5\text{ keV}T_\text{e}15.2\text{ keV}`$), $`y_\text{0}`$ varies by $`{}_{0.07}{}^{+0.09}\times 10^4`$. The peculiar velocity also has a weak temperature dependence; it rises with temperature at low $`T_\text{e}`$ ($`5\text{ keV}`$) as $`\tau `$ decreases (cf. equation 8), then falls at high $`T_e`$ ($`15\text{ keV}`$) due to relativistic corrections. The $`1\sigma `$ range of $`T_e`$ corresponds to a broad plateau in the $`T_e`$-$`v_\text{p}`$ curve, and we find $`v_\text{p}`$ varying by only $`+20\text{ km s}^1`$ and $`40\text{ km s}^1`$.

We have assumed an isothermal temperature distribution for the cluster gas, consistent with the H1997a and b treatments of Abell 2163. Although the assumed cluster temperature has little impact on $`y_\text{0}`$ and $`v_\text{p}`$, uncertainties in the thermal structure of the cluster may have a larger effect. In a recent study of Abell 2163 using the BeppoSAX X-ray satellite, Irwin & Bregman (2000) find the cluster temperature profile consistent with isothermality at the $`2\sigma `$ confidence level out to a radius of $`9\mathrm{}`$; this is consistent with the result of White (2000), who finds a temperature profile consistent with isothermality at the $`1\sigma `$ confidence level out to $`10\mathrm{}`$. Chandra or XMM data are necessary to place stronger constraints on the thermal profile of Abell 2163, and we postpone quantifying a systematic uncertainty until such data become available.

Finally, we consider the effect of the assumed $`\pm 8\%`$ absolute calibration uncertainty in each SuZIE band and the $`\pm 3\%`$ absolute calibration uncertainty in our interferometric 1 cm measurements. These combine to change $`y_\text{0}`$ by as much as $`\pm 0.17\times 10^4`$ and $`v_\text{p}`$ by up to $`\pm 170\text{ km s}^1`$.

## 4 CONCLUSIONS

We have used interferometric data obtained at the OVRO and BIMA observatories to determine the SZE intensity decrement in the cluster Abell 2163 at a wavelength of 1 cm. We find $`\mathrm{\Delta }I_{\text{30}}`$$`=0.048\pm 0.006`$ MJy sr-1 by modeling the ICM as an isothermal spherical $`\beta `$-model. We have thus expanded the spectral coverage of this cluster to include the Rayleigh-Jeans part of the CMB blackbody spectrum and we present the most complete SZE spectrum to date. Upon fitting the OVRO/BIMA and SuZIE data for Abell 2163 to an SZE spectral model, we find $`y_\text{0}=3.56_{0.41}^{+0.41}`$$`{}_{0.19}{}^{+0.27}\times 10^4`$ and $`v_\text{p}=410_{850}^{+1030}`$$`{}_{440}{}^{+460}\text{ km s}_{}^{1}`$ where we list statistical uncertainty followed by systematic uncertainty at 68% confidence. The systematic uncertainties in Table 3 have been combined in quadrature to determine the totals reported above. We have shown that adding the 1 cm measurement improves constraints on the Compton $`y`$-parameter, as expected, but has less effect on the peculiar velocity which is best constrained by measurements near the 218 GHz thermal null.

The spectral fit includes corrections for dust emission in the SuZIE bands; we find that the contamination level depends strongly on the observing scheme of the instrument due to the distribution of dust on the sky. The dust contamination level for a cluster will obviously depend on its location on the sky; the line of sight to Abell 2163 passes near the Galactic plane, so it is more susceptible to dust contamination than clusters at higher galactic latitudes.

Finally, we emphasize that the Abell 2163 spectrum shown in Figure 2 consists of data from two very different experiments; an array of six bolometers performing drift scans (SuZIE) and two interferometers (OVRO and BIMA). The remarkable consistency that we find amongst the two sets of measurements is further evidence that SZE observations are now a reliable probe of structure in the universe.

This work is supported by NASA LTSA grant NAG5–7986. Many thanks to Erik Leitch for the use of his CMB field-generating software and for many helpful discussions. Thanks also to Amber Miller for useful discussions and much help with point sources. Clem Pryke provided invaluable information on the flux distribution of point sources detected by the DASI experiment. We thank the staff of the OVRO and BIMA observatories for their support and hard work, in particular J. R. Forster, J. Lugten, S. Padin, R. Plambeck, S. Scott, and D. Woody. Radio astronomy at the BIMA millimeter array is supported by NSF grant AST 96–13998, and the OVRO millimeter array is supported by NSF grant AST 96–13717. E. R. acknowledges support from NASA GSRP Fellowship NGT5–50173, and J. C. acknowledges support from a NSF–YI grant and the David and Lucile Packard Foundation. The Digitized Sky Survey was produced at the Space Telescope Science Institute under U.S. Government grant NAG W-2166. The images of these surveys are based on photographic data obtained using the Oschin Schmidt Telescope on Palomar Mountain and the UK Schmidt Telescope. The plates were processed into the present compressed digital form with the permission of these institutions. We thank S. Hansen, S. Pastor, and D. Semikoz for politely pointing out an error in our initial spectral fit. Finally, we thank M. Birkinshaw and M. Giard for the useful information they have provided.


text:
# Sudden Collapse of a Granular Cluster

## Abstract

Single clusters in a vibro-fluidized granular gas in $`N`$ connected compartments become unstable at strong shaking. They are experimentally shown to collapse very abruptly. The observed cluster lifetime (as a function of the driving intensity) is analytically calculated within a flux model, making use of the self-similarity of the process. After collapse, the cluster diffuses out into the uniform distribution in a self-similar way, with an anomalous diffusion exponent $`1/3`$.

One of the key features of a granular gas, making it fundamentally different from ordinary molecular gases, is its tendency to spontaneously separate into dense and dilute regions first6 . This clustering originates from the dissipative nature of the particle collisions. It is an unwanted effect in many applications where granular material is brought into motion. Therefore we study (within a simple geometry) how declustering occurs. We find that the breakdown of a cluster can be very abrupt, making declustering very different from clustering in reverse time order.

The experimental system (see Fig. 1) consists of a row of $`N`$ equal compartments, separated by walls of height $`h=25.0`$ mm and filled with a few hundred steel beads of diameter 3.0 mm. We start out with all the particles in the middle compartment and bring them into a gaseous state by shaking the system vertically. For weak shaking the cluster is stable: after some initial spilling, a dynamical equilibrium is established between the outflux of slow particles from the cluster and the influx of fast particles from outside eggers99 ; vdweele01 ; vdmeer01 . For sufficiently strong shaking, however, the cluster breaks down.

Two different regimes are observed: (i) At very strong shaking the breakdown occurs immediately, and the cluster spreads out over the boxes with its profile widening as $`t^{1/3}`$ (instead of the standard $`t^{1/2}`$ diffusion law). (ii) At less violent shaking, the cluster seems to remain stable for a long time, showing only a small leakage to its neighbors. But suddenly it collapses and subsequently diffuses over all boxes. The sudden death of the cluster is depicted in Fig. 1.

One thing this figure shows is that the breakdown of a cluster is strikingly different from the reverse process of cluster formation, which is known to take place gradually and (for all $`N3`$) via transient states showing clusters in several boxes vdweele01 ; vdmeer01 . This lack of time-reversibility is yet another consequence of the dissipation in the system.

The abruptness of the collapse allows us to define a cluster lifetime $`\tau `$ (via $`\ddot{n}_{cl}(\tau )=0`$, see Fig. 1e). In Fig. 2 the measured lifetimes are plotted as a function of the inverse shaking strength $`\stackrel{~}{B}`$ (Eq. 2) for various values of $`N`$. The data lie on a universal envelope curve, until at some critical value $`\stackrel{~}{B}_{c,N}`$ (which grows with $`N`$) they diverge.

All the above experimental observations can be accounted for quantitatively by means of the flux model of refs. eggers99 ; vdweele01 ; vdmeer01 . At the heart of this model is a flux function $`F(n_k)`$, describing the outflow from the $`k`$th box to each of its neighbors. It is a non-monotonic function of $`n_k`$ (the particle fraction in the box): $`F(n_k)`$ first increases with $`n_k`$, but beyond a certain value of $`n_k`$ it decreases again, as the increasingly frequent inelastic collisions slow the particles down, so that they cannot make it over the wall to the neighboring compartments anymore. The precise form of $`F(n_k)`$ is not very critical, as long as it is a one-humped function. We will use Eggers’ flux function eggers99

$$F(n_k)=C\sqrt{\stackrel{~}{B}}n_k^2e^{\stackrel{~}{B}n_k^2},$$
(1)
with

$$\stackrel{~}{B}\frac{ghr^4(1e)^2P^2}{S^2(af)^2}.$$
(2)
The driving parameter $`\stackrel{~}{B}`$ depends on the total number of particles $`P`$ and their properties (radius $`r`$, restitution coefficient $`e`$ of the particle collisions), on the geometry of the system (height $`h`$ of the walls, ground area $`S`$ of each box), and on the frequency $`f`$ and amplitude $`a`$ of the shaking. The factor $`C`$ only determines the absolute rate of the flux, and can be incorporated in the time scale.

The equation of motion for the fraction in box $`k`$ is art3note1

$$\frac{dn_k}{dt}=F(n_{k1})2F(n_k)+F(n_{k+1}),$$
(3)
where $`k=1,2,..,N`$. Here we assume a nearest neighbor interaction, and a cyclic arrangement of the boxes ($`k=N+1`$ equals $`k=1`$). We further impose particle conservation, $`_kn_k=1`$.

The numerical results shown in Fig. 1e and Fig. 2 have been obtained using the above flux model, starting out with all particles in one box (labeled $`cl`$). They quantitatively agree with the experimental observations. The decaying cluster goes through three different stages:

The starting stage is a very short one, in which both $`n_{cl}`$ and $`F(n_{cl})`$ display a jump compared to $`n_i`$ and $`F(n_i)`$ in the surrounding boxes, $`i=1,2,..`$ (we have to consider one side only because of the symmetry in the system).

In the second stage, the flux has become continuous but the particle fraction remains discontinuous. However, its low-density counterpart $`n_0`$ (defined by $`F(n_0)=F(n_{cl})`$) does continuously connect to $`n_1`$. We will use this fact later in the analysis of the envelope curve. The flux gradually grows, and eventually $`F(n_{cl})`$ reaches its maximum value. This is accompanied by rapid density changes and the sudden death of the cluster at the lifetime $`\tau `$.

In the third and last stage, both $`n_{cl}`$ and $`F(n_{cl})`$ fit continuously to the other boxes (see Fig. 3a). The remains of the cluster diffuse over the whole system until the uniform distribution is reached.

In what follows we will analytically solve the flux model. First we focus on the third stage. We rewrite the problem into its continuum version, by setting $`n(x,t)n_k(t)`$ ($`x=kw`$ by definition, where the box width $`w`$ will be incorporated in the $`x`$-scale). Eq. 3 then becomes:

$`_tn`$ $`=_{xx}F\left(n(x,t)\right)`$ (4)
$`=C\sqrt{\stackrel{~}{B}}_{xx}\left(n(x,t)^2e^{\stackrel{~}{B}n(x,t)^2}\right),`$
and the conservation condition takes the form $`_{\mathrm{}}^{\mathrm{}}n(x,t)𝑑x=1`$.

For very strong shaking (regime (i), where $`\tau `$ is vanishingly small) the diffusive stage sets in almost immediately. Here $`\stackrel{~}{B}0`$, and Eq. 4 reduces to

$$_tn=C\sqrt{\stackrel{~}{B}}_{xx}(n^2)=2C\sqrt{\stackrel{~}{B}}((_xn)^2+n_{xx}n),$$
(5)
which is known as the porous media equation barenblatt96 ; barwit . The decay of the cluster in this limit is depicted in Fig. 3a. It is self-similar: all curves in Fig. 3a fall onto a single curve if we properly rescale the axes (Fig. 3b). The original partial differential equation (PDE) can thus be brought back to an ordinary differential equation (ODE) in terms of the self-similarity variable $`\eta =x/(C\stackrel{~}{B}^{1/2}t)^{1/3}`$. With $`n(x,t)=H(\eta )/(C\stackrel{~}{B}^{1/2}t)^{1/3}`$, Eq. 5 now takes the form:

$`_{\eta \eta }(H^2)+{\displaystyle \frac{1}{3}}_\eta (\eta H)=0.`$ (6)
Its symmetric solution is $`H(\eta )=H_0(1/12)\eta ^2`$ (with the constant $`H_0=(3^{1/3})/40.361`$ determined by $`_{\mathrm{}}^{\mathrm{}}H(\eta )𝑑\eta =1`$). This inverted parabola, depicted in Fig. 3b, represents in one curve all the stages of Fig. 3a. The scaling of the axes shows that the height of the cluster decreases as $`t^{1/3}`$, and its width grows as $`t^{1/3}`$. This anomalous diffusion (with exponent $`1/3`$) is also found in porous media barenblatt96 ; barwit . The slowed down diffusion of the front originates from the quadratic $`n`$-dependence in $`_{xx}n^2`$ (Eq. 5). In the context of our granular model this reflects that particles only diffuse to neighboring boxes through 2-particle collisions: In the strong shaking limit the absence of particles slows down further diffusion, the presence enhances it.

For non-zero $`\stackrel{~}{B}`$ (regime (ii), where the diffusive stage has to wait until after the sudden death) we have an additional dimensionless variable $`\chi =\stackrel{~}{B}/C(t\tau )`$. Its influence diminishes with time and the solutions converge to the inverted parabola of the case $`\stackrel{~}{B}0`$.

Next we turn to the second, semi-continuous stage. For moderate shaking, $`\stackrel{~}{B}\stackrel{~}{B}_{c,N}`$, this stage can take quite a long time. At the critical point $`\stackrel{~}{B}_{c,N}`$ the lifetime $`\tau `$ even diverges to infinity and the cluster becomes stable. Starting out from the initial ..,0,0,1,0,0,.. distribution, the system first very slowly approaches a distribution in which $`n_{cl}`$ is close to the cluster density at the saddle-node bifurcation (see Fig. 2c) and all other boxes contain equal fractions $`n_k=(1n_{cl})/(N1)`$. Only when it has passed this phantom equilibrium (i.e., when $`n_{cl}`$ is below the level of the saddle-node bifurcation), the system quickens its pace and the sudden collapse occurs.

This means that $`\tau `$ is the time it takes to pass the phantom equilibrium. It can be calculated either numerically or analytically (by integrating the Taylor expansion of Eq. 3), with the result $`\tau (\stackrel{~}{B}_{c,N}\stackrel{~}{B})^{1/2}`$. So $`\tau `$ diverges as the inverse square root of the distance to the critical point, which is the common (mean field) power-law behavior near a second order phase transition as we have here for $`\tau `$ landau80 .

Finally, we calculate the lifetime away from the critical point, i.e., for $`\stackrel{~}{B}`$ considerably smaller than $`\stackrel{~}{B}_{c,N}`$. This will give us an analytical expression for the envelope curve in Fig. 2a. For these $`\stackrel{~}{B}`$ values, the collapse occurs before the particles leaking out of the cluster have had time to fill the outermost boxes to any significant level. Therefore, the behavior does not depend on the value of $`N`$: The system does not feel its finite size during the cluster’s lifetime, so the number of boxes can be taken to be infinite.

The time-evolution of the cluster is described by Eq. 3, with $`F(n_1)=F(n_1)`$:

$$\frac{dn_{cl}}{dt}=2F(n_{cl})+2F(n_1)=2F(n_0)+2F(n_1).$$
(7)
This equation contains $`n_1`$, which is governed by a similar equation of motion (Eq. 3) containing $`n_2`$, etc. So we have to deal with an infinite number of coupled nonlinear ODE’s art3note2 . This is a problem that cannot be solved directly, so we attack it in five steps.

Step 1: We first rewrite the problem into its continuum version, and replace the cluster density $`n_{cl}`$ by its low-density counterpart $`n(0,t)n_0(t)n_{cl}`$exp$`(\stackrel{~}{B}n_{cl}^2/2)`$. Thus, without influencing the fluxes (since $`F(n_0)=F(n_{cl})`$) we make $`n(x,t)`$ continuous in $`x=0`$. The density $`n(x,t)`$ obeys Eq. 4, plus a conservation condition saying that the increase of material into the rest of the system equals the influx from $`x=0`$:

$`_t{\displaystyle _0^{\mathrm{}}}n(x,t)𝑑x`$ $`=(_x\left(F\left(n(x,t)\right)\right)_{x=0}.`$ (8)
Step 2: The cluster-equation Eq. 7 now becomes:

$`{\displaystyle \frac{dn_{cl}}{dt}}`$ $`=2\left(_xF\left(n(x,t)\right)\right)_{x=0}`$ (9)
$`=2F^{}\left(n(0,t)\right)\left(_xn(x,t)\right)_{x=0}.`$
Since $`F^{}`$ can be derived directly from Eq. 1, the problem reduces to determining $`_xn(x,t)`$ at $`x=0`$.

Step 3: In order to do so, we observe that changes in $`n(0,t)`$ happen on a much longer timescale than in the surrounding boxes art3note3 , so the cluster acts as a constant reservoir spilling granular material. This approximation is illustrated in Fig. 4a: the profile in the system builds up while $`n_0`$ remains constant. In fact, this build-up takes place in a self-similar way (see Fig. 4b). So, as before, the problem for $`n(x,t)`$ can be formulated in terms of one variable $`\xi =x/(n_0C\stackrel{~}{B}^{1/2}t)^{1/2}`$. Setting $`n(x,t)=n_0G(\xi )`$, Eq. 4 becomes an ODE for $`G(\xi )`$, and also the accompanying conservation condition (Eq. 8) contains $`\xi `$ only:

$`{\displaystyle \frac{1}{2}}\xi _\xi G`$ $`=_{\xi \xi }(G^2e^{\beta G^2}),`$ (10)
$`\left(_\xi G\right)_{\xi =0}`$ $`={\displaystyle \frac{e^\beta }{4(1\beta )}}{\displaystyle _0^{\mathrm{}}}G(\xi )𝑑\xi ,`$
where $`\beta =\stackrel{~}{B}n_0^2`$.

Step 4: The slope $`(_xn(x,t))_{x=0}`$ can now be approximated by $`n_0(t)(_\xi G(\xi ))_{\xi =0}(_x\xi )_{x=0}`$, where we have revived the (slow) time-dependence in $`n_0(t)`$. With $`_x\xi =(n_0(t)C\stackrel{~}{B}^{1/2}t)^{1/2}`$, Eq. 9 becomes:

$$\frac{dn_{cl}}{dt}=2\frac{\sqrt{n_0(t)}}{\sqrt{C\stackrel{~}{B}^{1/2}t}}F^{}\left(n_0(t)\right)\left(_\xi G(\xi )\right)_{\xi =0}.$$
(11)
All quantities on the right hand side are tractable. We re-express $`n_0(t)`$ in terms of $`n_{cl}(t)`$, derive $`F^{}`$ from Eq. 1, and solve Eq. 10 to determine $`\left(_\xi G(\xi )\right)_{\xi =0}`$. This last step still requires some work, because Eq. 10 does not allow an analytical solution and moreover contains $`\beta `$ (and hence $`n_0(t)`$) explicitly. Since $`\beta `$ is small, however, we may expand $`G(\xi )`$ and Eq. 10 in terms of $`\beta `$, and solve numerically. In leading order we find $`\left(_\xi G(\xi )\right)_{\xi =0}=K=0.3138`$. Inserting all this, Eq. 11 becomes:

$$\frac{dn_{cl}}{dt}=4\frac{Kn_{cl}^{3/2}e^{(3/4)\stackrel{~}{B}n_{cl}^2}}{\sqrt{C\stackrel{~}{B}^{1/2}t}}\left(1\stackrel{~}{B}n_{cl}^2e^{\stackrel{~}{B}n_{cl}^2}\right).$$
(12)
This is an ODE for $`n_{cl}`$ in closed form, which replaces the original problem (Eq. 7) consisting of an infinite number of coupled ODE’s.

Step 5: Finally, we integrate Eq. 12 over the cluster density (dropping the suffix $`cl`$) and find an analytical expression for the lifetime $`\tau `$ away from the critical points:

$`\tau =\left[{\displaystyle _{n_{thr}}^1}{\displaystyle \frac{\sqrt{C\stackrel{~}{B}^{1/2}}n^{3/2}e^{(3/4)\stackrel{~}{B}n^2}dn}{8K\left(1\stackrel{~}{B}n^2e^{\stackrel{~}{B}n^2}\right)}}\right]^2.`$ (13)
Here $`n_{thr}`$ is the value of $`n`$ at which the sudden death occurs. For the evaluation of the solid curve in Fig. 2a we used $`n_{thr}=0.5`$, but this value is not too critical (cf. Fig. 1e). The only free parameter is the constant $`C`$: if this is adjusted properly, the analytical $`\tau `$ curve agrees with the measured data over the whole range of $`\stackrel{~}{B}`$-values.

The above expression shows that $`\tau `$ roughly increases as exp$`(\stackrel{~}{B}^{3/2})`$. Recalling that $`\stackrel{~}{B}`$ is the inverse shaking strength, this underlines the experimental observation that even a small reduction in the shaking strength causes a tremendous increase of the cluster lifetime.

In conclusion, in the studied compartmentalized system clusters break down very abruptly, in contrast to their slow formation. As clustering itself, the lack of time reversibility originates from the dissipative nature of the particle collisions: The breakdown of the unstable cluster is delayed because most of the energy input is dissipated through collisions in the cluster. The dynamics is quantitatively described by a remarkably simple flux model, which can be analytically solved.

###### Acknowledgements.

Acknowledgments: We would like to thank Ernst van Nierop for measuring the experimental lifetimes. This work is part of the research program of the Stichting FOM, which is financially supported by NWO.


text:
# Susceptibility of a Magnetic Impurity in Two-Dimensional Disordered Electron Systems

## Abstract

We investigate the Kondo effect in two-dimensional disordered electron systems using a finite-temperature quantum Monte Carlo method. Depending on the position of a magnetic impurity, the local moment is screened or unscreened by the spin of the conduction electron. The results suggest that the Kondo temperature takes different values depending on the position of a magnetic impurity. We show that the distribution of the Kondo temperature becomes wide and the weight at $`T_K=0`$ becomes large as randomness increases. The average susceptibility shows a weak power-law divergence at low temperature in strongly disordered systems, indicating a non-Fermi-liquid behavior. We calculate the correlation function between the local moment and the spin of the conduction electron. The results are discussed in connection with Kondo screening.

In low-dimensional disordered systems, the conduction electron localizes within a characteristic length called the localization length And ; AALR . It was shown that the localization effect enhances the electron correlation effects significantly, leading to conspicuous features particularly in low-dimensional systems AA ; HF ; LR ; Fink ; DiC ; BK . The Kondo effect is a typical phenomenon caused by the electron correlation around a magnetic impurity. The Kondo effect in disordered systems was investigated using the perturbation expansion in the weakly localized regime Ohkawa ; VZ ; Suga1 ; Suga2 . On the basis of the $`sd`$ model with random potentials, Ohkawa et al. Ohkawa  showed that the Kondo logarithmic terms for conductivity and susceptibility are modified into anomalous terms. Within the most divergent approximation, it was shown that higher order corrections consist of the product of such anomalous terms and the Kondo logarithmic terms, and that the latter is scaled into the same Kondo temperature as that in an ordinary Kondo system Suga1 . The temperature dependence of the conductivity and the susceptibility around and below the Kondo temperature in the weakly localized regime was calculated using the $`1/N`$ expansion based on the Coqblin-Schrieffer model Suga1  and the perturbation expansion based on the Anderson model Suga2 , respectively.

The Kondo effect in strongly disordered systems was studied by Dobrosavljević et al. by taking account of the Coulomb interaction among conduction electrons DKK . On the basis of a slave-boson approach, they showed that, as a result of the localization and interaction effects, the Kondo temperature is modified, and that the magnetic susceptibility and the coefficient of the linear specific heat show divergences as temperature approaches zero. In fact, a strong spatial distribution in the Kondo temperature was suggested from experimental results of strong broadening of the Cu NMR line of $`\mathrm{UCu}_{5\mathrm{x}}\mathrm{Pd}_\mathrm{x}`$ NMR . In spite of these findings, the effects of strong randomness itself on the behavior of a magnetic impurity have not yet been fully investigated from a microscopic viewpoint.

In this letter, we study the Kondo effect in two-dimensional (2D) strongly disordered electron systems. Using a finite-temperature quantum Monte Carlo (QMC) method Hirsch ; GHS , we calculate the susceptibility of a magnetic impurity and the correlation function between the local moment and the spin of the conduction electron.

Let us consider the single-impurity Anderson model with on-site random potentials described by the Hamiltonian

$`H`$ $`=`$ $`{\displaystyle \underset{i\sigma }{}}ϵ_ic_{i\sigma }^{}c_{i\sigma }t{\displaystyle \underset{<ij>\sigma }{}}c_{i\sigma }^{}c_{j\sigma }+ϵ_d{\displaystyle \underset{\sigma }{}}n_{d\sigma }`$ (1)
$`+V{\displaystyle \underset{\sigma }{}}(d_\sigma ^{}c_{0\sigma }+H.c.)+Un_dn_d`$
where random on-site potentials $`ϵ_i`$ are chosen to be a flat distribution in the interval $`[W,W]`$ under the condition $`_iϵ_i=0`$, $`<i,j>`$ denotes the summation of the nearest-neighbor sites, and $`n_{d\sigma }=d_\sigma ^{}d_\sigma `$. The system is a $`41\times 41`$ square lattice and a magnetic impurity is located at $`i=0`$. For $`W3.0`$, the conduction electron is localized with the localization length $`\xi (W)37.5`$ MacKinnon . Thus, the system for $`W3.0`$ is probably in the strongly localized regime at low temperature. We set the condition $`ϵ_d+(1/2)U=0`$. In the case without a random potential, this condition represents the electron-hole symmetry. We use the parameters $`U=2.0`$ and $`V=1.0`$ in units of $`t`$. The resonance width can be estimated as $`\mathrm{\Delta }\pi \rho _0V^20.39`$, where $`\rho _0`$ is the density of states at the Fermi energy and can be of the order of the inverse of the bandwidth ($`1/8`$). Therefore, the system is in the Kondo regime. The Kondo temperature without randomness can be estimated as TW  $`T_K^0=[V^2U/\pi ^2]^{1/2}\mathrm{exp}[(\pi /4)(U/V^2V^2/U)]0.14`$.

Using the QMC procedure reported in refs. 13 and 14, we calculate the susceptibility of a magnetic impurity,

$`\chi (T)=`$ (2)
$`\mathrm{\Delta }\tau {\displaystyle \underset{l=0}{\overset{L1}{}}}\left[n_d(\tau _l)n_d(\tau _l)\right]\left[n_d(\tau _0)n_d(\tau _0)\right],`$
where $`\mathrm{\Delta }\tau =\beta /L`$ with $`\beta =1/(k_\mathrm{B}T)`$ and $`L`$ is the Trotter number. We set $`g\mu _B/2=1`$. We first numerically diagonalize the Hamiltonian without local correlation effects between $`d`$ electrons to obtain the initial Green functions. We then calculate the local correlation effects by the QMC procedure, and obtain the Green functions needed to calculate physical quantities. We use the Trotter time-slice size $`\mathrm{\Delta }\tau =0.25`$. An open boundary condition is employed. In order to see the effect of the boundary condition on the low-energy properties, we have carried out QMC calculation, applying also the periodic boundary condition. The difference between the two results is negligibly small. We have run 50000 Monte Carlo sweeps for measurements after 1000 sweeps in the warming-up run, and have taken an average over the Monte Carlo sweeps. Note that the effect of a negative sign problem is negligibly small because of the conditions $`ϵ_d+(1/2)U=0`$ and $`_iϵ_i=0`$.

In Fig. 1(a), $`\chi (T)`$ versus $`\mathrm{log}T`$ is shown as a function of $`W`$. A magnetic impurity is located at the center of the $`41\times 41`$ square lattice. In the case without randomness, we observe an ordinary Kondo behavior. In this realization of a random potential, the susceptibility of a magnetic impurity shows a divergence for $`1.0W3.5`$ at low temperatures, while it shows a local Fermi-liquid behavior for $`W=0.5`$. Note that we have confirmed that the data in the latter case can be well fitted by $`\chi (T)=aT^2+b`$. We evaluate the divergent exponent $`\alpha `$ for $`\chi (T)T^\alpha `$ using data below $`T_K^00.14`$ by the least-squares method. The results are summarized in Table I. For $`W=1.0`$ and $`1.5`$, the susceptibility of a magnetic impurity diverges with a power weaker than that of a free spin, while for $`W2.0`$, a magnetic impurity behaves like a free spin. The results indicate that the impurity spin at this site is not completely screened by the spin of the conduction electron for $`W1.0`$.

Shifting the position of a magnetic impurity twenty-four times around the center in the same realization of the random potential, we calculate the susceptibility of a magnetic impurity. At seven positions, the susceptibility shows a local Fermi-liquid behavior irrespective of $`W`$, while at the remaining seventeen positions, the susceptibility shows a divergence behavior with increasing $`W`$. In the former case, a typical example is shown in Fig. 1(b). In contrast to the results in Fig. 1(a), the susceptibility at this site shows a local Fermi-liquid behavior for $`0W3.5`$ with monotonic decrease in $`\chi (0)`$ as a function of $`W`$. In the latter case, a logarithmic divergence can be seen at three positions, as well as a weak power-law divergence at nine positions. The results for weak power-law and logarithmic divergences are depicted in Fig. 2 for $`W=3.5`$. The divergent exponents are evaluated in the same way as in Fig. 1(a) and are presented in Table II. In Table III, the behavior of $`\chi `$ and $`ϵ_i`$ for $`W=3.5`$

at each position of a magnetic impurity is summarized. From the results, we conclude that the local moment in 2D disordered electron systems is screened or unscreened by the spin of the conduction electron, depending on the position of a magnetic impurity. Since $`\chi (0)=(2\pi T_K)^1`$ with $`T_K`$ being the Kondo temperature TW , this conclusion indicates that the Kondo temperature has a spatial distribution down to $`T_K=0`$.

On the basis of the results for the susceptibility of a magnetic impurity, we estimate the distribution of the Kondo temperature. When the susceptibility shows a Fermi-liquid behavior, we extrapolate $`\chi (0)`$ by the least-squares method using $`\chi (T)=aT^2+b`$. The results are summarized in Fig. 3. At $`W=1.0`$, the Kondo temperature is distributed around $`T_K=0.9T_K^0`$, where the largest weight exists, and a small weight lies at $`T_K=0`$. At $`W=3.5`$, the distribution of $`T_K`$ becomes wider and the weight at $`T_K=0`$ increases considerably with reduced weights for $`T_K0`$. We have shown that the spatial distribution of the Kondo temperature can be caused only by the effects of a random potential.

In the disordered system, we take an ensemble average over the distribution of random potentials in order to extract the universal properties. Judging from our numerical results, it seems that the susceptibility of a magnetic impurity does not show the universal behavior at low temperature. Therefore, in 2D strongly disordered electron systems with dilute magnetic impurities, where each magnetic impurity acts as a single magnetic impurity, the observable susceptibility at low temperature may be obtained as the arithmetic mean over the susceptibility of each magnetic impurity. The average susceptibility thus obtained in the random potential for $`W=3.5`$ realized in Fig. 1 is shown in Fig. 4. The average susceptibility for $`0.05<T<T_K^0`$ shows a weak power-law divergence $`\chi _{\mathrm{av}}T^\alpha `$ with $`\alpha 0.744`$, indicating a non-Fermi-liquid behavior.

To obtain local information about Kondo screening, we calculate the correlation function between the local moment at the $`0`$ site and the spin of the conduction electron at the $`i`$ site,

$`<S_0^zs_i^z>={\displaystyle \frac{1}{4}}<[n_dn_d][n_in_i]>,`$ (3)
where $`n_{i\sigma }=c_{i\sigma }^{}c_{i\sigma }`$. In Figs. 5(a) and 5(b), we show the spin correlation functions along the vertical and horizontal axes of the square lattice, respectively. The realized random potential and the position of the local moment are the same as those used in Figs. 1(a) and 1(b). Parameters used were $`W=3.0`$ and $`\beta =20`$. Thus, the results shown in Fig. 1(a) (open circles) correspond to the unscreened case, while those shown in Fig. 1(b) (full triangles) correspond to the screened case. In the former case, the antiferromagnetic correlation at the magnetic impurity is suppressed due to randomness. The suppression of the antiferromagnetic correlation at the magnetic impurity is commonly seen in other locations of a magnetic impurity where the susceptibility shows a divergence. This suppression probably causes the incomplete screening. Although the antiferromagnetic correlation at the magnetic impurity is not suppressed in the screened case, the behavior of the spin correlation function around a magnetic impurity is different from that for $`W=0`$; the ferromagnetic correlation is reduced.

In summary, we have investigated the Kondo effect in 2D disordered electron systems using a finite-temperature QMC method. We have shown that depending on the position of a magnetic impurity, its local moment can be screened or unscreened by the spin of the conduction electron. In the latter case, the antiferromagnetic spin correlation between the local moment and the conduction electron is suppressed around a magnetic impurity. We have demonstrated that the effects of randomness induce a distribution of the Kondo temperature down to $`T_K=0`$, and that the average susceptibility diverges with a weak power law at low temperature in strongly disordered systems.

We would like to thank K. Naito and T. Saso for useful comments and valuable discussions. Numerical computation was partly carried out at the Supercomputer Center, the Institute for Solid State Physics, University of Tokyo. This work was partly supported by a Grant-in-Aid for Scientific Research from the Ministry of Education, Culture, Sports, Science and Technology, Japan.


text:
# Back Reaction And Local Cosmological Expansion Rate

## I Introduction

Because of the nonlinear nature of the Einstein equations, linear fluctuations about a homogeneous and isotropic FRW cosmology will back-react on the background on which they live. The back-reaction of short wavelength gravitational waves was studied a long time ago by Brill, Hartle and Isaacson . In this approach, the back-reaction of fluctuations on the spatially averaged metric can be described by an effective energy-momentum tensor which contains the spatially averaged terms in the Einstein tensor which are second order in the amplitude of the fluctuations. More recently, Tsamis and Woodard  have initiated a detailed study of the back-reaction of long wavelength gravitational waves in a de Sitter background. Scalar-type metric fluctuations (called “cosmological fluctuations” in what follows) are also expected to back-react on the background space-time. Since in the context of inflationary and post-inflationary cosmology the scalar metric fluctuations are believed to dominate over the effects of gravitational waves, the back-reaction of these cosmological perturbations are expected to give the dominating back-reaction effect.

In , the formalism of  was generalized to describe the back-reaction of cosmological perturbations on the spatially averaged metric. On this basis, it was argued  that gravitational back-reaction in scalar field-driven inflationary models, calculated up to quadratic order in perturbations and to leading order in the long wavelength expansion and in the slow roll approximation could decrease the expansion rate of the universe and potentially solve the cosmological constant problem  (see also  for similar results in the context of back-reaction studies of gravitational waves). The formalism of  is covariant under first order space-time diffeomorphisms ***The back-reaction of small-scale (i.e. smaller than the Hubble radius) cosmological perturbations has been considered (without taking into account the issues of gauge freedom) in . This problem has also been considered in the context of Newtonian cosmological perturbation theory in . More recently, Nambu  has initiated a program to compute back-reaction effects on the spatially averaged metric using the renormalization group method..

However, as emphasized by Unruh  We are also grateful to Andrei Linde and Alan Guth for detailed private discussions on these points., the approach of  is deficient in several respects. First of all, due to the nonlinear nature of the Einstein equations, calculating an “observable” from the spatially averaged metric will not in general give the same result as calculating the spatially averaged value of the observable. More importantly, the spatially averaged metric is not a local physical observable. Thus, to take into account the deficiencies of the previous work on gravitational back-reaction, we must identify a local physical variable which describes the expansion rate of the Universe, calculate the back-reaction of cosmological perturbations on this quantity, and then take the spatially averaged value of the result. It is important to fix the hypersurface of averaging by a clear physical prescription in order to remove the possibility of being misled by effects which are second order gauge artifacts.

In this paper, we propose an implementation of this approach. We focus on the variable which yields the general relativistic definition of the local expansion rate, calculate this quantity to second order in the amplitude of the cosmological fluctuations, in terms of a time variable defined by an unambiguous physical prescription. For simplicity, we assume that matter is described by a single field (either a hydrodynamical field or a single scalar field).

We study two examples, first a matter dominated Universe, and second the inflationary phase of a cosmology dominated by a single scalar field. In both cases we find that the leading infrared contributions to the back-reaction on the local expansion rate of the Universe vanish, in contrast to the findings of the initial work on gravitational back-reaction of cosmological fluctuations , and confirming the analysis of . We thus confirm the conclusions reached recently in  where a different variable related to the local expansion rate is proposed, and different techniques to evaluate this variable are used.

Note that when evaluated at a fixed background time, our leading infrared back-reaction terms give a non-vanishing contribution. This leads us to the conjecture that in more realistic models in which a second field is present to determine time (e.g. the microwave background), the leading infrared back-reaction terms will not vanish.

## II A local observable

For a general perfect fluid flow in a curved space-time we consider the velocity four-vector field $`u^\alpha `$ tangential to a family of world lines. In the context of cosmology, we can always define a preferred family of world lines representing the motion of a set of comoving observers. In the case of hydrodynamical matter, this is easy since the energy-momentum tensor is already defined in terms of a velocity four-vector field. Also in the case of scalar field matter, a corresponding four-vector field can be defined, although a bit more care is required to obtain a consistent definition. In both cases we have

$$u^\alpha u_\alpha =\mathrm{\hspace{0.17em}1},$$
(1)
where $`\alpha `$ runs over the space-time indices. The projection tensor onto tangential three-surfaces orthogonal to $`u^\alpha `$ is

$$h_{\alpha \beta }=g_{\alpha \beta }u_\alpha u_\beta .$$
(2)
The first covariant derivative of the four-velocity can be decomposed as (see e.g.  for details )

$$u_{\alpha ;\beta }=\omega _{\alpha \beta }+\sigma _{\alpha \beta }+\frac{1}{3}\mathrm{\Theta }h_{\alpha \beta }\dot{u}_\alpha u_\beta .$$
(3)
Here

$$\mathrm{\Theta }u_{;\alpha }^\alpha $$
(4)
is the local expansion rate of the tangential surfaces orthogonal to the fluid flow, $`\omega `$ is the vorticity tensor (with $`\omega _{\alpha \beta }u^\beta =0`$), and $`\sigma `$ is the shear tensor (satisfying $`\sigma _{\alpha \beta }u^\beta =0,\sigma _\alpha ^\alpha =0`$).

For a homogeneous Universe with scale factor $`a(t)`$ the Hubble expansion rate $`H`$ is

$$H\dot{a}/a=\frac{1}{3}\mathrm{\Theta }.$$
(5)
For a cosmological model with fluctuations, $`\mathrm{\Theta }`$ is local in space and time. We will use $`\mathrm{\Theta }`$ to define the local expansion rate $`\dot{S}/S`$ via the local analog of (5), namely via the equation

$$\dot{S}/S\frac{1}{3}\mathrm{\Theta }.$$
(6)
The quantity $`\dot{S}/S`$ is a much better measure of the locally measured expansion rate in a Universe with fluctuations than the Hubble expansion rate used in  determined via the spatially averaged metric including back-reaction. It is a mathematically simpler object than the variable recently introduced in  which involves the integral along the past light cone of the observation point. If we are interested in evaluating the expansion rate for a typical observer, we propose to take the spatial average of the local expansion rate defined via (6).

Now that we have defined the observable we are interested in, the procedure will be as follows. First, we must determine the velocity four-vector field $`u^\alpha `$ for the models we are interested in. Then, we use the Einstein equations to express $`u^\alpha `$ in terms of the metric perturbation. Taking the relative amplitude of the metric fluctuations as the expansion parameter, we then calculate $`\mathrm{\Theta }`$, our local measure of the Hubble expansion rate, to second order. After evaluating the result on a physically determined hypersurface we can then study the back-reaction of cosmological fluctuations on the locally measured Hubble expansion rate. In this paper we will focus on the leading infrared contributions to back-reaction, the terms found to dominate the back-reaction effects in .

## III Deriving the Expansion Rate for Scalar Metric Perturbations

In this section we consider a model with hydrodynamical matter. Starting from the expression for the metric to linear order in the fluctuations $`\mathrm{\Phi }`$ (see  for a detailed review), we determine the velocity four vector field $`u^\alpha `$ to the order required to analyze the leading infrared terms in the back-reaction to quadratic order. To obtain the full back-reaction terms (including terms which dominate in the ultraviolet but are negligible in the infrared) we should calculate $`u^\alpha `$ consistently up to second order. However, if we are only interested in the leading infrared terms, it is sufficient to keep all the terms quadratic in $`\mathrm{\Phi }`$ but not containing any spatial gradients.

In order to obtain the complete result for gravitational back-reaction we would have to look at the Einstein equations for a perfect fluid with energy density $`\rho `$ and pressure $`P`$,

$$G_{\mu \nu }=(P+\rho )u_\mu u_\nu Pg_{\mu \nu }$$
(7)
(in units in which $`8\pi G=1`$), which, since $`G_\mu ^\mu =R`$, will yield

$`R`$ $`=`$ $`\rho 3P`$ (8)
$`\rho `$ $`=`$ $`u^\mu G_{\mu \nu }u^\nu ,`$ (9)
and lead to an equation that can be solved perturbatively to any desired order for $`u_i`$:

$$G_i^0=\frac{4}{3}u^\mu G_{\mu \nu }u^\nu u^0u_i+Ru^0u_i$$
(10)
However, as mentioned above, here we just use the results for $`u^i`$ which are of linear order. Since we will calculate the divergence of $`u^\mu `$, our prescription implies that we are ignoring some of the extra second order gradient terms . The result can also be used for scalar fields if we define the $`u`$ vector field in a proper way.

For an unperturbed Robertson-Walker metric, the four-velocity field $`u`$ in comoving coordinates would be

$$u^\mu =(1,0,0,0)$$
(11)
In linear perturbation theory, and in the case of simple forms of matter (such as a single fluid or a single scalar field) for which there is to linear order no anisotropic stress, the metric (in longitudinal gauge) can be written as

$`ds^2`$ $`=`$ $`a(\eta )^2\left((1+2\mathrm{\Phi })d\eta ^2(12\mathrm{\Psi })\gamma _{ij}dx^idx^j\right),`$ (12)
$`\gamma _{ij}`$ $`=`$ $`\delta _{ij}[1+{\displaystyle \frac{1}{4}}𝒦(x^2+y^2+z^2)]`$ (13)
where $`𝒦=0,1,1`$ depending on whether the three-dimensional space corresponding to the hypersurface $`t=\mathrm{const}.`$ is flat, closed or open. In this paper we will take it to be zero in order to simplify the calculations. The time variable $`\eta `$ appearing in (12) is conformal time and is related to the physical time $`t`$ via $`d\eta =a^1dt`$. For the forms of matter considered here, $`\mathrm{\Psi }=\mathrm{\Phi }`$ at linear order Even if we did not make the assumption $`\mathrm{\Psi }=\mathrm{\Phi }`$, it turns out that at second order all infrared terms depending on $`\mathrm{\Psi }\mathrm{\Phi }`$ will drop out.. As discussed e.g. in , in longitudinal gauge the spatial components of the four-velocity vector field are related to $`\mathrm{\Phi }`$ via

$$\delta u_i=a^2(^2^{}+𝒦)^1(a\mathrm{\Phi })_{,i}^{}$$
(14)
where a prime denotes differentiation with respect to $`\eta `$ and $`=a^{}/a`$. Using equation (1) we can derive the expression for the time component of $`u^\alpha `$ in terms of $`\mathrm{\Phi }`$:

$$u^0(\eta )=a^1(1\mathrm{\Phi }+\frac{3}{2}\mathrm{\Phi }^2)$$
(15)
Now that we have all components of $`u^\alpha `$, we take the covariant derivative of it and retain all $`\mathrm{\Phi }`$ dependence up to second order §§§This is in the philosophy of the general back-reaction approach in which it is assumed that the fluctuations of the metric and matter satisfy the linear perturbation equations, and we compute their back-reaction on physical quantities to second order. It is not a consistent second order perturbative formalism. These second order terms principally come from the Christoffel symbols. Other second order terms (which as mentioned before are gradient terms) could be added to the ones computed here if we were to solve the Einstein equations in the form (10) beyond linear order. A straightforward calculation yields

$`\mathrm{\Theta }`$ $`=`$ $`\mathrm{\hspace{0.17em}3}{\displaystyle \frac{a^{}}{a^2}}(1\mathrm{\Phi }+{\displaystyle \frac{3}{2}}\mathrm{\Phi }^2)3{\displaystyle \frac{\mathrm{\Phi }^{}}{a}}`$ (16)
$`+`$ $`{\displaystyle \frac{(a_i\mathrm{\Phi })^{}(_i\mathrm{\Phi })+(a_i^2\mathrm{\Phi })^{}}{a^2(^{}^2)}}.`$ (17)
Since $`\dot{a}/a=a^{}/a^2`$, we can immediately read off the extra terms contributing to the local expansion rate which result from the presence of cosmological fluctuations. Upon spatial averaging at a fixed conformal time, the terms linear in $`\mathrm{\Phi }`$ drop out. Hence, it follows that if evaluated at a fixed conformal time, infrared modes on average lead to an increase in the expansion rate compared to what would be obtained at the same conformal time in the absence of metric fluctuations. Whether this is a physically measurable effect from an observational point of view will be discussed in more depth in following sections.

## IV Expansion Rate for a Matter Dominated Universe

Now let us use the result of the last section to (as an example) calculate the local Hubble expansion rate for a matter-dominated Universe. In this case, the scale factor $`a`$ and the scalar metric perturbation $`\mathrm{\Phi }`$ have the following dependence on the conformal time $`\eta `$ :

$`a(\eta )`$ $`=`$ $`a_m\eta ^2/2`$ (18)
$`\mathrm{\Phi }(\eta ,x)`$ $`=`$ $`C_1(x)+C_2(x)\eta ^5,`$ (19)
where $`a_m`$ is a constant, and $`C_1`$ and $`C_2`$ are time-independent. The second equation is valid in the long wavelength (super-Hubble-scale) limit and has been explicitly derived in . Now using equations (16), (18) and (19) we obtain $`\mathrm{\Theta }`$ in terms of conformal time:

$`\mathrm{\Theta }`$ $`=`$ $`{\displaystyle \frac{12}{a_m}}\eta ^3(1C_1(x)+{\displaystyle \frac{3}{2}}C_1^2)`$ (20)
$``$ $`{\displaystyle \frac{2}{3}}{\displaystyle \frac{1}{a_m}}\eta ^1((_iC_1(x))^2+_i^2C_1(x))`$ (21)
$`+`$ $`{\displaystyle \frac{1}{a_m}}\eta ^6({\displaystyle \frac{1}{3}}(_iC_1(x))(_iC_2(x))+_i^2C_2(x))`$ (22)
$`+`$ $`{\displaystyle \frac{1}{a_m}}\eta ^8(3C_2(x)+36C_1(x)C_2(x))`$ (23)
$`+`$ $`{\displaystyle \frac{1}{a_m}}\eta ^{11}(_iC_2(x))^2+{\displaystyle \frac{1}{a_m}}\eta ^{13}(C_2(x))^2`$ (24)
Some of the terms are decreasing very fast as a function of time and thus we can ignore them. If we just keep the terms with the powers $`3`$ and $`1`$ of $`\eta `$ then we get:

$`\mathrm{\Theta }`$ $`=`$ $`{\displaystyle \frac{12}{a_m}}\eta ^3(1C_1(x)+{\displaystyle \frac{3}{2}}C_1^2)`$ (25)
$``$ $`{\displaystyle \frac{2}{3}}{\displaystyle \frac{1}{a_m}}\eta ^1[(_iC_1(x))^2+_i^2C_1(x)]`$ (26)
If we take the average of $`\mathrm{\Theta }`$ on a constant $`\eta `$ hypersurface, only the terms quadratic in the fluctuation variables survive. Thus, considering large values of $`\eta `$ and focusing on the second order terms, it appears from (25) that infrared modes give a positive contribution to $`\mathrm{\Theta }`$ and thus lead to a speeding-up of the expansion, whereas ultraviolet modes enter with a negative sign and thus yield a slowing effect, the latter becoming more significant (relative to the unperturbed expansion rate) for larger values of $`\eta `$.

However, before drawing definite physical conclusions from our analysis, we must take into account that the background time $`\eta `$ is not an observable quantity. To obtain results for back-reaction of any real physical significance we have to find an observable variable like proper time and evaluate the expansion rate in terms of this variable, so that we can discuss its evolution from an observer’s point of view.

If we use equation (12) for the metric, equation (19) for $`\mathrm{\Phi }`$, and equation (18) for the scale factor, we can find the expression for the proper time $`\tau `$ in terms of conformal time. Since

$$d\tau ^2=a(\eta )^2(1+2\mathrm{\Phi })d\eta ^2,$$
(27)
a simple integration yields

$`\tau `$ $`=`$ $`{\displaystyle \frac{a_m}{6}}(1+C_1{\displaystyle \frac{1}{2}}C_1^2)\eta ^3`$ (28)
$``$ $`{\displaystyle \frac{a_m}{4}}(C_2C_1C_2)\eta ^2+{\displaystyle \frac{a_m}{28}}C_2^2\eta ^7`$ (29)
In the approximation of large values of $`\eta `$ we can ignore the second and the third term of this equation and thus obtain:

$$\eta ^3=\frac{a_m}{6}(1+C_1\frac{1}{2}C_1^2)\tau ^1.$$
(30)
Now we can use this relation and substitute into the first term of equation (25). We see the effects of the dominant infrared terms on the local expansion rate cancel exactly up to second order in perturbations, when evaluating $`\mathrm{\Theta }`$ at a fixed proper time:

$$\mathrm{\Theta }_{IR}=2\tau ^1.$$
(31)
This implies that at least in the approximation of keeping only the leading infrared terms, there is no local gravitational back-reaction of cosmological fluctuations on the local Hubble expansion rate in this matter-dominated universe.

## V Expansion Rate in Terms of Scalar Field as an Observable

We now move on to the example more relevant to the work of , namely a Universe dominated by a single real scalar field $`\phi `$,which is a toy model for inflationary cosmology. During inflation, fluctuations which are generated on sub-Hubble scales early on during the inflationary phase are red-shifted to scales much larger than the Hubble radius. Thus, in this context it is of great interest to consider the back-reaction of infrared modes.

In the following, we will generalize the previous analysis to be applicable to matter consisting of a single scalar field. In this case one can treat the scalar field as a perfect fluid and derive the velocity four-vector field. To do this, we need to write the energy-momentum tensor $`T_{\mu \nu }`$ of the scalar field (right hand side of the following equation) in the form of an energy-momentum tensor for a perfect fluid (left hand side of the following equation):

$$(\rho +P)u_\mu u_\nu Pg_{\mu \nu }=_\mu \phi _\nu \phi .g_{\mu \nu }$$
(32)
At the level of the background fields, the two expressions are identical if we take $`P=`$ and $`u_\mu =A_\mu \phi `$ with

$$A=\left(^\nu \phi _\nu \phi \right)^{1/2}.$$
(33)
Now that we have shown that the energy-momentum tensor of a scalar field can be written in the form of that of a perfect fluid, we can use the expression (16) which gives $`\mathrm{\Theta }`$ in terms of the metric fluctuations to evaluate the local effect of gravitational back-reaction of cosmological fluctuations. Let us first for convenience rewrite $`\mathrm{\Theta }`$ of (16) in terms of the physical time $`t`$:

$$\mathrm{\Theta }=\mathrm{\hspace{0.17em}3}\frac{\dot{a}}{a}(1\mathrm{\Phi }+\frac{3}{2}\mathrm{\Phi }^2)3\dot{\mathrm{\Phi }}+\frac{(_i\dot{\mathrm{\Phi }})(_i\mathrm{\Phi })+_i^2\dot{\mathrm{\Phi }}}{a\ddot{a}\dot{a}^2}$$
(34)
Theoretically, the scalar field $`\phi `$ is an observable. In fact, in a system with a single matter field $`\phi `$, it is this field which must be used as a clock. Hence, to obtain results with physical meaning, we must evaluate $`\mathrm{\Theta }`$ on a surface of constant $`\phi `$ and not constant $`t`$.

As discussed above, in the context of inflationary cosmology it is important to study the effects of infrared modes up to second order. We will also assume that we are in the slow rolling regime of inflation. In our case (as well as for a more general case), the prescription is to calculate $`t`$, $`\mathrm{\Phi }`$, $`a`$ and $`\frac{}{t}`$ in terms of $`\phi `$, and to insert the results into the general expression (34) for $`\mathrm{\Theta }`$. The relation between $`t`$ and $`\phi `$ can be derived starting from:

$$\phi (t)=\phi _0(t)+\delta \phi _1(t),$$
(35)
which can be written as:

$$t=\phi _0^1[\phi \delta \phi _1(t)].$$
(36)
Thus, the equation

$`t`$ $`=`$ $`\phi _0^1(\phi ){\displaystyle \frac{\phi _0^1(\phi )}{\phi }}\delta \phi _1(\phi _0^1(\phi ))`$ (37)
$`+`$ $`({\displaystyle \frac{\phi _0^1(\phi )}{\phi }})^2{\displaystyle \frac{\delta \phi _1(\phi _0^1(\phi ))}{t}}\delta \phi _1(\phi _0^1(\phi ))`$ (38)
$`+`$ $`{\displaystyle \frac{1}{2}}{\displaystyle \frac{^2\phi _0^1(\phi )}{\phi ^2}}\delta \phi _1^2(\phi _0^1(\phi ))`$ (39)
relates $`t`$ and $`\phi `$.

We now wish to express the value of $`\mathrm{\Theta }`$ in terms of $`\phi `$ (note that we are considering the local value of $`\mathrm{\Theta }`$ in this analysis, and there is no need to perform a spatial averaging). In order to relate the metric, its fluctuations and the scale factor to the scalar field and its fluctuations we need to take an explicit form for the potential and make use of the Einstein constraint equations. The $`G^{0i}`$ and $`G^{00}`$ equations relate $`\mathrm{\Phi }`$ to $`\phi `$ (at the level of the first order fluctuations) and the Hubble parameter of the background to the background scalar field $`\phi _0`$ (at the level of the unperturbed Friedmann equations). In our simple system we do not need to go through these calculations explicitly (as long as we are interested only in the leading infrared terms) since by dropping the gradient terms from the $`G^{0i}`$ and $`G^{00}`$ equations, it can be shown that

$$\frac{H}{\sqrt{1+2\mathrm{\Phi }}}=\frac{1}{\sqrt{3}}\sqrt{V(\phi )}.$$
(40)
If we expand the right hand side in terms of $`\mathrm{\Phi }`$ we get the result

$$H(1\mathrm{\Phi }+\frac{3}{2}\mathrm{\Phi }^2)=\frac{1}{\sqrt{3}}\sqrt{V(\phi )},$$
(41)
and substituting this result into equation (34) and neglecting terms containing $`\dot{\mathrm{\Phi }}`$ and $`_i\mathrm{\Phi }`$ (which are sub-dominant compared to other terms in the infrared and slow roll limit) leads to the final result for the local expansion rate:

$$\mathrm{\Theta }=\sqrt{3}\sqrt{V(\phi )}$$
(42)
which as a function of $`\phi `$ is the same as the relation for an unperturbed background.

Thus, again we do not see any back-reaction of cosmological perturbations on the local expansion rate in this approximation. In retrospect it is easy to understand this result, since  by neglecting gradient terms we have a Friedmann Universe whose expansion rate satisfies equation (42).

## VI Conclusions

In this paper, we have studied the back-reaction effects on a local observable which measures the local expansion rate of the Universe. The observable gives the rate at which neighboring comoving observers separate and coincides with the usual definition of expansion in the context of the fluid approach to cosmology. In order to obtain a physical quantity, we evaluated the observable at a fixed value of the scalar field.

We evaluated our observable, the local physical expansion rate, in a simple toy model of chaotic inflation consisting of a single scalar matter field coupled to gravity. We found that the leading infrared terms, the terms which dominate the effects discussed in  and , cancel if we evaluate the observable at a fixed value of the scalar field, whereas they do not vanish if we evaluate them at a fixed value of the background time. The former result is a physical result since it corresponds to a physical observable evaluated at a space-time point specified by a physical prescription, whereas the latter result (obtained by evaluating at a fixed background time) does not have a diffeomorphism-invariant meaning. Our analysis thus confirms the concern of  that the results obtained in  and  are not invariant under second order gauge transformations. Our results confirm the conclusions of  reached by means of a different method of analysis applied to a different physical observable.

Our result does not imply that there is no back-reaction of the infrared modes of cosmological perturbations. There is no reason to expect that the next to leading infrared terms in our result will cancel (they do not cancel in the analysis of ). One of the advantages of our technique is that they can be evaluated without too much trouble. This is left to a future publication. So, even in single field models of inflation there might be some non-vanishing back-reaction of infrared modes.

We expect that back-reaction of infrared modes will be much more important in two field models of inflation. Let us assume that the matter sector of the theory contains both an inflation field $`\phi `$ and a regular matter field $`\chi `$ (with nonvanishing and time-dependent spatial average) which, for example, could represent the cosmic microwave background. In this case, it is no longer true that long wavelength fluctuations have no physical effects on local observables. If the measurement point is (in an unambiguous physical way) determined by a fixed value of the field $`\chi `$, then the local expansion rate may sensitively depend on the amplitude of the long wavelength fluctuations in $`\phi `$. Thus, the leading infrared terms may not cancel when evaluated according to the abovementioned physical prescription in the same way that they do not cancel in the analysis of this paper when the observable is evaluated at a fixed value of the background time.

There is a close analogy with the analysis of the parametric amplification of super-Hubble-scale cosmological fluctuations during inflationary reheating. From the point of view of the background space-time coordinates, it appears  that the parametric amplification of matter fluctuations on super-Hubble scales in an unperturbed cosmological background (see e.g.  for a discussion of parametric resonance during reheating) would imply the parametric amplification of the cosmological fluctuations on these scales. However, it can be shown that in single field models physical observables measuring the amplitude of cosmological fluctuations do not feel any resonance . In contrast, in two field models of inflation there is  parametric amplification of super-Hubble-scale cosmological fluctuations. In this case, there is a fluctuation mode corresponding to entropy fluctuations which cannot locally be gauged away. This mode is (in certain theories) parametrically amplified during reheating, and in turn drives the parametric resonance of the super-Hubble scale curvature fluctuations.

Our techniques allow us to calculate the back-reaction of cosmological fluctuations in two field models in a very similar way to what is presented here. Results will be presented in a followup publication .

Acknowledgments

We would like to thank Bill Unruh, Alan Guth and Andrei Linde for discussions which stimulated this research. We are grateful to Raul Abramo, Niayesh Afshordi, Yasusada Nambu and Richard Woodard for many useful insights during the course of this research. This research is supported in part (at Brown) by the US Department of Energy under Contract DE-FG0291ER40688, Task A. One of us (GG) is grateful to the Department of Astrophysical Sciences of Princeton University for hospitality during the course of this work. The other of us (RB) wishes to thank the CERN Theory Division and the Institut d’Astrophysique de Paris for hospitality and financial support.