text:
# Angular Clustering with Photometric Redshifts in the Sloan Digital Sky Survey: Bimodality in the Clustering Properties of Galaxies

## 1. Introduction

One of the primary tools for studying the evolution and formation of structure within the universe has been the angular correlation function (Totsuji & Kihara, 1969). Possibly the simplest of these point process statistics is the angular 2-point function which measures the excess number of pairs of galaxies, as a function of separation, when compared to a random distribution. If the universe can be described by a Gaussian random process then the 2-point function will fully describe the clustering of galaxies. While this is clearly not the case and higher order correlation functions play a significant role in understanding the clustering of structure and its evolution, the 2-point function remains an important statistical tool. In this paper we will utilize the angular 2-point function to determine the type and luminosity dependence of the clustering of galaxies within the Sloan Digital Sky Survey (SDSS; York et al., 2000).

Studying the angular correlation function has a natural advantage over the spatial or redshift correlation function. By only requiring positional information we can derive the clustering signal from photometric surveys alone (i.e. without requiring spectroscopic followup observations). Given the relative efficiencies of photometric surveys over their spectroscopic counterparts, this enables the correlation function to be estimated for wide-angle surveys covering statistically representative volumes of the universe and without being limited by discreteness error. The disadvantage of the angular correlation function has been that it is the projection of the spatial correlation function over the redshift distribution of the galaxy sample. For bright magnitude limits the redshift distribution is well known (and relatively narrow) and therefore deprojecting the angular clustering to estimate the clustering length is relatively straightforward. At fainter magnitudes the redshift distribution becomes broader and the details of the clustering signal can be washed out.

We can overcome many of the disadvantages of the angular clustering if we utilize photometric redshifts. Photometric redshifts provide a statistical estimate of the redshift, luminosity and type of a galaxy based on its broadband colors. As we can control the redshift interval from which we select the galaxies (and the distribution of galaxy types and luminosities) we can determine how the clustering signal evolves with redshift and invert it accurately to estimate the real space clustering length, $`r_0`$, for galaxies. The ability to utilize large, multicolor photometric surveys as opposed to the smaller spectroscopic samples means that we can subdivide the galaxy distributions by luminosity and type without being limited by the size of the resulting subsample (i.e. most of our analyzes will not be limited by shot noise). As we expect the dependence of the clustering signal to vary smoothly with luminosity, type and redshift, it is not expected that the statistical uncertainties in the redshift estimates will significantly bias our resulting measures.

The utility of photometric redshifts for measuring the clustering of galaxies as a function of redshift and type has been recognized when studying high redshift galaxies (Connolly et al., 1998; Magliocchetti & Maddox, 1999; Arnouts et al., 1999; Roukema et al., 1999; Brunner, Szalay & Connolly, 2000; Teplitz et al., 2001; Firth et al., 2002). In this paper we focus on studying the angular clustering of intermediate redshift galaxies $`z<0.3`$ and the dependence of the clustering length on luminosity and galaxy type. This represents one of the first applications of the photometric redshifts to the clustering of intermediate redshift galaxies for which we have a large, homogeneously and statistically significant sample of galaxies. This paper is divided into five sections. In Section 2 we describe the data set used in this analysis and the selection of a volume limited sample of galaxies. In Section 3 we apply a novel approach for estimating the 2-point angular correlation function using Fast Fourier Transforms and we show the dependence of the slope and amplitude of the correlation function on luminosity and galaxy type. In Section 4 we invert the projected angular correlation function and derive the correlation lengths. In Section 5 we discuss the bimodal behavior of clustering properties.

## 2. Defining a Photometric Sample

The SDSS is a photometric and spectroscopic survey designed to map the distribution of stars and galaxies in the local and intermediate redshift universe (SDSS; York et al., 2000). On completion the SDSS will have observed the majority of the northern sky ($`\pi `$ steradians) and approximately 1000 square degrees in the southern hemisphere. These observations are undertaken in a drift-scan mode where a dedicated 2.5m telescope scans along great circles, imaging 2.5 degree wide stripes of the sky. The imaging data is taken using a mosaic camera (Gunn et al., 1998) through the five photometric passbands $`u^{}`$, $`g^{}`$, $`r^{}`$, $`i^{}`$ and $`z^{}`$ (covering the ultraviolet through to the near infrared) as defined in Fukugita et al. (1996). The photometric system is described in detail by Smith et al. (2002) and the photometric monitor by Hogg et al. (2001). All data are reduced by an automated software pipeline (Lupton et al., 2003; Pier et al., 2002) and the outputs loaded in to a commercial SQL database. In our analysis, we will include data from runs with the longest contiguous scans (typically in excess of 50 degrees). These data comprise a subset of the data that will be released to the public as part of data release (DR1; Abazajian et al., 2003). In comparison to the Early Data Release (Stoughton et al., 2002) the area analyzed in this paper is approximately $`10`$ times larger, or approximately 20% of the entire survey area.

Given the five band photometry from the SDSS imaging data we estimate the photometric redshifts of the galaxies using the techniques outlined in Budavári et al. (2000, 1999); Csabai et al. (2000); Connolly et al. (1999). The details of the estimation techniques employed together with the expected uncertainties within the redshift estimators are given in Csabai et al. (2002). In this paper we will just note the effective rms error of the photometric redshifts (typically $`\mathrm{\Delta }z_{\mathrm{rms}}=0.04`$ at $`r^{}<18`$). For all sources within a sample the redshift and its uncertainty are calculated together with a measure of the spectral type of the galaxy and its variance and covariance (with redshift). From these measures we estimate the luminosity distance to each galaxy and calculate its r-band absolute magnitude.

### 2.1. Building the Sample Database

From the current photometric data in the SDSS archives we extract eight stripes for our analysis. These stripes combine to form approximately five coherent regions on the sky which range from approximately 90 degrees to about 120 degrees in length. In the nomenclature of the SDSS these stripes are designated the numbers 10–12, 35–37 and 76 and 86. The last two stripes come from the southern component of the survey. All data from these stripes have been designated as having survey quality photometric observations and astrometry. In total these stripes add up close to 20 million galaxies and are accessible through the SDSS Science Archive (Thakar et al., 2001).

Currently, there are two versions of the SDSS science archive running at Fermilab and remotely accessible to the collaboration. The “chunk” database contains stripes that have passed through the target selection process for identifying candidates for spectroscopic followup but with only photometry that was available when the spectroscopic target selection was run on the region (i.e. photometry for which the calibrations were not necessarily optimal). The “staging” database has the latest photometric data (for this paper we use the ver. 5.2.8 of the photometric pipeline) but without the full target selection information. For our purposes the quality of the photometric measurements is important but the target selection is completely irrelevant and thus we take our sample from the staging database.

In order to be able to efficiently store, search and select galaxies from the catalog, we create a local database using Microsoft’s SQL Server. The relevant properties (position, redshift, galaxy type, absolute and apparent magnitudes) of all galaxies in the staging database were stored in this database server. The regions of the sky surveyed by these data, the seeing of the observations as a function of position on the sky and the position of bright stars that must be masked out when defining the survey geometry were all calculated internally from this data set.

From these data we restrict our sample to galaxies brighter than $`r^{}=21`$. At this magnitude limit the star-galaxy separation is sufficiently accurate that it will not affect the angular clustering measures (Scranton et al., 2003) and the photometric redshift errors are typically less than $`\sigma =0.06`$ (Csabai et al., 2002). Applying this magnitude limit yields approximately 13 million galaxies from which to estimate the clustering signal. The sample was further restricted by excluding those regions of the stripes affected by the wings of bright stars or that were observed with poor seeing. Figure 1 shows the density of galaxies in stripe 11 with the masks over-plotted. Fields with seeing worse than $`1.7^{\prime \prime }`$ and a $`3^{}\times 3^{}`$ neighborhood around all bright stars with $`r^{}<14`$ were discarded.

We note that these selections were all accomplished by applying a series of SQL queries to the database rather than progressively pruning a catalog of galaxies. The boundaries of the stripes, the seeing on a field-by-field basis and all bright stars were stored in the local database, so that masks could be generated on the fly over the area to be analyzed. As such the selection criteria that were applied to the database could be optimized in a relatively short period of time.

### 2.2. Clustering from a Volume Limited Sample

Often the clustering evolution of galaxies, particularly that defined by angular clustering studies, is characterized as a function of limiting magnitude. While observationally this is simple to determine, the results of these analyzes are often difficult to interpret because in a magnitude limited sample the mix of the spectral types and absolute luminosities of galaxies is redshift dependent. We are, essentially, looking at the clustering properties of different types of galaxies as a function of limiting magnitude. The models to account for the clustering signal must, therefore, also be able to describe the evolution of the distribution of galaxies. Ideally we would separate out the effects of population mixing and study the evolution of angular clustering in terms of the intrinsic properties of galaxies (i.e. rest-frame color and luminosity), along with their distances. We can accomplish this if we use the photometric redshifts to select a volume limited sample of galaxies (i.e. one with a fixed absolute magnitude range as a function of redshift).

The SDSS Early Data Release photometric redshift catalog by Csabai et al. (2002) is based on techniques (Budavári et al., 2000, 1999; Csabai et al., 2000; Connolly et al., 1999) that estimate the physical parameters describing the galaxy samples in a self-consistent way. The relationship for the EDR data set is applied to the photometric data selected from the imaging stripes and the spectral types, absolute magnitudes and k-corrections are stored within a database together with the photometric and positional information. All derived quantities assume an $`\mathrm{\Lambda }`$ cosmology with $`h=0.7`$, $`\mathrm{\Omega }_\mathrm{M}=0.3`$ and $`\mathrm{\Omega }_\mathrm{\Lambda }=0.7`$.

Figure 2 illustrates how the absolute magnitude limits in the $`r^{}`$ band vary as a function of redshift. These absolute magnitude boundaries are well defined, as a function of redshift, by the galaxies within our sample. In this paper we study a volume-limited sample that extends out to a redshift $`z=0.3`$ with limiting absolute magnitude $`M_r^{}=19.97`$. We further restrict the data set to those galaxies more distant than redshift $`z=0.1`$. The reason for this lower redshift limit is that the main spectroscopic galaxy sample will contain many of the low redshift galaxies, so using photometric redshift estimates is not necessary. Also the fractional error in $`z`$ becomes large at lower redshifts due to uncertainty in the phototometric redshifts. The final catalog size is over 2 million galaxies.

### 2.3. Rest-frame Selected Subsamples

To analyze how angular clustering changes with luminosity in Section 3, the volume limited sample is divided into 3 absolute magnitude bins. These subsamples represented by $`M_1`$, $`M_2`$ and $`M_3`$ have limiting absolute magnitudes $`M_r^{}>21`$, $`21>M_r^{}>22`$ and $`22>M_r^{}>23`$ respectively. The size of these subsets decrease by approximately a factor of two as a function of increasing luminosity.

For the type dependent selection we utilize the continuous spectral type parameter $`t`$ from the photometric redshift estimation. This essentially encodes the rest-frame colors of the galaxies as there is a direct one-to-one mapping between the type $`t`$ and the spectral energy distribution (SED) of the galaxy. The value $`t=0`$ represents a template galaxy that is as red as the elliptical spectrum of Coleman, Wu & Weedman (1980); as $`t`$ increases the galaxy type becomes progressively later. In Section 3, we subdivide by spectral class breaking the luminosity classes into four subgroups (each with comparable numbers of galaxies). The cuts in the spectral type parameter $`t`$ from red to blue are $`t<0.02`$, $`0.02<t<0.3`$, $`0.3<t<0.65`$ and $`t>0.65`$. The cuts are defined as the $`T_1`$, $`T_2`$, $`T_3`$ and $`T_4`$ subsamples respectively. Our selection is motivated by the spectral energy distributions of Coleman, Wu & Weedman (1980). The first class consists of galaxies with SEDs similar to the CWW elliptical template (Ell), the second, third and fourth classes contain a broader distribution of galaxy types approximately corresponding to Sbc, Scd and Irr types, respectively. The distribution of types and our classification are shown in Figure 3 and 4.

## 3. The Angular Correlation Function

The properties of the angular correlation function and the estimators used to measure it from photometric catalogs have been extensively discussed in the astronomical literature (Kerscher et al., 2000). The probability of finding galaxy within a solid angle $`\delta \mathrm{\Omega }`$ on the celestial plane of the sky at distance $`\theta `$ from a randomly chosen object is given by (Peebles, 1980)

$$\delta P=n[1+w(\theta )]\delta \mathrm{\Omega },$$
(1)
where $`n`$ is the mean number of objects per unit solid angle. The angular two-point correlation function $`w(\theta )`$ basically gives the excess probability of finding an object compared to a uniform Poisson random point process.

Traditionally the observations are compared to random catalogs that match the geometry of the survey. The computation usually consists of counting pairs of objects drawn from the actual and random catalogs and applying a minimum variance estimator such as that defined by Landy & Szalay (1993) or Hamilton (1993). In this study, we use the Landy-Szalay estimator (Landy & Szalay, 1993) as

$$w_{LS}=\frac{\mathrm{DD}2\mathrm{D}\mathrm{R}+\mathrm{RR}}{\mathrm{RR}},$$
(2)
where DD, RR and DR represent a count of the data-data, random-random and data-random pairs with $`\theta `$ angular separation summed over the entire survey area.

### 3.1. Estimating $`w(\theta )`$ with a Fast Fourier Transform

Even though distances on the sky are easy to compute mathematically, measuring the correlation function is not a trivial task, especially when it comes to large surveys. It’s easy to see that any naive algorithm implementing the estimator in Eq. (2) scales with the square of the number of objects $`N`$ in the survey, $`𝒪(N^2)`$. It, therefore, becomes progressively more expensive in computational power to apply these techniques to data the size of the SDSS samples.

We propose a novel method to estimate $`w(\theta )`$ using the Fast Fourier Transform (e.g. in Press et al., 1992, hereafter FFT), which scales as $`𝒪(N\mathrm{log}N)`$, a significant improvement over the naive approach. The principle behind this is to group all galaxies into small cells within a grid and analyze this matrix instead of the point catalog. The implementation of this FFT approach called eSpICE is a Euclidean version of SpICE by Szapudi, Prunet & Colombi (2001). A discussion on its properties and scalings is given elsewhere (Szapudi et al., 2003). Here we present an application of eSpICE to the SDSS galaxy angular clustering and limit our discussion of the algorithm to only an outline of how the method works.

To apply a standard FFT analysis to the problem we operate in Euclidean space. In other words, Euclidean distances are computed instead of the correct angular separation. Given that we only examine separations of less than 2 degrees, the maximum relative distance error introduced by the use of the Euclidean approximation is only 0.00005. This has a negligible effect on the accuracy of our analysis, and the use of Euclidean distances increases the speed of the algorithms significantly. We use SDSS survey coordinates ($`\mu ,\nu `$) to define the position of an object within a stripe. These coordinates are defined locally for each stripe. As the individual stripes comprise great circles on a sphere the survey coordinates along these great circles are close to Euclidean (i.e. in this coordinate system every stripe looks as if it were equatorial). In fact for equatorial stripes 10 and 82, the ($`\mu ,\nu `$) coordinates are the same as (RA, Dec). Beyond the above data grid of galaxies, we need to describe the geometry of the survey area. This is done by a second grid, with the same dimensions that describes the survey boundaries. We call this the window grid because it is 1 if the pixel is entirely inside the boundaries and 0 otherwise. In this way we can also incorporate arbitrarily complex masks (e.g. for excluding regions around bright stars) as they are simply applied to the window grid. The data and window grids are then padded with zeros up to the maximum angular scale to avoid aliasing in the Fast Fourier Transform. Finally, eSpICE is used to calculate the two-point correlation function directly.

Our approach is to apply this algorithm to one stripe at a time. The data and window matrices are constructed by querying the SQL database to select the appropriate galaxies and masks. From these data we calculate the two-dimensional angular correlation function for all eight stripes. Figure 5 shows the 2D correlation function for all stripes used in this analysis. We find that the 2D correlation function is an extremely sensitive diagnostic tool for identifying systematics within the photometric data. The correlation function is expected to be isotropic. Artifacts within the data or survey geometry distort this symmetry. This is seen, to varying degrees, in each of the eight 2D correlations functions. We find that within the 2D correlation function there is an elongated streak, at zero lag, along the scan direction which has structure on scales in excess of a few degrees. This arises due to errors in the flat-field vector.

In a drift-scan survey the flat field is a one-dimensional vector (orthogonal to the direction of the scan). Errors within the flat field tend, therefore, to be correlated along the scan direction (i.e. along the columns of a stripe). This effect is seen within all of the individual stripe correlation functions. As it is, by definition, a zero lag effect we can exclude it from the analysis by censoring this region of the 2D correlation function before azimuthally averaging the signal to get a 1D correlation function. We note, however, that even if we do not censor the data to remove this effect the result of averaging the correlation function azimuthally (and the fact that this elongated streak only affects a small fraction of the 2D correlation function) the results discussed in the following sections are not affected.

In total, the computation of the one-dimensional angular correlation function $`w(\theta )`$ takes less than 3 minutes for a stripe, which is several orders of magnitude faster than traditional two-point estimators. Having computed the correlation function for all stripes, we co-add the results by properly weighting with the number of galaxy pairs. At the same time the covariance of the signal is also estimated from the scatter between the different stripes Figure 6 shows the covariance matrix for the volume limited sample ($`z<0.3`$). The image of the matrix is normalized so that the diagonal elements are always white and the (off-diagonal) gray scale values represent the correlation between the bins. Because of the logarithmic sampling of $`w(\theta )`$, the neighboring bins are farther apart and hence correlate progressively less as we go to larger scales and thus to larger bins.

### 3.2. The Clustering Scale Length

With the measured correlation amplitudes as function of angular separation in hand, we can obtain a parametric form of the scaling. The amplitude and power of the correlation are calculated by fitting the usual formula,

$$w(\theta )=A_w\left(\frac{\theta }{\theta _0}\right)^\delta .$$
(3)
Normalizing $`\theta `$ by $`\theta _0=0.1^{}`$ enables the amplitude $`A_w`$ be directly compared with the figures showing $`w(\theta )`$ measurements; $`A_w`$ is essentially the value of the correlation function at $`\theta _0`$, since $`A_ww(\theta _0)`$. The parameters $`A_w`$ and $`\delta `$ are estimated by minimizing the cost function

$$\chi ^2(A_w,\delta )=\frac{1}{N2}\underset{i,j=1}{\overset{N}{}}\mathrm{\Delta }w_iC_{ij}^1\mathrm{\Delta }w_j,$$
(4)
where $`\mathrm{\Delta }w_i=w_iw(\theta _i|A_w,\delta )`$, $`N`$ is the number of bins and $`C_{ij}`$ is the covariance matrix. Although, $`\chi ^2(A_w,\delta )`$ is not quadratic in the parameter $`\delta `$, it is in $`A_w`$, thus the equation

$$\frac{\chi ^2}{A_w}=0$$
(5)
can be solved analytically for $`A_w`$ to reduce the dimensionality of the problem. We then use a method by Brent (1973) (see Press et al., 1992) to search for the optimal value of $`\delta `$. The $`\chi ^2`$ fit also provides information about the covariances of the estimated parameters that are shown as error ellipses in the next section.

Figure 7 illustrates the angular clustering in our fiducial sample of all galaxies within the volume limited sample. In the left panel, the correlation function $`w(\theta )`$ is shown along with the best power-law fit. The range in angular separation varies from 1 arcminute to 2 degrees. At a mean redshift for the volume-limited sample of $`z=0.2`$, $`2^{}`$ corresponds to $`17h^1\mathrm{Mpc}`$. The error bars on the measurements are computed as one over the square root of the diagonal elements of the covariance matrix (e.g. shown in Figure 6). In the right panel of Figure 7, we plot the best fit parameters and their error ellipses. For the full volume-limited sample the best fit to the data has a slope of $`\delta =0.84\pm 0.02`$ with an amplitude of $`A_w=0.078\pm 0.001`$.

### 3.3. Clustering as a Function of Luminosity

The angular clustering of the galaxies in the three luminosity bins described in Section 3 are compared in Figure 8. As noted earlier, the left panel shows the correlation function, and the right panel gives the parameters of the power-law fits. The first section (I) of Table 1 presents the values and errors on these measured parameters. As expected, the more luminous galaxies are clustered more strongly: the amplitudes are roughly larger by a factor of 1.5 from one sample to the next and are measured to be $`0.061\pm 0.001`$, $`0.092\pm 0.002`$ and $`0.138\pm 0.004`$. The slope of the correlation functions are consistent for all luminosity bins and with the fiducial value derived for the entire volume (the estimated slope parameters scatter around $`\delta =0.84`$).

### 3.4. Bimodality of $`w(\theta )`$ as a Function of Rest-frame Color

The type dependence of the correlation function is not as simple as that found for the luminosity classes. It is well known that different types of galaxies have different clustering behavior (Giovanelli et al., 1986). Red elliptical galaxies are more likely to be found in higher density regions than spirals. Here, we study the evolution of the angular correlation function with spectral type in two absolute magnitude ranges, $`M_r^{}>21`$ and $`21>M_r^{}>23`$. Both of these yield approximately 1 million galaxies within our volume. Figure 9 shows how the clustering changes as a function of spectral type in the lower luminosity sample. The 2 reddest classes of the galaxy population ($`T_1`$ and $`T_2`$) are essentially indistinguishable, their clustering is significantly stronger than for the other classes or our fiducial results. The bluest 2 classes ($`T_3`$ and $`T_4`$) have approximately the same power-law exponents but with different amplitudes. Section II of Table 1 summarizes the results of parameter fitting, which is also seen in the right panel of Figure 9. The higher luminosity classes show the same basic trends but with a stronger correlation amplitude. The results of the type dependence of power law fits to the high luminosity class are given in Figure 10 and Table 1.

## 4. The Correlation Length for Galaxies

### 4.1. Limber’s Equation

From the angular clustering we can derive the spatial correlation length, $`r_0`$, given the redshift distribution of the data (Peebles, 1980). This is accomplished by integrating over the comoving coordinate along two lines of sight $`r_1`$ and $`r_2`$, separated by the angle $`\theta `$, to calculate the projected angular correlations from the real-space correlation function $`\xi (r)=(r/r_0)^\gamma `$,

$$w(\theta )=\frac{r_1^2\mathrm{\Phi }(r_1)dr_1}{F(r_1)}\frac{r_2^2\mathrm{\Phi }(r_2)dr_2}{F(r_2)}\xi _{12},$$
(6)
where $`\mathrm{\Phi }(r)`$ is the selection function and the factor $`F(r)=1kr^2`$ accounts for the curvature of space.<sup>1</sup><sup>1</sup>1In our case $`F(r)`$ would be constant $`1`$ because the parameters $`\mathrm{\Omega }_M=0.3`$ and $`\mathrm{\Omega }_\mathrm{\Lambda }=0.7`$ assume a flat universe, however, we will do the integral in redshift (see later).  In the small angle approximation, the correlation function $`\xi _{12}=\xi (r_1^2+r_2^22r_1r_2\mathrm{cos}\theta )`$ becomes $`\xi _{12}=\left(r/r_0\right)^\gamma \left(\theta ^2+u^2/r^2\right)^{\gamma /2}`$ using variables $`2r=r_1+r_2`$ and $`u=r_1r_2`$. This simplifies the above integral, because one can separate out the part that has $`u`$ and arrive at

$$w(\theta )=r_0^\gamma H_\gamma \theta ^{1\gamma }\frac{r^{5\gamma }\mathrm{\Phi }^2(r)}{F^2(r)}𝑑r,$$
(7)
where

$$H_\gamma =\frac{\mathrm{\Gamma }(\frac{1}{2})\mathrm{\Gamma }(\frac{\gamma 1}{2})}{\mathrm{\Gamma }(\frac{\gamma }{2})}.$$
(8)
We can rewrite the integral using redshift $`z`$ and its distribution. Substituting $`r^2\frac{\mathrm{\Phi }}{F}dr=\frac{dN}{dz}dz`$ one gets the final integral of

$$w(\theta )=r_0^\gamma H_\gamma \theta ^{1\gamma }_0^{\mathrm{}}r^{1\gamma }\left(\frac{dN}{dz}\right)^2\left(\frac{dr}{dz}\right)^1𝑑z$$
(9)
that may be compared to the measured quantities directly. We have assumed that the selection function and the redshift distribution are normalized to $`\frac{dN}{dz}𝑑z=1`$, which can be easily achieved numerically.

### 4.2. Estimating $`dN/dz`$

To determine the correlation length, we need to estimate $`dN/dz`$. Given a photometric redshift and its random error, there is a conditional probability of having a galaxy at a certain true redshift. In addition, we need to incorporate the apparent magnitude limit and photometric redshift selection criteria in the estimate of the redshift distribution. The real redshift $`z`$ and the photometric redshift $`s`$ of a galaxy are different. We assume that $`s=z+\nu `$ where $`\nu `$ is the error and drawn from a normal distribution. Thus

$$P(s|z)=\frac{1}{\sqrt{2\pi \sigma ^2}}\mathrm{exp}\left\{\frac{(sz)^2}{2\sigma ^2}\right\},$$
(10)
where $`\sigma `$ determines the precision of the estimates. We need the inverse: what is the true redshift, given the photometric estimate. This may be obtained from Bayes’ theorem,

$$P(z|s)=\frac{P(z)P(s|z)}{P(s)},$$
(11)
where $`P(z)`$ is the true redshift distribution calculated from the LF that also depends on the apparent magnitude cuts.

Our volume limited sample is selected by a window function of photometric redshifts,

$$W(s)=\{\begin{array}{cc}1\hfill & \text{if }s\text{ between 0.1 and 0.3,}\hfill \\ 0\hfill & \text{otherwise.}\hfill \end{array}$$
(12)
The conditional probability of having a galaxy in a sample selected by this window function is calculated by the integral over the distribution $`P(s)`$,

$`P(z|W)`$ $`={\displaystyle 𝑑sP(s)W(s)P(z|s)}`$ (13)
$`=P(z){\displaystyle 𝑑sW(s)P(s|z)}`$ (14)
$`=P(z)W_\sigma (z),`$ (15)
where $`W_\sigma (z)`$ is the photometric redshift selection function convolved with the photometric redshift uncertainty.

The precision of photometric redshifts is a strong function of the apparent magnitude. From the photometric redshift catalog, we compute the mean redshift errors $`\sigma _i`$ and the galaxy counts $`n_i`$ (for proper weighting) in $`\mathrm{\Delta }r^{}=0.2`$ wide magnitude bins. Using the LF by Blanton et al. (2003), we derive the redshift distributions $`P_i(z)`$ for the same magnitude bins up to $`r^{}=21`$, which is the limiting magnitude in the sample. The final redshift distribution is the weighted average of these probabilities,

$$\overline{P}(z|W)\underset{i}{}n_iP_i(z)W_{\sigma _i}(z).$$
(16)
The top panel of Figure 11 illustrates the redshift distribution derived from the LF for galaxies with $`18<r^{}<18.2`$ and a smoothed window $`W_\sigma (z)`$ with $`\sigma =0.05`$. The middle panel shows the effective contribution of this magnitude bin to the final $`dN/dz`$, which is plotted in the last panel.

### 4.3. Results for $`r_0`$

Applying the above probabilistic redshift distribution and substituting the power-law fits derived earlier, we find a correlation length for the full volume limited sample of $`5.77h^1\text{Mpc}`$. Estimating the uncertainty of this value is not trivial. The statistical errors on the parametric fits to $`w(\theta )`$ are known and may be used to estimate the errors on the correlation length. These are estimated by calculating the scatter of the predicted $`r_0`$ measurements for 20,000 Monte-Carlo realizations of the fitting parameters, $`A_w`$ and $`\delta `$, based on their covariance matrix. For the full volume limited sample, we obtain $`\mathrm{\Delta }_{\text{rms}}=0.05h^1\text{Mpc}`$. However, the uncertainty on $`r_0`$ is also affected by the uncertainty in the redshift distribution. The primary source of change in the $`dN/dz`$ is the uncertainty in the LF parameters. We compute the partial derivatives $`(r_0/\alpha )`$, $`(r_0/M_{})`$ and $`(r_0/Q)`$ numerically and propagate the quoted errors of Blanton et al. (2003). We find that the evolutionary parameter $`Q`$ makes the largest difference, the $`r_0`$ errors from the uncertainty of $`M_{}`$ and $`\alpha `$ are negligible. For the fiducial $`r_0`$ value, we estimate an error of $`\mathrm{\Delta }_{\text{LF}}=0.09h^1\text{Mpc}`$. The dependence of $`\sigma `$ on the apparent magnitude was determined empirically using the actual measurements, thus the results are not affected by Malmquist bias.

In Figure 12, the correlation length is plotted as a function of luminosity (left panel) and SED type (right panel). The relation between $`r_0`$ and luminosity and spectral type are consistent with that observed directly from the angular data and from measures of the clustering length from spectroscopic surveys (Zehavi et al., 2002). Over an absolute magnitude range of $`19.97`$ to $`22`$, $`r_0`$ increases with luminosity from a value of $`5.04\pm 0.09`$ to $`7.87\pm 0.24h^1\text{Mpc}`$, which is consistent with the increase observed by Zehavi et al. (2003).

The color dependence of $`r_0`$ (i.e. clustering as a function of spectral type) shows the expected increase in clustering length for early type galaxies. The values of the correlation length for the spectral type subsamples are given in Table 1. In the lower luminosity bin red, $`T_1`$, galaxies have a correlation length of $`6.59\pm 0.17h^1\text{Mpc}`$ for the $`M_r^{}>21`$ sample and the bluest galaxies, $`T_4`$, have a correlation length of $`4.51\pm 0.19h^1\text{Mpc}`$. The trend in this relation is again consistent with the observed dependence of the correlation length as a function of spectral and morphological type (Giovanelli et al., 1986).

## 5. Discussion

The interesting aspect of the luminosity and color dependent clustering is not that the observed clustering length scales with luminosity and color as this has been demonstrated from many different surveys. It arises from how the shape of the correlation function depends on luminosity and color. It is remarkable that the luminosity simply effects the amplitude of the correlation function and not the slope, whereas the type selection affects both the slope and amplitude. This is particularly intriguing when we note that we would expect an intrinsic correlation between the luminosity and spectral type of a galaxy.

Observationally, early type galaxies tend to reside in clusters of galaxies whereas later type galaxies are more often found in the field. We would expect, therefore, that early type galaxy samples would have more small separation pairs than samples selected for late type galaxies and that their resulting correlation functions would be steeper. This is consistent with the data except that we would expect there to be a smooth transition from early to late type galaxies and that the correlation function slope should smoothly change from a steep value of $`0.96`$ for the early types to the more shallow value of $`0.68`$ for the late types. What we observe, however, is that red galaxies (types $`T_1`$ and $`T_2`$) have a common slope of $`0.1`$ and blue galaxies (types $`T_3`$ and $`T_4`$) have a slope of $`0.07`$ (i.e. there does not appear to be a smooth transition).

We can, however, explain this behavior if we consider a simple model for the distribution of galaxy types. From Figure 3 we see that the type histogram for the galaxies is almost bimodal (Strateva et al., 2001; Hogg et al., 2003) with the distribution being well fit by two Gaussians. For simplicity we will denote these subclasses as “red” and “blue”. If the “red” and “blue” populations have distinct correlation functions (i.e. with different slopes), then any observed correlation function should simply come from mixing these populations. As we change the mix of “red” and “blue” galaxies then the resulting slope of the correlation function will also change. This is exactly what we observe with the correlation function as we move from the $`T_2`$ to $`T_3`$ selections.

If we selected galaxies only from either the “red” or “blue” sub-populations we would expect no change in the correlation function slope as all of the “red” or “blue” galaxies have a common correlation function. Again, this what we observe from the data. The slopes of the correlation functions for the $`T_1`$ and $`T_2`$ red samples are identical as are the correlation functions for the blue $`T_3`$ and $`T_4`$ samples. In reality, the color selection that we applied to the SDSS data (types $`T_1`$ through $`T_4`$) was not chosen to optimally separate two distinct populations of galaxies but rather to provide a simple subdivision of galaxies based on the CWW spectral energy distributions. We might expect there to remain some population mixing in our $`T_1`$, $`T_2`$, $`T_3`$ and $`T_4`$ color cuts. We observe this effect where the amplitude of the correlation function for the $`T_3`$ and $`T_4`$ classes are close but not identical; this would imply that the $`T_3`$ still contains a subset of the “red” galaxies.

The luminosity dependence can be explained if we note that the luminosity functions of the “red” and “blue” classes are identical for magnitudes brighter than $`M_r^{}=20`$ (Baldry et al., 2003). They deviate only for the faint end of the luminosity function (“blue” galaxies having a steeper faint end slope). Varying the luminosity cuts should not change the mix of the galaxy populations (unless we sample galaxies with $`M_r^{}>20`$). We would, therefore, expect the shape of the correlation function to be independent of the luminosity cuts (as is found from the data). Given this hypothesis if we selected a volume limited sample for galaxies with luminosities less that $`M_r^{}=20`$ we would expect to find a dependence on the slope with luminosity.

It is, therefore, remarkable that with such a simple model for the distribution of galaxies (i.e. just two classes with differing correlation functions) we can qualitatively describe the behavior of the correlation functions with color and luminosity. What is difficult to understand is why there would be a simple scaling of the amplitude of the correlation function with intrinsic luminosity as the spatial scales we sample are in the non-linear regime (i.e. a simple linear bias model is not necessarily appropriate). Identifying the physical mechanism that could give rise to the luminosity and type dependent bias that we observe remains an open question.

We would like to thank Michael Strauss and Michael Blanton for their valuable comments. AS was supported by NSF grants AST-9802 980, KDI/PHY-9980044, NASA LTSA NAG-53503 and NASA NAG-58590. TB acknowledges support from NSF KDI/PHY-9980044. IC, AS and TB acknowledge partial support from MTA-NSF grant no. 124 and the Hungarian National Scientific Reasearch Foundation (OTKA) grant T030836. AJC acknowledges partial support from an NSF career award AST99-84924 an ITR NSF award 0121671. RS acknowledges partial support form ITR NSF award 0121671. IS was supported by NASA through AISR grants NAG5-10750, NAG5-11996, and ATP grant NASA NAG5-12101 as well as by NSF grants AST02-06243. Funding for the creation and distribution of the SDSS Archive has been provided by the Alfred P. Sloan Foundation, the Participating Institutions, the National Aeronautics and Space Administration, the National Science Foundation, the U.S. Department of Energy, the Japanese Monbukagakusho, and the Max Planck Society. The SDSS Web site is http://www.sdss.org/. The SDSS is managed by the Astrophysical Research Consortium (ARC) for the Participating Institutions. The Participating Institutions are The University of Chicago, Fermilab, the Institute for Advanced Study, the Japan Participation Group, The Johns Hopkins University, Los Alamos National Laboratory, the Max-Planck-Institute for Astronomy (MPIA), the Max-Planck-Institute for Astrophysics (MPA), New Mexico State University, Princeton University, the United States Naval Observatory, University of Pittsburgh and the University of Washington.


format: text

text:
# Molecular dynamic simulation of the fluid concentration fluctuations in two dimensional porous medium

## Conclusion

In our computer experiments the decay of the concentration of the particles in separate pore is well described by exponents. This leads to the validity of relaxation theories . The time of relaxation does not depend on the size of the system. The increasing of the repelling forces between the fluid and matrix particles tends to increase the time of relaxation of fluid particles in pores, what corresponds to the growing of filtration resistance.

## Figure captions

Fig.1. The time dependence of the particles number in pore (1,1). The density of the particles in pore (1,1) more than in others 3 times (curve 1), 15 times (curve 2), 81 times (curve 3).

Fig.2. The time dependence of the number of particles in pores (1,2), (2,2) and (1,1).

Fig.3. The dependence of the characteristic time of particles density decays in pore (1,1) from effective width of the leaking field.


text:
# References

A note on the thermodynamics of gravitational radiation

T. Padmanabhana111e-mail address: nabhan@iucaa.ernet.in and T. P. Singhb222e-mail address: tpsingh@tifr.res.in

aInter-University Center for Astronomy and Astrophysics
Post Bag 4, Ganeshkhind, Pune 411007, India
bTata Institute of Fundamental Research,
Homi Bhabha Road, Mumbai 400 005, India.

## Abstract

It is shown that linearized gravitational radiation confined in a cavity can achieve thermal equilibrium if the mean density of the radiation and the size of the cavity satisfy certain constraints.

Is it possible, in principle, to confine gravitational radiation in a box, for a time long enough for it to achieve thermal equilibrium? This issue has been investigated before by various authors , not all of whom have arrived at the same conclusion. The purpose of this note is to pursue a particular line of thought which does not seem to have been looked at in earlier work. We argue that it is possible to thermalize gravitational radiation, starting from a non-equilibrium configuration, so long as certain constraints (to be discussed below) are satisfied.

Let us begin by summarizing the arguments given by Smolin . He showed, by considering various mechanisms of absorption of radiation by matter, that no realistic material is an efficient absorber of gravitational radiation, except over a narrow band of frequencies. These mechanisms include absorption of classical gravitational radiation by classical matter, ionization of atoms by gravitons, and phonon excitation. He concluded that a state of thermal equilibrium between radiation and matter cannot be reached in a finite time.

Garfinkle and Wald  presented a counter-example to the result by using a shell of charged matter balanced just outside its Schwarzschild radius. However Dell  later argued that such a shell is unstable in the presence of electromagnetic radiation.

Whatever the nature of the confining cavity, the absorption coefficient $`f`$ will always be smaller than one. If the box is made thick enough to raise $`f`$ to one, it undergoes gravitational collapse.

However, it is not absolutely essential, for the purpose of attaining equilibrium, to have $`f`$ close to unity. For any cavity, it is possible to define a leakage time-scale $`t_L`$ over which, say, a fraction $`1/e`$ of the initial radiation escapes. Also, for an initial configuration $`\rho (\nu )`$ \- the energy-density of gravitational radiation in the cavity - an equilibrium time-scale $`t_E`$ can be defined, over which the radiation may achieve equilibrium. What needs to be settled is whether $`t_E`$ can be smaller than $`t_L`$, for a general initial configuration $`\rho (\nu )`$. We argue below that this is possible, and $`t_E<t_L`$ need not imply gravitational collapse of the cavity or of the radiation inside it. Moreover, we have no need to concern ourselves with any specific mechanism for confining radiation.

Our argument is analogous to the standard discussion of Brownian motion of a molecule in a fluid, where one shows that as a consequence of the random collisions and viscous drag, the molecule attains a Maxwellian velocity distribution. A variant of this argument was used by Einstein  to derive the Planckian distribution for electromagnetic radiation in equilibrium in a cavity. The argument went as follows. He assumed that the radiation was in equilibrium with a molecule (say, with two energy levels). The molecule absorbs and re-emits at a characteristic frequency $`\nu _0`$. The momentum transfers to the molecule take place through random collisions with the quanta of radiation, and through the systematic drag force it experiences because it sees a Doppler-shifted radiation while in motion. Einstein showed that if the molecule has a Maxwellian velocity distribution, then the exchange of momenta with radiation preserves this distribution only if the energy density of the radiation obeys the Planck law.

In his paper, Einstein concerned himself with the equilibrium situation. The discussion of momentum transfers to the molecule assumes greater significance in a study of $`\mathrm{𝑎𝑝𝑝𝑟𝑜𝑎𝑐ℎ}`$ to equilibrium. Now, there is a systematic evolution of the r.m.s. velocity of the molecule, which may be approximately described by the Langevin equation

$$M\frac{d𝐯}{dt}=\alpha 𝐯+\beta (t),$$
(1)
$`\alpha `$ being the drag coefficient and $`\beta `$ the fluctuating force ($`\alpha `$ is positive). Standard discussions  of Brownian motion show that the molecule attains a Maxwellian distribution over a time-scale

$$t_E\frac{M}{\alpha }.$$
(2)
In particular, the ‘fluid’ pushing the molecule around could be electromagnetic radiation, and (2) will continue to hold. \[The assumptions which go into using the Langevin equation as well as subtleties related to the dependence of the drag coefficient $`\alpha `$ deserve discussion; this is presented in the last section since they are unrelated to the main theme of this paper.\] Once the Maxwellian distribution has been reached for the molecule, Einstein’s argument can be applied to arrive at Planck’s law for the radiation.

We now return to the case of linearized gravitational radiation in a cavity, interacting with a two-level atom.

To begin with, we compare this system with the one consisting of an atom in electromagnetic radiation. The only difference between the two kinds of radiation is that the former is described by a spin-2 field, while the latter is spin one. As regards the two-level atom it is sufficient for us to assume that it has a non-zero cross section $`\sigma `$ for the absorption of gravitons at a frequency $`\nu _0`$ in its rest frame. The next point concerns the walls of the confining cavity. For the case of electromagnetic radiation it is safe to assume that one can construct walls with arbitrarily high absorptivity, and that the leakage time $`t_L`$ can be pushed to infinity. For the gravitational case, we assume that the cavity has an absorption coefficient $`f(\nu )1\eta (\nu )`$, so that for a spherical cavity of radius $`R`$

$$\frac{d(\nu )/dt}{(\nu )}=\frac{\eta (\nu )\rho (\nu ,t)c.4\pi R^2}{4\pi R^3\rho (\nu ,t)/3}$$
(3)
where $`(\nu )`$ is the energy in the cavity at frequency $`\nu `$.

This may be used to define a leakage time-scale

$$t_L(\nu )R/c\eta (\nu ),$$
(4)
and a frequency averaged scale $`\overline{t}_L=R/c\eta `$, $`\eta `$ being the mean ‘leakage’ coefficient. Clearly, any reference to the thickness or density of the absorbing material will be through $`\eta (\nu )`$. As expected, leakage takes more time for larger cavities.

Consider next, the motion of the atom of mass $`M`$, when it is pushed about by the gravitational radiation, which has a time-dependent spectrum $`\rho (\nu ,t)`$. We assume that this motion can be described by the Langevin equation (2). This implies that the velocity of the atom will acquire a Maxwellian distribution over a time-scale $`M/\alpha (\nu _0)`$, irrespective of its initial velocity. When this has happened for all possible frequencies, we can say that the radiation has also reached equilibrium. The time-scale $`t_E`$ should thus be defined as $`M/\alpha _{min}`$, where $`\alpha _{min}`$ is the minimum value of $`\alpha (\nu )`$, obtained by varying $`\nu `$. A major task is to relate $`\alpha _{min}`$ to an initial distribution $`\rho _0(\nu )`$ and we address this question below. Right now, we compare $`t_E`$ with the average time-scale $`\overline{t}_L`$ defined earlier. Requiring $`t_E<\overline{t}_L`$ gives

$$R>Mc\eta /\alpha _{min}$$
(5)
as the necessary condition for equilibrium.

We must ensure that the box is not dense enough to undergo gravitational collapse. In other words, $`\eta `$ is assumed to be close to one. We also need to ensure that when $`R`$ is chosen as in (5), the radiation in the cavity does not collapse to form a black hole. This requires

$$2G.\frac{4\pi }{3}R^3\frac{\overline{\rho }}{c^2}<Rc^2$$
(6)
i.e., (ignoring numerical coefficients),

$$R^2<c^4/G\overline{\rho }$$
(7)
$`\overline{\rho }`$ being the mean density of radiation. The compatibility of (5) and (7) requires

$$\frac{M^2c^2\eta ^2}{\alpha _{min}^2}<R^2<c^4/G\overline{\rho }$$
(8)
which may be taken to imply that $`\rho `$ should satisfy

$$\overline{\rho }<\alpha _{min}^2c^2/GM^2\eta ^2.$$
(9)
We interpret the equations (5)-(9) as follows. If the mean density satisfies the bound (9), then the radiation can be made to achieve equilibrium by choosing $`R`$ within the bounds given by (8). The situation is different from the one for electromagnetism because $`\eta `$ cannot be made one for electromagnetic radiation, and because we do not want $`t_E<\overline{t}_L`$ to imply gravitational collapse. Such a constraint on $`\eta `$ is only to be expected, because unlike for the electromagnetic case, the principle of equivalence a priori fixes the ‘charge to mass’ ratio for gravity.

It is evident that perfect absorption $`(f=1)`$ is a sufficient but not a necessary condition for thermalization. This allows the construction of the above mechanism for achieving equilibrium. On the other hand, Smolin was mainly concerned with measurement of the quantum state of the gravitational field, and that certainly requires $`f1`$. Thus our results should be considered complimentary, rather than contradictory, to those of Smolin.

We end this discussion by a comment on the coefficient $`\alpha _{min}`$ in (9). Near equilibrium, $`\alpha _{min}(\nu )`$ is proportional to $`\nu `$, and vanishes for $`\nu =0`$, which only means that no momentum is absorbed when the energy levels coincide. So it is reasonable to define $`\alpha _{min}`$ as $`\alpha (\nu _{min})`$ where $`\nu _{min}`$ is the smallest frequency for which we are interested in studying the approach to equilibrium. Clearly, this issue arises for electromagnetic radiation also, and is not peculiar to gravity.

For the equilibrium case, $`\alpha (\nu )`$ is simple to calculate; but for an arbitrary configuration an exact calculation is not possible. A very crude estimate for $`\alpha _{min}`$ may be obtained as follows. Consider an atom which is moving with a speed $`V`$ in an isotropic radiation with mean density $`\overline{\rho }`$. The average momentum transferred to the atom per unit time will be of the order

$$\frac{\overline{\rho }}{c}.\frac{V}{c}.c\sigma =\overline{\rho }\sigma .\frac{V}{c}$$
(10)
$`\sigma `$ being the mean cross-section of the atoms to absorption of gravitons. This gives $`\alpha _{min}\overline{\rho }\sigma /c`$, and from (9) it follows that

$$\overline{\rho }>GM^2\eta ^2/\sigma ^2.$$
(11)
Numerical estimates of the allowed ranges for $`R`$ and $`\overline{\rho }`$ will depend on $`\sigma `$, and on the nature of the confining box (through $`\eta `$). In principle, however, gravitational radiation with $`\overline{\rho }`$ exceeding the bound in (11) will reach equilibrium if placed in a suitable box.

Having concluded the central part of our discussion, we would like to comment on several aspects of the analysis which deserve careful understanding.

We begin with the key technique used in the analysis, viz. the Langevin equation and its validity. Langevin equation arises in many situations dealing with stochastic forces (see e.g., for a discussion of its validity and comparison with Fokker-Planck equation) and is a good assumption if the scattering process is uncorrelated and homogeneous. In this particular case we assume that the system is interacting with a bath of gravitons which can be described classically as a stochastic background of gravitational waves. While one could think of more complicated scenarios, the description based on Langevin equation is adequate for our purpose. Of course, Langevin, equation is never “fundamental” and is more like a phenomenological equation such as, say, Hooke’s law. But since we are discussing a specific situation well within the domain of its validity, our results are trustworthy.

The second, and somewhat more subtle issue is whether $`\alpha `$ depends on $`M`$ and if so, how. This is interesting because, if $`\alpha `$ is proportional to $`M`$ with, say $`\alpha =M\alpha _0(\nu )`$ where $`\alpha _0(\nu )`$ is independent of $`M`$, then some of our bounds (which depend on $`\alpha /M`$) will be independent of $`M`$. In fact, one may naively imagine that principle of equivalence will require $`\alpha M`$ making final results independent of $`M`$. In reality, however, $`\alpha M^2`$ and there are different ways of seeing this. The first, and the simplest, is to note that in the electromagnetic case $`\alpha q^2`$ (where $`q`$ is the charge) and not $`\alpha q`$; this is obvious since Langevin equation will treat both positive and negative charged particles at the same footing and the sign of $`\alpha `$ cannot depend on the sign of $`q`$. By analogy, the $`\alpha `$ for gravity should scale as $`(GM)^2`$. The second argument is based on the result in equation (10) where we showed that $`\alpha \sigma `$. Now the lowest order scattering of gravitational waves by matter consists of matter oscillating under the infulence of the incident gravitational wave and radiating out gravitational waves due to its motion. This is completely analogous to the Thomson scattering of electromagnetic radiation and the scattering crosssection for this process can be worked out in a straight forward manner (Padmanabhan, 1978; unpublished). The result is as expected with $`\sigma (GM/c^2)^2;`$ i.e, the Schwarzschild radius acts as the analogue of classical electron radius. Once again, we find that $`\alpha \sigma M^2`$. Finally, note that the averaging of the stochastic force acting on the particle will involve the two point correlation function of the force at the lowest order. That is why $`\alpha `$ will scale in proportion with the power spectrum \[energy density\] of the stochastic background (quadratic in the gravitational wave amplitude) and — by the same token — will scale as $`M^2`$. This imples that the time taken to reach equilibrium depends on the mass $`M`$ of the atom. One should not use principle of equivalence naively to argue otherwise; principle of equivalence has no direct role to play in this issue. This is the same result as in the treatment of Brownian motion via the Langevin equation. As a result, the various bounds obtained above depend on $`M`$. Of course, the mass dependence of $`\alpha `$ is not germane to the main issue we are addressing. Our main concern here is to demonstrate that equilibrium can in principle be reached, even though it implies that the bounds will depend on $`M`$.

It may appear that Eqns. (9) and (11) give conflicting bounds, the former being an upper bound and the latter a lower bound, on the density. What is happening here of course is that one does not in general know what $`\alpha _{min}`$ is, and (11) results from using a very rough estimate for $`\alpha _{min}`$. It is also important to note that, we have tried to obtain fairly general conclusions with most plausible assumptions. In specific contexts, one may be able to strengthen the bounds or supercede them with other considerations but these will not be generic. For example, we have treated gravitational waves as a small perturbation on a flat background. While the gravity waves being a perturbation defines the context of the problem, one can think of the effects of background being curved. If the scale length of background curvature is comparable to the mean free path of gravitons, there will be further effects. We suspect that any such effect, which involves scattering of gravitons by the curvature will only go towards randomising the processes and will not affect our conclusions but this claim is hard to establish.

Finally, we would like to speculate on some implication our results might have for the issue of information loss in black hole evaporation. It is well-known that a black hole radiates a Planck spectrum for linearized gravitons also, besides the Planck spectrum for other quantum fields of lower spin. For all practical purposes the radiated particles behave like a set of particles thermalised by some physical process. On the other hand, complete evaporation of a black hole resulting in purely thermal radiation as relic will lead to well-known information loss problems and possibly problems with unitary evolution of quantum mechanical systems.

In the context of our work one could not help speculating whether the restrictions on thermalization of gravitational radiation by ”normal” systems have any implications for the radiation of thermal gravitons by black holes. If the black-hole graviton system can be treated like the ‘gravitational radiation in a box’ system discussed above, and if the system does not satisfy the conditions necessary for thermalization, the escaping non-thermal gravitational radiation could carry information.

It is not clear whether one could interpret the black hole radiation as arising out of a ”thermalization” process; but if it is, then the possibility opens up that the gravitons can contain at least part of the information lost when the event horizon has formed and may lead to some interesting implications for the information loss problem. This idea, though rather speculative at this stage, deserves to be investigated further.


text:
# Questions of Stability near Black Hole Critical Points

## I Acknowledgments

The authors graciously thank V.P. Frolov for helpful conversations. GG is also grateful for the Killam Trust for its financial support.


text:
# Exchange Driven Growth

## I Introduction

A multitude of growth phenomena in physical processes are driven by exchange of particles between clusters. Examples include droplet growth via evaporation and re-condensation meakin , island growth in deposition processes az , and phase ordering ls ; ajb ; sm . Exchange processes have been also used to model social and economical systems including segregation of heterogeneous populations ts , the distribution of wealth in a society ikr , and growth of urban populations lr ; klz .

In exchange processes, clusters are composed of ‘atoms’ (monomers). Monomers detach from one cluster and re-attach to another cluster. We shall consider the detachment controlled limit where the time scale for transport between clusters is much faster then the time scale for detachment. Exchange processes incorporate both reversible and irreversible features. Clusters may grow or shrink, yet when a monomer attaches to another cluster, its respective cluster disappears. This irreversible step provides the mechanism for cluster growth. Therefore, exchange-driven processes are fundamentally different from irreversible growth processes, particularly aggregation mvs ; sc ; skf .

Such mass transfer processes are governed by an exchange kernel $`K(i,j)`$ that represents the rate of transfer of monomers from a cluster of size $`i`$ to a cluster of size $`j`$. Generally, the rate of monomer exchange between two clusters depends on their sizes. Moreover, we consider the case where there is no preferable direction for exchanges, i.e., symmetric exchange kernels, $`K(i,j)=K(j,i)`$. This is unlike migration processes where the exchange is preferential (“big gets bigger” or “rich gets richer”). Migration underlies certain physical processes (e.g., coarsening with conserved order parameter ls ; ajb ) as well as social and economical processes ikr ; lr ; klz .

We investigate homogeneous exchange kernels, $`K(ai,aj)=a^{2\lambda }K(i,j)`$. In particular, we consider the product kernel $`K(i,j)=(ij)^\lambda `$ and its generalization $`K(i,j)=i^\nu j^\mu +i^\mu j^\nu `$ with $`\nu +\mu =2\lambda `$ and $`\nu \mu `$. We obtain a complete description of the problem in the asymptotic scaling regime. The overall range of possible behaviors and the emergence of self-similar size distributions are as in aggregation and migration processes. However, there are quantitative and qualitative differences. Unlike aggregation, the gelation transition is complete, and unlike migration, the size distributions are extended rather than compact.

The behavior falls into three categories. I) Growth: When $`\nu <2`$ and $`\lambda <3/2`$, clusters grow indefinitely. The typical cluster size grows algebraically with time, $`kt^{1/(32\lambda )}`$, and the cluster size distribution is given by a self-similar distribution with a stretched exponential tail. II) Gelation: When $`\nu <2`$ and $`\lambda >3/2`$, the entire mass in the system is transformed into an infinite gel in a finite time. The cluster mass diverges algebraically near the gelation point, $`k(t_ct)^{1/(32\lambda )}`$ and a scaling behavior similar to the one underlying the growth phase is found. In the borderline case $`\nu =2`$ the scaling function has an algebraic tail with a universal exponent $`\mathrm{\Phi }(x)x^5`$. Scaling breaks down in the special point $`\nu =\mu =2`$ where the distribution is log-normal. III) Instant gelation: When $`\nu >2`$, the gelation time vanishes logarithmically with the system size, $`t_c[\mathrm{ln}N]^{(\nu 2)}`$. In particular, for an infinite system, gelation is instantaneous.

This paper is organized as follows. In the next section, we define the exchange process. The governing equations are analyzed using scaling techniques and exact solutions for the moments. We first analyze the product kernel (Sec. III) and then, the generalized kernel (Sec. IV). The gelation time in finite systems is investigated in Sec. V using heuristic arguments and numerical simulations. Applications to coarsening in the Ising model with infinite range Kawasaki dynamics and in electrostatically driven granular layers are briefly discussed in Sec. VI, and conclusions are given in Sec. VII.

## II Exchange Processes

We consider the following elementary exchange process. The system consists of an ensemble of clusters and clusters evolve via transfer of a single monomer from one cluster to another. Symbolically,

$$(i,j)\stackrel{K(i,j)}{}(i\pm 1,j1)$$
(1)
with $`i`$ and $`j`$ the number of particles in each cluster and $`K(i,j)`$ the exchange kernel. In an exchange event, a cluster is equally likely to gain or to lose a particle. Since the exchange process is unbiased, the matrix of transition rates is symmetric: $`K(i,j)=K(j,i)`$.

Let $`A_k(t)`$ be the density of clusters containing $`k`$ monomers at time $`t`$. It evolves according to the following rate equation

$$\frac{dA_k}{dt}=\underset{i,j}{}A_iA_jK(i,j)\left[\delta _{k,i+1}+\delta _{k,i1}2\delta _{k,i}\right].$$
(2)
This equation assumes perfect mixing, or equivalently, absence of spatial correlations. We restrict our attention to monodisperse initial conditions, $`A_k(0)=\delta _{k,1}`$. The exchange process has a single conservation law. As reflected by the evolution equations, the total mass is conserved, $`M_1=1`$ with $`M_a=_nk^aA_k(t)`$ the moments of the size distribution. It is natural to consider homogeneous kernels, $`K(ai,aj)=a^{2\lambda }K(i,j)`$ with $`2\lambda `$ the homogeneity degree and we present results for the product kernel $`K(i,j)=(ij)^\lambda `$ and the generalized homogeneous kernel $`K(i,j)=i^\nu j^\mu +i^\mu j^\nu `$ with $`\mu +\nu =2\lambda `$.

## III The product kernel

For the product kernel, $`K(i,j)=(ij)^\lambda `$, the rate equations (2) read

$`{\displaystyle \frac{dA_k}{dt}}=M_\lambda \left[(k+1)^\lambda A_{k+1}+(k1)^\lambda A_{k1}2k^\lambda A_k\right]`$
with the boundary condition $`A_00`$. These evolution equations demonstrate the diffusive character of the exchange process. Absorbing the factor $`M_\lambda `$ into the time variable

$$\tau =_0^t𝑑t^{}M_\lambda (t^{})$$
(3)
we recast the governing equations into

$$\frac{dA_k}{d\tau }=(k+1)^\lambda A_{k+1}+(k1)^\lambda A_{k1}2k^\lambda A_k.$$
(4)
Alternatively, one can study integer moments of the size distribution. The total density obeys $`\frac{d}{d\tau }M_0=A_1`$, the total mass is conserved $`\frac{d}{d\tau }M_1=0`$, and higher integer moments satisfy the following hierarchy of equations

$$\frac{d}{d\tau }M_n=2\underset{l=1}{\overset{[n/2]}{}}\left(\genfrac{}{}{0pt}{}{n}{2l}\right)M_{n2l+\lambda }.$$
(5)
Only for integer values of the homogeneity index is this hierarchy closed. We employ different approaches for different $`\lambda `$’s. For $`\lambda <2`$, we perform a scaling analysis of the rate equations and for $`\lambda 2`$, we analyze the moment equations. This general analysis is augmented by exact solutions for the integer values $`\lambda =0`$, $`1`$, and $`2`$.

### III.1 Scaling ($`\lambda <2`$)

When $`\lambda <2`$, dimensional analysis of Eq. (4) shows that the typical cluster size grows as

$$k\tau ^\alpha ,\mathrm{with}\alpha =\frac{1}{2\lambda }.$$
(6)
Using $`\frac{d\tau }{dt}=M_\lambda \tau ^{\alpha (\lambda 1)}`$, the growth of the typical scale is expressed in terms of the physical time

$$k\{\begin{array}{cc}t^\beta \hfill & \lambda <3/2\text{,}\hfill \\ \mathrm{exp}(\mathrm{const}.\times t)\hfill & \lambda =3/2\text{,}\hfill \\ (t_ct)^\beta \hfill & 3/2<\lambda <2\text{.}\hfill \end{array}$$
(7)
The dynamical exponent is $`\beta =(32\lambda )^1`$. As long as $`\lambda <3/2`$, clusters grow indefinitely and the characteristic size grows algebraically with time. For $`\lambda >3/2`$, a gelation transition occurs, i.e., the system develops a giant cluster in a finite time $`t_c`$.

We seek a scaling solution of the rate equations

$$A_k(\tau )\tau ^{2\alpha }\mathrm{\Phi }\left(k\tau ^\alpha \right).$$
(8)
Mass conservation dictates the normalization $`J_1=1`$ where $`J_a=𝑑xx^a\mathrm{\Phi }(x)`$ is the $`a^{\mathrm{th}}`$ moment of the scaling distribution. Technically, the scaling function describes the behavior in the limits $`k\mathrm{}`$, $`\tau \mathrm{}`$ with the variable $`x=k\tau ^\alpha `$ fixed. Thus, we consider the continuum limit of the rate equation $`\frac{}{\tau }A(k,\tau )=\frac{^2}{k^2}\left[k^\lambda A(k,\tau )\right]`$. The scaling function satisfies the second order linear differential equation

$$(2\lambda )\frac{d^2}{dx^2}\left[x^\lambda \mathrm{\Phi }(x)\right]+x\frac{d}{dx}\mathrm{\Phi }(x)+2\mathrm{\Phi }(x)=0.$$
(9)
Multiplying this equation by $`x`$, employing the identities $`x^2\mathrm{\Phi }^{}+2x\mathrm{\Phi }=(x^2\mathrm{\Phi })^{}`$, $`x\mathrm{\Psi }^{\prime \prime }=(x\mathrm{\Psi })^{\prime \prime }2\mathrm{\Psi }^{}`$, and integrating once yields $`(2\lambda )\left[(x^{\lambda +1}\mathrm{\Phi })^{}2x^\lambda \mathrm{\Phi }\right]+x^2\mathrm{\Phi }(x)=0`$. Integrating a second time gives the scaling function:

$$\mathrm{\Phi }(x)=Cx^{1\lambda }\mathrm{exp}\left[\frac{x^{2\lambda }}{(2\lambda )^2}\right]$$
(10)
with $`C=(2\lambda )^{2/(2\lambda )}/\mathrm{\Gamma }\left(\frac{1}{2\lambda }\right)`$ found from the condition $`J_1=1`$. The nature of the scaling function differs from that found for migration where $`K(l,m)=0`$ for $`l<m`$ lr : Exchange is characterized by extended distributions, while migration is characterized by compact distributions.

There are two physically relevant cases for which the rate equations can be solved exactly. When the exchange kernel is independent of the cluster size ($`\lambda =0`$), the rate equation is $`\frac{d}{d\tau }A_k=A_{k+1}+A_{k1}2A_k`$ and the cluster size distribution is ikr

$$A_k=e^{2\tau }\left[I_{k1}(2\tau )I_{k+1}(2\tau )\right],$$
(11)
where $`I_n`$ is the modified Bessel function of order $`n`$ bo . In agreement with the general scaling analysis, the typical scale grows diffusively, $`k\tau ^{1/2}`$, and the scaling function is given by $`\mathrm{\Phi }(x)=(4\pi )^{1/2}x\mathrm{exp}(x^2/4)`$.

For the pure product kernel ($`\lambda =1`$), the rate equations read $`\frac{d}{dt}A_k=(k+1)A_{k+1}+(k1)A_{k1}2kA_k`$ (in this case $`t=\tau `$). Substituting the mass-conserving ansatz $`A_k=(1u)^2u^{k1}`$ reduces the infinite set of rate equations into a single ordinary differential equation $`\frac{d}{dt}u=(1u)^2`$ subject the initial condition $`u(0)=0`$. The size distribution in this case is

$$A_k=\frac{t^{k1}}{(1+t)^{k+1}}.$$
(12)
The typical cluster size grows ballistically, $`kt`$, and the scaling function is purely exponential, $`\mathrm{\Phi }(x)=e^x`$, again in agreement with the above scaling results.

When $`3/2<\lambda <2`$, an infinite cluster is formed at some finite time $`t_c`$, termed the gelation time. The gelation time depends on the initial condition and its determination requires the full time dependent behavior. Even without knowing the gelation time exactly, one can describe the behavior in the pregel stage since the size distribution still admits the scaling form (10). Thus, for all $`\lambda <2`$ we have

$$A_k(\tau )Ck^{1\lambda }\tau ^{\frac{3\lambda }{2\lambda }}\mathrm{exp}\left[\frac{k^{2\lambda }\tau ^1}{(2\lambda )^2}\right]$$
(13)
implying that $`A_k0`$ in the limit $`\tau \mathrm{}`$ ($`tt_c`$). In other words, the gelation is complete, $`A_k(t)=0`$, for $`tt_c`$. This behavior differs from aggregation processes where at the gelation time, all mass in the system is contained in finite size clusters vE ; van ; FL .

Complete gelation can be alternatively shown as follows. Let us assume that the cluster size distribution approaches a constant $`A_kA_k^{}>0`$ as $`\tau \mathrm{}`$. From Eqs. (4), the quantities $`B_kk^\lambda A_k^{}`$ satisfy the discrete Laplace equation $`B_{k+1}+B_{k1}2B_k=0`$ for $`k>1`$ and $`B_2=2B_1`$. Solving recursively yields $`B_k=kB_1=kA_1^{}`$ or $`A_k^{}=k^{1\lambda }A_1^{}`$. Mass conservation, $`_kkA_k^{}=1`$, implies $`A_1^{}=0`$, and thence $`A_k^{}=0`$ for all $`k`$, i.e., complete gelation.

### III.2 Multiscaling ($`\lambda =2`$)

In this special case, the moment equations (5) are closed for $`n2`$. For example

$`{\displaystyle \frac{dM_2}{d\tau }}`$ $`=`$ $`2M_2,`$
$`{\displaystyle \frac{dM_3}{d\tau }}`$ $`=`$ $`6M_3,`$ (14)
$`{\displaystyle \frac{dM_4}{d\tau }}`$ $`=`$ $`12M_4+2M_2.`$
The solutions to these equations are combinations of exponentials: $`M_2=e^{2\tau }`$, $`M_3=e^{6\tau }`$, $`M_4=\frac{6}{5}e^{12\tau }\frac{1}{5}e^{2\tau }`$, etc. The physical time $`t=\frac{1}{2}[1e^{2\tau }]`$ is found from $`t=_0^\tau 𝑑sM_2^1(s)`$ so

$`M_2`$ $`=`$ $`(12t)^1,`$
$`M_3`$ $`=`$ $`(12t)^3,`$ (15)
$`M_4`$ $`=`$ $`{\displaystyle \frac{6}{5}}(12t)^6{\displaystyle \frac{1}{5}}(12t)^1.`$
Therefore, the gelation time is $`t_c=1/2`$. Asymptotically, the first term in (5) dominates: $`\frac{d}{d\tau }M_nn(n1)M_n`$ implying $`M_n\mathrm{exp}[n(n1)\tau ]`$ for $`n>1`$. Close to the gelation time ($`tt_c`$), the moments diverge according to

$$M_n(t)(t_ct)^{n(n1)/2}.$$
(16)
In this special case, moments exhibit multiscaling asymptotic behavior, namely, the normalized moments $`M_n^{1/n}/M_1`$ diverge.

To determine the asymptotic form of the size distribution we treat $`k`$ as a continuous variable. For $`\lambda =2`$, Eq. (4) becomes $`\frac{}{\tau }A_k=\frac{^2}{k^2}[k^2A_k]`$. This equation is equidimensional in $`k`$ bo  thereby suggesting use of the variable $`\xi =\mathrm{ln}k`$ instead of $`k`$. Making the transformation from $`A_k(t)`$ to $`A(\xi ,\tau )`$ defined via $`A_kdk=A(\xi )d\xi `$, we recast above equation for $`A_k(t)`$ into the following constant coefficients diffusion-convection equation

$$\left(\frac{}{\tau }\frac{}{\xi }\right)A(\xi ,\tau )=\frac{^2}{\xi ^2}A(\xi ,\tau ).$$
(17)
With the initial conditions $`A(\xi ,0)=\delta (\xi )`$, the solution reads $`A(\xi ,\tau )=(4\pi \tau )^{1/2}\mathrm{exp}[(\xi +\tau )^2/(4\tau )]`$. The original distribution $`A_k=k^1A(\xi )`$ is log-normal

$$A_k(\tau )(4\pi \tau )^{1/2}e^{\tau /4}k^{3/2}\mathrm{exp}\left[\frac{(\mathrm{ln}k)^2}{4\tau }\right].$$
(18)
Again, the distribution vanishes at the transition point, i.e., the gelation transition is complete. Moreover, the mass distribution is algebraic, $`A_k(t)M_0(t)k^{3/2}`$ for sufficiently small masses, $`k\sqrt{\mathrm{ln}\frac{1}{12t}}`$. The total density vanishes quite slowly near the transition point

$$M_0(t)(12t)^{1/8}\left(\mathrm{ln}\frac{1}{12t}\right)^{1/2}.$$
(19)
We note that the density follows a different law than the one characterizing higher than first moments (16).

The size distribution does not follow a scaling behavior asymptotically and the log-normal distribution is responsible for the multiscaling behavior (16) of the moments. This differs from aggregation processes where the moments diverge as $`M_n(t)(t_ct)^{\alpha _n}`$ vE  with the exponent $`\alpha _n`$ linear in $`n`$.

### III.3 Instant Gelation ($`\lambda >2`$)

Gelation is now instantaneous and complete, that is $`A_k(t)=0`$ for all $`k`$ when $`t>0`$. To prove this assertion we assume the opposite and arrive at a contradiction. Our analysis follows van Dongen’s approach developed in the context of aggregation processes van .

The moments $`M_n`$ with integer $`n2`$ evolve according to (5). The first term in the summation yields a lower bound for their growth rate

$$\frac{dM_n}{d\tau }n(n1)M_{n2+\lambda }n(n1)(M_n)^{1+\mathrm{\Lambda }}$$
(20)
with $`\mathrm{\Lambda }=\frac{\lambda 2}{n1}`$. The second inequality follows from the Jensen’s inequality as shown below. Consider the auxiliary functions $`_n`$, evolving according to

$`{\displaystyle \frac{d_n}{d\tau }}=n(n1)(_n)^{1+\mathrm{\Lambda }}.`$ (21)
Solving this equation subject to the initial condition $`_n(0)=1`$ yields $`_n=[1n(\lambda 2)\tau ]^{1/\mathrm{\Lambda }}`$. Therefore, $`_n\mathrm{}`$ as $`\tau \tau _n=n^1(\lambda 2)^1`$. Since $`M_n_n`$, the moment $`M_n`$ diverges at least at $`\tau _n`$. The series of times $`\tau _n`$ set an upper bound for the gelation time $`\tau _c`$ since all moments should be finite for $`\tau <\tau _c`$. As $`\tau _n0`$ when $`n\mathrm{}`$, we conclude that $`\tau _c=0`$ and thence, the gelation time vanishes $`t_c=0`$.

The inequality $`M_{n2+\lambda }(M_n)^{1+\mathrm{\Lambda }}`$ with $`\mathrm{\Lambda }=\frac{\lambda 2}{n1}`$ is derived as follows. Let the parameters $`p_j0`$ satisfy $`_jp_j=1`$ and let $`\mathrm{\Phi }(x)`$ be a convex function. A convex function satisfies the Jensen inequality

$`{\displaystyle \underset{j=1}{\overset{\mathrm{}}{}}}p_j\mathrm{\Phi }(x_j)\mathrm{\Phi }\left({\displaystyle \underset{j=1}{\overset{\mathrm{}}{}}}p_jx_j\right).`$ (22)
First, we substitute the coefficients $`p_j=jA_j`$ (from mass conservation $`_jjA_j=1`$) and the convex function $`\mathrm{\Phi }(x)=x^{1+\mathrm{\Lambda }}`$ ($`\mathrm{\Lambda }>0`$ for $`\lambda >2`$) into the Jensen inequality. Then, choosing $`x_j=j^{n1}`$ and using $`p_jx_j=j^nA_j=M_n`$ and $`p_j\mathrm{\Phi }(x_j)=M_{n2+\lambda }`$ yields the above inequality.

## IV Generalized kernels

The rates $`K(i,j)`$ underlying exchange processes are typically homogeneous functions of $`i`$ and $`j`$ (at least for large $`i`$ and $`j`$). We restrict ourselves to such kernels. Apart from the homogeneity degree $`2\lambda `$, homogeneous kernels are characterized by an additional exponent $`\nu `$ defined through the asymptotic $`K(1,j)j^\nu `$ as $`j1`$. For $`ij`$ the exchange kernel scales as $`K(i,j)=i^{2\lambda }K(1,j/i)i^\mu j^\nu `$ with $`2\lambda =\nu +\mu `$. Therefore, we consider a specific generalization of the product kernel that exhibits these homogeneity properties

$$K(i,j)=i^\nu j^\mu +i^\mu j^\nu .$$
(23)
More precisely, the asymptotics $`K(i,j)i^\mu j^\nu `$ occurs for $`ij`$ if $`\nu \mu `$; since the kernel is symmetric, we can assume that $`\nu \mu `$ without loss of generality. We expect that the homogeneity indices govern the overall qualitative behavior (growth, gelation, instant gelation), while the precise form of the kernel controls quantitative characteristics such as the size distribution.

For this exchange kernel, the rate equations (2) become

$`{\displaystyle \frac{dA_k}{dt}}`$ $`=`$ $`M_\mu \left[(k+1)^\nu A_{k+1}+(k1)^\nu A_{k1}2k^\nu A_k\right]`$
$`+`$ $`M_\nu \left[(k+1)^\mu A_{k+1}+(k1)^\mu A_{k1}2k^\mu A_k\right].`$
The following generalization of the modified time variable

$$\tau =_0^t𝑑t^{}\sqrt{M_\nu (t^{})M_\mu (t^{})}$$
(24)
handles the two indices symmetrically. In terms of this time, the evolution equations are

$`{\displaystyle \frac{dA_k}{d\tau }}`$ $`=`$ $`R\left[(k+1)^\nu A_{k+1}+(k1)^\nu A_{k1}2k^\nu A_k\right]`$
$`+`$ $`R^1\left[(k+1)^\mu A_{k+1}+(k1)^\mu A_{k1}2k^\mu A_k\right],`$
with $`R=\sqrt{M_\mu /M_\nu }`$. Of course, the dynamics conserve mass: $`\frac{d}{d\tau }M_1=0`$. Higher integer moments evolve according to

$$\frac{dM_n}{d\tau }=2\underset{l=1}{\overset{[n/2]}{}}\left(\genfrac{}{}{0pt}{}{n}{2l}\right)\left[RM_{n2l+\nu }+R^1M_{n2l+\mu }\right].$$
(25)
When $`\nu <2`$, the scaling analysis follows closely the product kernel case. The overall growth laws (6) and (7) remain unchanged and the homogeneity degree $`\lambda `$ characterizes the scaling behavior. However, we shall see that the individual indices $`\nu `$ and $`\mu `$ play an important role since they dictate the range for which this law holds.

We seek a scaling solution of the form (8). The scaling function $`\mathrm{\Phi }(x)`$ satisfies

$$\frac{d^2}{dx^2}\left[\left(Ux^\nu +Vx^\mu \right)\mathrm{\Phi }(x)\right]+x\frac{d}{dx}\mathrm{\Phi }(x)+2\mathrm{\Phi }(x)=0$$
(26)
with the constants $`U=\alpha ^1A`$ and $`V=\alpha ^1A^1`$ determined from the ratio $`A=\sqrt{J_\mu /J_\nu }`$. The scaling function reads

$$\mathrm{\Phi }(x)=C\frac{x}{Ux^\nu +Vx^\mu }\mathrm{exp}\left[_0^x𝑑y\frac{y}{Uy^\nu +Vy^\mu }\right].$$
(27)
The scaling solution involves three parameters $`U`$, $`V`$, and $`C`$. Substituting $`U=\alpha ^1A`$ and $`V=\alpha ^1A^1`$ into the equality $`A=\sqrt{J_\mu /J_\nu }`$ yields a closed equation for the parameter $`A`$. Once $`A`$ is determined, the parameters $`U`$ and $`V`$ follow, and finally, the amplitude $`C`$ is found from the normalization $`J_1=1`$.

We illustrate this procedure for the special case $`(\nu ,\mu )=(1,0)`$, i.e., for the pure sum kernel $`K(i,j)=i+j`$. In this case, the integral on the right-hand side of Eq. (27) is readily computed. Using $`U=\frac{3A}{2}`$ and $`V=\frac{3}{2A}`$ we arrive at

$$\mathrm{\Phi }(x)=Cx\left(1+A^2x\right)^{a1}\mathrm{exp}\left[aA^2x\right]$$
(28)
with $`a=\frac{2}{3}A^3`$. We now substitute this solution into the right-hand side of the equality $`A=\sqrt{J_0/J_1}`$ and transform it into the transcendental equation $`(\frac{e}{a})^a\mathrm{\Gamma }(a,a)+a^1=1`$ involving the incomplete gamma function (see Appendix A). The amplitude is then explicitly evaluated to give $`C=aA^6`$. Numerically, we find $`a2.82649`$, $`A0.428397`$, and $`C0.0174713`$. Interestingly, there is a nontrivial algebraic correction to the leading exponential behavior, $`\mathrm{\Phi }(x)x^a\mathrm{exp}(aAx)`$ for large $`x`$.

On the boundary $`\nu =2`$ separating regime III from the two other regimes, the solution of Eq. (27) significantly simplifies. We find $`A=1/[2(2\mu )]`$, $`U=1/4`$, $`V=1/(2\mu )^2`$, and consequently, the scaling function

$$\mathrm{\Phi }(x)=Cx^{1\mu }\left[1+\frac{x^{2\mu }}{4(2\mu )^2}\right]^{1\frac{4}{2\mu }}.$$
(29)
The constant $`C=2\left[2(2\mu )\right]^{1\frac{2}{2\mu }}[B(\frac{1}{2\mu },\frac{3}{2\mu })]^1`$ is expressed in terms of the beta function. Remarkably, the scaling function (29) exhibits a universal large-$`x`$ asymptotic behavior

$$\mathrm{\Phi }(x)x^5.$$
(30)
In other words, the size distribution is algebraic, $`A_k(\tau )k^5\tau ^{3\alpha }`$ where $`\alpha =(2\lambda )^1=2/(2\mu )`$. With this algebraic divergence, sufficiently small moments are characterized by ordinary scaling behavior while higher moments exhibit multiscaling behavior:

$$M_n\{\begin{array}{cc}\tau ^{\alpha (n1)}\hfill & n<4\text{,}\hfill \\ \tau ^{\alpha n(n1)/4}\hfill & n>4\text{.}\hfill \end{array}$$
(31)
This behavior follows from the leading term in the moment equation (25), viz. $`\frac{d}{d\tau }M_n=n(n1)M_nR`$. With $`R=\sqrt{M_\mu /M_2}A\tau ^1`$ and $`A=\alpha /4`$, this equation becomes $`\frac{d}{d\tau }M_n=\frac{\alpha n(n1)}{4\tau }M_n`$, leading to the multiscaling behavior (31).

The determination of $`A`$ in the general situation requires numerical evaluation, yet the form and nature of the size distribution is clear. For example, the minimal (maximal) index governs the distribution of small (large) clusters. Indeed, from Eq. (27), the extremal behaviors are

$$\mathrm{\Phi }(x)\{\begin{array}{cc}x^{1\mu }\hfill & x1\text{,}\hfill \\ \mathrm{exp}[x^{2\nu }]\hfill & x1\text{.}\hfill \end{array}$$
(32)
Apart from the point $`(\nu ,\mu )=(2,2)`$, the scaling solution holds for all $`\nu 2`$. As in the product kernel case, growth occurs when $`\lambda <3/2`$ and gelation occurs when $`3/2<\lambda <2`$.

For $`\nu >2`$, the scaling solution (27) predicts $`\mathrm{\Phi }x^{1\nu }`$. Such behavior is inconsistent since the moment $`J_\mu `$ diverges and instead, instantaneous gelation occurs. The moments $`M_n`$ with $`n>1`$ satisfy Eq. (25) and the first term in the summation yields a lower bound for the moment growth $`\frac{d}{d\tau }M_nRn(n1)M_{n2+\nu }`$. Keeping only this term and absorbing the factor $`R`$ into the time variable, the previous proof applies. Thus, gelation is instantaneous.

To summarize, there are 3 types of behaviors, determined by the homogeneity degrees $`\mu `$ and $`\nu `$ (Fig. 1):

1. Growth. The cluster size grows indefinitely, and the size distribution obeys scaling.
2. Gelation. The cluster size diverges in a finite time and the size distribution follows a scaling solution near the gelation time. Gelation is complete.
3. Instant Gelation. The cluster size distribution vanishes for all $`t>0`$.

The cluster size distribution exhibits a scaling behavior in regimes I and II. Scaling behavior underlies the system everywhere except for regime III and the point $`(2,2)`$. In the bulk of regimes I and II the size distribution is a stretched exponential, while in the boundary with region III, the cluster size distribution has an algebraic tail. Last, at the point $`(\mu ,\nu )=(2,2)`$ scaling breaks down and the distribution is log-normal.

## V The gelation time

Instantaneous gelation is certainly counter intuitive: A finite time singularity that occurs at time $`t=+0`$ ! Instantaneous gelation was investigated exclusively in the context of aggregation FL ; dom ; hez ; spouge ; van ; jeon ; mg . For infinite systems, it is impossible to quantify the difference between two instant gelling. Finite systems, on the other hand, naturally quantify how fast a system gels.

Consider a system consisting initially of $`N`$ monomers. In a finite time $`t_N`$, all mass in the system condenses into a single ‘runaway’ cluster. How does the average gelation time $`T_N=t_N`$ depend on $`N`$? When growth or ordinary gelation occurs, the answer follows from our previous analysis. In the scaling regime, the growth law (7) indicates that the condensation time grows algebraically with the system size, $`T_NN^{1/\beta }`$. In the case of ordinary gelation, the average gelation time saturates at an $`N`$-independent value: $`T_Nt_c`$. The interesting case is instant gelation where the gelation time vanishes in the thermodynamic limit, $`T_N0`$ as $`N\mathrm{}`$.

For simplicity, we discuss the product kernel. The vanishing gelation time is ultimately related to the short time behavior. Early on, loss terms in the rate equation (4) are negligible and to leading order $`\frac{d}{dt}A_j(j+1)^\lambda A_{j+1}`$, where we tacitly assumed $`\tau t`$. For the initial condition $`A_j(0)=\delta _{j,1}`$, the leading order behavior of the density is

$$A_{j+1}(j!)^{\lambda 1}t^j.$$
(33)
In a finite system consisting initially of $`N`$ monomers, a $`j`$-mer first appears at time $`t_j(j!)^{(\lambda 1)/j}N^{1/j}`$, estimated from $`NA_j(t_j)=1`$. For example, the first dimer and trimer appear at times $`t_2=N^1`$ and $`t_3=2^{(\lambda 1)/2}N^{1/2}`$, respectively. By definition, the times increase monotonically, $`t_{j+1}>t_j`$, yet the above estimates increase monotonically only for sufficiently small $`j<j_{}`$. From $`t_j_{}=t_{j_{}+1}`$, we obtain the extremum $`j_{}=(\lambda 1)^1\mathrm{ln}N`$ using the Stirling formula. The corresponding time $`T_{}t_j_{}`$ is

$$T_{}\left(\frac{\lambda 1}{\mathrm{ln}N}\right)^{\lambda 1}.$$
(34)
For later times, $`tT_{}`$, the rate equations should be modified to account for the finiteness of the system (see e.g. aal ; tn ) since significant statistical fluctuations are induced by large runaway clusters that take over (eventually only one such cluster remains). The critical size of such clusters is $`j_{}(\lambda 1)^1\mathrm{ln}N`$. As a complete analytical solution seems out of reach, we proceed heuristically by focusing on the leading cluster that eventually grows to be the gel. Since it exchanges monomers back and forth with other clusters its growth mechanism is diffusive. For an ordinary diffusive process, $`\frac{d}{dt}k=0`$, while $`\frac{d}{dt}k^2=D`$. In our case, $`D=k^\lambda `$ with the typical size $`k^2k^2`$. Therefore, the typical size of the runaway cluster grows according to

$$\frac{dk}{dt}=k^{\lambda 1}.$$
(35)
Integrating this rate equation from the critical size $`k=j_{}`$ to the system size $`k=N`$ gives the gelation time

$$T_N=T_{}+\frac{1}{\lambda 2}\left[\frac{1}{j_{}^{\lambda 2}}\frac{1}{N^{\lambda 2}}\right].$$
(36)
Since $`j_{}=(\lambda 1)^1\mathrm{ln}N`$, the duration of the latter growth phase is much larger than that of the nucleation phase, $`T_NT_{}`$. Therefore, the gelation time vanishes logarithmically

$$T_N(\mathrm{ln}N)^{(\lambda 2)},$$
(37)
in the thermodynamic limit. A straightforward extension of the above argument to the generalized exchange kernel (23) gives $`T_N(\mathrm{ln}N)^{(\nu 2)}`$.

Therefore, in a finite system it may be difficult to distinguish instantaneous gelation from the ordinary one. We verified the logarithmic law (37) numerically for $`\lambda =3`$ (Fig. 2). Probing fluctuations in the gelation time, we observe that the normalized variance $`\sigma _N^2=t_N^2/t_N^21`$ vanishes logarithmically in the thermodynamic limit (Fig. 2). The distribution of normalized gelation times becomes trivial $`P(t_N/T_N)\delta (z1)`$ implying that the gelation time is a self-averaging quantity.

We also examined the gelation time in two other growth processes, namely aggregation mvs ; sc ; skf  and addition bk ; l . The above heuristic picture yields a similar logarithmic law albeit with a different exponent mg . Self-averaging is observed numerically as well, and we conclude that the behavior found for exchange processes is generic.

## VI Applications to Coarsening

In exchange processes, a monomer detaches from a cluster and subsequently re-attaches to another cluster. This elementary mechanism underlies a number of growth and coarsening processes. We apply our general theory to two coarsening processes.

### VI.1 Infinite range Ising-Kawasaki model

The zero temperature Ising model with infinite range Kawasaki dynamics evolves via exchanges of spins of opposite signs belonging to domain walls kk ; tk ; ab . This model is equivalent to the Shelling’s segregation model ts . In the limit of a vanishing volume fraction of one of the two phases, the domains of the minority phase are isolated and the process is essentially an exchange process with a product kernel $`K(i,j)=(ij)^\lambda `$. The dependence of the number of exchange candidates (i.e., spins in domain walls that can lower their energy by hoping to a different cluster) on the cluster size dictates the homogeneity degree. For spherical clusters, only perimeter spins may exchange. Since the island size and the surface size grow with the radius according to $`kR^d`$ and $`\sigma R^{d1}k^{(d1)/d}`$, respectively, one has $`\lambda =(d1)/d`$ and $`\beta =1/(32\lambda )=d/(d+2)`$. The dynamical exponent (defined through $`Rt^z`$) is therefore $`z=1/(d+2)`$. If the islands are polygons, a distinct possibility on a lattice, then only corner spins are active so $`\lambda =0`$, and consequently $`\beta =1/3`$ and $`z=1/(3d)`$. Both estimates agree in one dimension, consistent with the exact solution $`z=1/3`$ bk1 .

### VI.2 Coarsening of thin granular layers

In electrostatically driven granular layers, clusters nucleate around large grains ia . Then, as charged grains oscillate back and forth between the two bounding plates due to the oscillating electric field, they may scatter of the plate or collide with other particle. Consequently, individual grains may transfer from one cluster to another. The rate of hoping into and out of a cluster is proportional to its area. Therefore, the homogeneity degree is unity, $`\lambda =1`$, implying $`\beta =1`$ and a dynamical exponent of $`z=1/d`$. In two dimensions, this prediction is consistent with the experimental observations $`z=1/2`$ ia . The corresponding size distribution $`\mathrm{\Phi }(x)=e^x`$ provides a reasonable approximation for the experimental size distribution, obtained from a relatively small number of clusters (see Fig. 3). The normalized variance, defined via $`\sigma ^2=J_2/J_1^21`$, is experimentally found to be $`\sigma =0.80\pm 0.05`$ compared with the theoretical value $`\sigma =1`$. Comparing ”spatial” exchange processes where grains are exchanged only between neighboring clusters with further experimental data may help elucidate the relevance of spatial correlations.

## VII Conclusions

We have shown that kinetics of exchange processes are classified by the homogeneity indices of the governing rates. There are three possible regimes including indefinite growth, gelation in a finite time, and instant gelation. Scaling behavior underlies the first two regimes. The size distributions are generally extended, decaying exponentially or algebraically for large sizes, in contrast with migration processes.

We also studied the gelation time in finite systems and found that it decays rather slowly, following an inverse logarithmic law. It would be interesting to determine the full time dependent behavior of moments of the size distribution in the instant gelation regime. In the realm of aggregation, instantaneous gelation is relevant in astrophysics, e.g., to the coalescence of planetesimals to form planets or of stars to form super-massive black holes (see mg  and references therein). The physical relevance of instant gelation in exchange processes has yet to be demonstrated.

Our description was on a mean-field level where all pairs of clusters in the system are equally likely to interact. It will be interesting to incorporate spatial fluctuations into this description. The nature of the spatial fluctuations depends on the mechanism for transporting monomers from one cluster to the other. For diffusive transport, one can incorporate effective fluxes into clusters, using the standard techniques developed for reaction-diffusion processes.

###### Acknowledgements.

We thank Igor Aronson and Avner Peleg for useful discussions. We are grateful to Igor Aronson for providing the experimental data. This research was supported by DOE (W-7405-ENG-36).

## Appendix A The case $`(\nu ,\mu )=(1,0)`$

Substituting Eq. (28) into $`A^2=J_0/J_1`$ yields

$`1`$ $`=`$ $`{\displaystyle \frac{_0^{\mathrm{}}𝑑xx^2(1+x)^{b1}e^{ax}}{_0^{\mathrm{}}𝑑xx(1+x)^{b1}e^{ax}}}`$
with $`b=a`$. Evaluation of the ratio of the integrals is performed as follows

$`1`$ $`=`$ $`{\displaystyle \frac{d}{da}}\mathrm{ln}\left[{\displaystyle _0^{\mathrm{}}}𝑑xx(1+x)^{b1}e^{ax}\right]|_{a=b}`$
$`=`$ $`{\displaystyle \frac{d}{da}}\mathrm{ln}\left[{\displaystyle \frac{d}{da}}\left({\displaystyle _0^{\mathrm{}}}𝑑x(1+x)^{b1}e^{ax}\right)\right]|_{a=b}`$
$`=`$ $`{\displaystyle \frac{d}{da}}\mathrm{ln}\left[{\displaystyle \frac{d}{da}}\left(e^aa^b\mathrm{\Gamma }(b,a)\right)\right]|_{a=b}`$
$`=`$ $`{\displaystyle \frac{a^2+e^aa^{a1}\mathrm{\Gamma }(a,a)}{a^1}}.`$