text:
# The distance to the nearest star-forming clouds: MBM12 and MBM20.

## 1 Introduction

The two nearest known star forming clouds to the Sun are MBM12 and MBM20. The MBM12 complex consists of clouds 11, 12, and 13 from the catalog of Magnani et al. (1985) and is located at (l,b) $``$ (159$`\stackrel{}{.}`$4,$``$34$`\stackrel{}{.}`$3) (we will refer to the entire complex as MBM12). It was first identified by Lynds (1962) and appears as objects L1453-L1454, L1457, L1458 in her catalog of dark nebulae. The mass of the entire complex is estimated to be $``$ 30–200 M based on radio maps of the region in 12CO, 13CO and C18O (Pound et al. 1990; Zimmermann & Ungerechts 1990). The cloud MBM20 was first identified as L1642 in the Lynds (1962) catalog of dark nebulae. The mass of MBM20 is $``$ 84 M and it is located at (l,b) $``$ (210$`\stackrel{}{.}`$9,$``$36$`\stackrel{}{.}`$6) southwest of the Orion star forming complex (Magnani et al. 1985).

Both MBM12 and MBM20 are star forming clouds. Several T Tauri stars have been identified in MBM12 via H$`\alpha `$ emission line and X-ray surveys (Herbig & Bell 1988; Stephenson 1986; Hearty et al. 2000). Two binary classical T Tauri stars have been identified in the central region of MBM20 (Sandell et al. 1987). Subsequent high-resolution 12CO(J=1-0) observations near the center of the cloud have shown that the infrared point source IRAS04325-1419 corresponding to one of the binaries is associated with a red and blue-shifted bipolar outflow (Liljeström et al. 1989).

Using ROSAT observations, Kuntz et al. (1997) find evidence based on the possible detection of 0.25 keV X-ray shadows of these clouds that both are located either within or at the edge of the Local Bubble. They estimate that distances, $`d`$, of $`60d90`$ pc to MBM12 and $`d`$ $`<`$ 100 pc to MBM20 are consistent with the foreground 0.25 keV emission seen in the direction of each cloud. The distances are also consistent with previous distance estimates (see Sect. 2).

Since MBM12 and MBM20 are both star forming clouds which are probably located within or at the edge of the Local Bubble, determining an accurate distance to these clouds is important for studies of stars and of the interstellar medium. In addition, observations of MBM12 are already scheduled for XMM (50 ks) and Chandra (100 ks) to investigate the diffuse X-ray emission of the Local Bubble. Therefore we present observations to improve the distance estimates to both clouds. In Sect. 2 we review previous distance estimates for these clouds and the Hipparcos results that revise them. In Sect. 3 we present our spectroscopic observations of Hipparcos stars along the line of sight to both of the clouds which we use to improve the previous distance estimates. In Sect. 4 we summarize our investigations and suggest future observations.

## 2 Previous distance estimates

The distances to MBM12 and MBM20 have previously been derived using both photometric and spectroscopic techniques. The spectroscopic technique used in this paper was first used by Hobbs et al. (1986) and Hobbs et al. (1988) to estimate the distance to MBM12 and by Penprase (1993) to estimate the distance to MBM20. These studies looked for interstellar Na I D lines at $`\lambda 5889.950`$ Å and $`\lambda 5895.924`$ Å in bright stars in the direction of each cloud. Since the distance to each star was known via its spectroscopic parallax, the distance to the cloud could be determined. The stars with interstellar Na I absorption features were presumed to be behind the cloud and those without interstellar features were presumed to be in front of the cloud or not located behind a sufficient column density to detect the interstellar line. Since the Hipparcos satellite measured the trigonometric parallax for most of the stars used to determine the distance to MBM12 and MBM20, it is no longer necessary to assume a spectral type or intrinsic luminosity (as is necessary for a spectroscopic parallax) to measure their distance.

The stars used by Hobbs et al. (1986), Hobbs et al. (1988), and Penprase (1993) to establish the distance to MBM12 and MBM20 are listed in Table 1 with their apparent magnitude, spectral type, distance based on spectroscopic parallax, distance based on the Hipparcos parallax111All Hipparcos distance errors, $`\sigma _d`$, in this paper are approximated based on the parallax, $`\pi `$, and standard parallax error, $`\sigma _\pi `$, from the Hipparcos catalog using the error propagation equation given by Bevington & Robinson (1992). Thus, the error in distance is given by $`\sigma _d=\sigma _\pi /\pi ^2`$, where $`\sigma _d`$ is in parcecs, and $`\pi `$ and $`\sigma _\pi `$ are in arcseconds. and whether the spectrum presented in Hobbs et al. (1986) and Hobbs et al. (1988) showed interstellar Na I absorption lines.

### 2.1 MBM12

Previous distance estimates to the cloud MBM12 have already been discussed by Hearty et al. (2000). We just note the Hipparcos results here. In their search for interstellar Na I D lines in the spectra of stars along the line of sight to MBM12, Hobbs et al. (1986) found that the star HD18404 (distance $``$ 60 pc) showed no interstellar absorption features and is therefore presumably in front of the cloud and the star HD18519/20 (distance $``$ 70 pc) did show interstellar absorption features and is therefore behind the cloud. According to Hipparcos, the distance to HD18404 is $``$ 32$`\pm `$1 pc and the distance to HD18519/20 is $``$ 90$`\pm `$12 pc. Although the Hipparcos results indicate the distance to MBM12 is not as well constrained, it is consistent with previous estimates (i.e., 32$`\pm `$1 pc $`<d<`$ 90$`\pm `$12 pc).

Hobbs et al. (1986) note that it remains possible, although unlikely, that a Na I cloud is located in the foreground of a more distant CO cloud at virtually the same radial velocity. Therefore it is important to note that there are at least two stars located behind the cloud which show a large extinction (see Fig. 1). However, both of these stars are too faint to have been observed with Hipparcos. One of the stars, H0253+193, is an eclipsing cataclysmic variable (Koyama et al. 1991; Zuckerman et al. 1992). Zuckerman et al. (1992) estimate a distance of $``$ 200 pc and a visual extinction $`A_v`$ $``$ 11.5 mag based on its near infrared colors. In addition, Hearty et al. (2000) showed that the optical spectrum of another star, DC48 (Duerr & Craine 1982), corresponds to a G9 star. Comparing the magnitude and color measured by Duerr & Craine (1982) for this star ($`V`$ = 18.7 and V-I = 5.6 mag) to the intrinsic values of a G9 star, it could be a main sequence star with $`A_v`$ $``$ 8.9 mag at a distance of $``$ 63 pc or a giant star with $`A_v`$ $``$ 8.4 mag at a distance of $``$ 950 pc. Future observations to better constrain the distance, the luminosity class, and the spectral type of this star may prove that the CO and the Na I clouds along this line of sight are one and the same and provide a new upper limit on the distance to MBM12.

### 2.2 MBM20

The cloud MBM20 is located southwest of several cometary globules probably associated with the Orion star forming region. However, previous photometric and spectroscopic observations of MBM20 and the morphology of the cometary globules suggest that MBM20 is much closer to the Sun than the cometary globules in the region. Studies of star formation in Bok globules in the Gum Nebula have identified Bok globules with comet-like tails pointed away from the central OB association (e.g., Reipurth 1983). This type of large scale morphology which Bally et al. (1991) argue is caused by the energy injected from massive stars can also be seen in the cometary clouds surrounding the Orion OB association where many clouds have dense heads with tails pointing away from the center of the Orion OB association. Whatever the mechanisms are that generate the cometary clouds, there are no signs that MBM20 is interacting with the bright stars in Orion (i.e., there is no cometary tail pointed away from the OB association).

Magnani & de Vries (1986) used star counts to estimate an upper limit for the distance to this cloud to be $``$ 125 pc. More recently, Penprase (1992) derived a photometric distance, $`d`$, to MBM20 of $`100d120`$ pc.

Using the same technique we use in this paper, Penprase (1993) arrived at an upper limit for the distance to MBM20 ($``$ 110 pc) by observing Na I D interstellar absorption features is stars for which a spectroscopic parallax could be determined (cf. Hobbs et al. 1986). The nearest star which showed interstellar Na I absorption, HD29851, is at a distance of 110 pc based on its spectroscopic parallax. Since none of the stars observed by Penprase (1993) were superimposed on the cloud contours, a lower limit for the distance to the cloud could not be determined. However, even without a lower limit on the distance, these observations showed that MBM20 is among the nearest star forming molecular clouds to the Sun. The Hipparcos parallax for HD29851 places it at $``$ $`161\pm 21`$ pc. Therefore, the Hipparcos results indicate the distance to the cloud could be larger than previous estimates.

## 3 Spectroscopic Observations

Since there are additional Hipparcos stars at intermediate distances to the previous estimates in the direction of both MBM12 and MBM20, we obtained high-resolution spectroscopic observations near the Na I D lines of some of these stars to refine the distances to both clouds (cf. Hobbs et al 1986). The stars used by Hobbs et al (1986), Hobbs et al. (1988), and by us to determine the distance to MBM12 are displayed on an IRAS 100$`\mu `$m image of the region in Fig. 1 and the stars used by Penprase (1993) and by us to determine the distance to MBM20 are displayed on an IRAS 100$`\mu `$m image of the region around MBM20 in Fig. 2

We obtained the spectra with the FOCES Echelle Spectrograph at the 2.2-m telescope in Calar Alto, Spain in August 1998. The Tek$`\mathrm{\#}`$13c 1K$`\times `$1K chip with the FOCES Echelle Spectrograph provided a spectral resolution, $`R`$ = $`\lambda /\delta \lambda `$ $``$ 49,000, derived from measurements of the FWHM (FWHM $``$ 0.12 Å) of several well isolated emission lines of the comparison spectra near the Na I D lines. As a rule, the smallest measurable equivalent width in units of mÅ is approximately the reciprocal dispersion in units Å mm-1 (Jaschek & Jaschek 1987). Therefore, given the 24$`\mu `$m pixels of the CCD we used for the observations, the smallest measurable equivalent width for these spectra is $``$ 5 mÅ. This corresponds to a minimum detectable column density N(Na I) $``$ 6 $`\times `$ 1010 cm-2 in the optically thin case. All spectra were given an initial inspection at the telescope. If a particular integration produced fewer than $``$ 1000 cts pixel-1, at least one additional integration was performed.

The observed stars are listed in Table 2 with their coordinates, apparent magnitude, spectral type, Hipparcos distance, and whether the observed spectrum contains interstellar Na I features. We observed HIP13914, although it was already known to show interstellar Na I absorption, to ensure that we could detect the lines, if present. Although we made no attempt to correct for telluric features, our measured equivalent widths of the interstellar Na I D1 and D2 lines 129 mÅ and 167 mÅ, respectively, in HIP13914 are in good agreement with the measurements of Hobbs et al. (1986), 127 mÅ and 149 mÅ which were performed at the same spectral resolution. All of the observed spectra are displayed in Fig. 3.

None of the observed stars along the line of sight to MBM12 except HIP13914 showed interstellar Na I absorption, therefore the upper limit on the distance to the cloud remains $``$ $`90\pm 12`$ pc. Since HIP13631 is the only star we observed that is superimposed directly on the cloud and it does not show interstellar Na I absorption, we can place a lower limit on the distance to the cloud of $``$ $`58\pm 5`$ pc. Therefore the distance to MBM12 is $`58\pm 5`$ pc $`<`$ $`d`$ $`<`$ $`90\pm 12`$ pc.

The observed stars along the line of sight to MBM20 are all projected on the central region of that cloud. Since we did not detect interstellar Na I D lines in any of the observed stars, we can place a lower limit on the distance to MBM20 of $``$ $`112\pm 15`$ pc. Therefore the distance to MBM20 is $`112\pm 15`$ pc $`<`$ $`d`$ $`<`$ $`161\pm 21`$ pc.

## 4 Conclusions

We have investigated the Hipparcos parallax measurements of stars to estimate the distance to the two nearest star forming clouds to the Sun: MBM12 and MBM20. The Hipparcos observations of previously observed stars suggest that the distances to both of these clouds are not as well constrained as previously thought, although they are consistent with the previous estimates. Therefore, we obtained high resolution spectra of additional Hipparcos stars at intermediate distances and improved the distance estimates to both clouds. The distance to MBM12 is 58$`\pm `$5 pc $`<d<`$ 90$`\pm `$12 pc and the distance to MBM20 is $``$ $`112\pm 15`$ pc $`<d<`$ $`161\pm 21`$ pc. Since there are additional stars that were observed with Hipparcos that are projected near each cloud (though farther off axis than the stars we observed), future spectroscopic observations like those presented here may further refine our distance estimates.

###### Acknowledgements.

We thank C. Bertout and L.M. Hobbs for comments which helped us to improve this paper. TH wishes to acknowledge a stipendium from the Max-Planck-Gesellschaft and the hospitality of Professor J. Trümper for support of this research.


format: text

text:
# Accepted for publication in Monthly Notices Linear Redshift Distortions and Power in the PSCz Survey

## 1 Introduction

The IRAS Point Source Catalog Redshift Survey (PSCz), made public on 9 January 2000 (Saunders et al. 2000), covers more volume of the Universe than any other publically available redshift survey. This, together with the considerable care taken by its authors to ensure uniformity of selection over a prescribed volume, makes the PSCz the finest available galaxy survey for measuring the power spectrum and its redshift distortions at large, linear scales.

The purpose of this Letter is to report a state-of-the-art linear redshift distortion analysis of the PSCz. Our approach has its roots in the work of Fisher, Scharf & Lahav (1994), who were the first to apply a maximum likelihood approach to the analysis of linear redshift distortions, and in the work of Heavens & Taylor (1995, hereafter HT), who may be credited with accomplishing the first linear likelihood analysis designed to retain as much information as possible at the largest, linear scales. Fisher et al. and HT applied their approaches to the the IRAS 1.2 Jy survey (Fisher et al. 1995), obtaining values of the redshift distortion parameter $`\beta =0.96_{0.18}^{+0.20}`$ and $`\beta =1.0\pm 0.3`$ respectively, where $`\beta f/b`$ is the dimensionless linear growth rate $`f\mathrm{\Omega }_\mathrm{m}^{0.6}`$ divided by the linear galaxy-to-mass bias factor $`b`$. A comprehensive review of these and other measurements of linear redshift distortions, complete up to mid-1997, may be found in Hamilton (1998).

Most recently, Tadros et al. (1999) applied HT’s method to the PSCz survey, conservatively cut to 0.75 Jy. To deal with the more complicated PSCz mask, Tadros et al. used an improved treatment of the angular mask, allowing them to include monopole and dipole modes. Tadros et al. obtained $`\beta =0.58\pm 0.26`$ (marginalized error) if the shape of the power spectrum was fixed to be that of a $`\mathrm{\Gamma }=0.2`$ CDM model, or $`\beta =0.47\pm 0.16`$ (conditional error, with the power spectrum fixed at its best fit) if the power spectrum was simultaneously measured in 6 bins.

The goal of the present paper is to measure the linear power spectrum and its redshift distortions at high resolution with minimal loss of information. The paper brings together several ideas from the literature, chief amongst which are: (1) use of the logarithmic spherical wave basis (Hamilton & Culhane 1996); (2) linear compression into Karhunen-Loève (KL), or signal-to-noise, eigenmodes (Bond 1995; Bunn & Sugiyama 1995; Vogeley & Szalay 1996; Matsubara, Szalay & Landy 2000); (3) quadratic compression into high resolution powers (Tegmark 1997, 1998, hereafter T97, T98; Tegmark, Taylor & Heavens 1997, hereafter TTH; Tegmark et al. 1998; Padmanabhan, Tegmark & Hamilton 2000); (4) inclusion of stochastic bias in the analysis (Pen 1998; Tegmark & Peebles 1998; Dekel & Lahav 1999); (5) decorrelation (Hamilton 1997; Tegmark & Hamilton 1998; Hamilton & Tegmark 2000); (6) the elimination of bias resulting from the ‘pair-integral constraint’, by the trick of making all KL modes orthogonal to the mean mode (Fisher et al. 1993; Tegmark et al. 1998).

## 2 The 2-point correlation function in redshift space

Figure 1 shows the redshift space 2-point correlation function of PSCz. The expected large scale squashing effect caused by infall toward overdense regions is plainly visible. Nonlinear fingers-of-god show up as $`4h^1\mathrm{Mpc}`$ wide by $`20h^1\mathrm{Mpc}`$ long enhancements along the line-of-sight axis.

There is little sign of the large scale ‘arm-of-god’ affliction beyond $`80h^1\mathrm{Mpc}`$ found by Hamilton (1995) in the QDOT survey, the 1-in-6 precursor to PSCz (Lawrence et al. 1999). The contour plots in Figure 1 are constructed in a manner essentially identical to those in Hamilton (1995) and Hamilton (1998). The only difference is that, thanks to the lower noise level of the PSCz, the contour plots in Figure 1 extend to larger scales and to lower contour levels than those of Hamilton (1995, 1998).

## 3 Analysis

The analysis will be described in full in a subsequent paper. Here we summarize the main points.

As emphasized by Tadros et al. (1999), in studies of weak clustering at large scales, uniformity of selection is paramount, since non-uniformity will masquerade as spurious excess power. Following the advice of Saunders et al. (2000) and Tadros et al., we adopt the high-latitude angular mask (hibpsczmask.dat, part of the PSCz package).

It is important to measure the radial selection function of the survey as accurately as possible, since errors in the selection function translate into spurious large scale power. We use the maximum likelihood method of Sandage, Tammann & Yahil (1979), and fit the selection function to a smooth analytic function with enough free parameters to yield an excellent fit. The assumed geometry is that of a flat $`\mathrm{\Omega }_\mathrm{m}=0.3`$, $`\mathrm{\Omega }_\mathrm{\Lambda }=0.7`$ Universe. The measurement yields evidence for what appears to be strong evolution, which we model as luminosity evolution with luminosity $`(1+z)^{3.4}`$. Similarly large evolution was previously reported in the QDOT survey by Saunders et al. (1990). Figure 2 shows the observed number density of galaxies, divided by the measured selection function, both with and without evolution.

Figure 2 also shows what appears to be incompleteness beyond $`420h^1\mathrm{Mpc}`$. This may be presumed to be the incompleteness at high redshift described in §4.2 of Saunders et al. (2000), associated with the policy not to pursue redshifts of galaxies optically fainter than $`b_J=19.5^\mathrm{m}`$. Since this incompleteness is greater in regions of higher optical extinction, and is systematic rather than random over the sky (Fig 4. of Saunders et al. 2000), we choose to cut the survey at $`10^{2.625}h^1\mathrm{Mpc}420h^1\mathrm{Mpc}`$. We also choose to cut the survey at a near distance of $`1h^1\mathrm{Mpc}`$, in order to eliminate the local nonlinear region, the Local Group of galaxies. This leaves 12504 galaxies in the sample.

As a working basis from which Karhunen-Loève modes (see below) are constructed, we use logarithmic spherical waves (Hamilton & Culhane 1996), which are orthonormal eigenfunctions $`Z_{\omega \mathrm{}m}(𝒓)=(2\pi )^{1/2}e^{(3/2+\mathrm{i}\omega )\mathrm{ln}r}Y_\mathrm{}m(\widehat{𝒓})`$ of the complete set of commuting Hermitian operators

$$\mathrm{i}\left(\frac{}{\mathrm{ln}r}+\frac{3}{2}\right)=\mathrm{i}\left(\frac{}{\mathrm{ln}k}+\frac{3}{2}\right),L^2,L_z.$$
(1)
The advantages of this basis are: (1) the linear redshift distortion operator is diagonal in this basis (Hamilton & Culhane 1996); (2) radial modes can be transformed rapidly between real, $`\omega `$, and Fourier space using Fast Fourier Transforms; (3) the radial modes discretize naturally on to a grid of radial depths $`r`$ that is uniformly spaced in the log, which makes it well suited to surveys like PSCz, where the flux limit causes the survey to be finely sampled nearby, and coarsely sampled far away.

The principal numerical limitation to the linear likelihood method is the number $`N`$ of modes that can be included in the Gaussian likelihood function. Solving the likelihood formulae involves manipulating $`N\times N`$ matrices, which is a $`N^3`$ process. Thus the modes should be crafted so as to include as much information as possible about the quantities of interest, the parameters to be measured (HT; TTH).

Here we use Karhunen-Loève (KL), or signal-to-noise, modes, as first suggested in the context of galaxy clustering by Vogeley and Szalay (1996). The covariance matrix of galaxy densities is a sum $`𝐂=𝐒+𝐍`$ of cosmic signal $`𝐒`$ and Poisson sampling noise $`𝐍`$. KL modes are constructed by diagonalizing the signal-to-noise matrix $`𝐍^{1/2}\mathrm{𝐒𝐍}^{1/2}`$, and selecting the eigenmodes with the largest eigenvalues.

The main difficulty with the KL procedure is that in order to construct $`N`$ modes it is necessary to diagonalize a matrix of dimension $`N`$. But this would seem to defeat the original goal, to avoid manipulating huge matrices. We get around this difficulty by constructing KL modes in a two stage process, constructing first angular KL modes, then radial KL modes within each angular mode. The resulting KL modes are not perfect, since they ignore covariances between radial modes belonging to different angular modes. But the KL modes need not be perfect; it is enough that they should contain virtually all the information of interest.

In practice, we use 1600 angular modes, and 192 radial modes, so there is a potential pool of $`1600\times 1923\times 10^5`$ modes. This is $`754^3`$ times the number 4096 of KL modes retained in the analysis. In effect, the KL modes are spatially over-resolved by a factor of about 4 in each dimension.

In order to isolate the effects of the selection function, and of the motion of the Local Group (LG) through the Cosmic Microwave Background (CMB), we choose the first four angular ‘KL’ modes to be the cut (i.e. masked) monopole, and three cut dipole modes. The cut dipole modes contain admixtures of the cut monopole, to make them orthogonal to the latter, and all other angular KL modes are orthogonal to the cut monopole and dipole. Within the cut monopole mode, the first two radial modes are the mean (selection function) mode and the Local Group mode (equation 4.42 of Hamilton 1998), the latter containing a small admixture of the mean to make it orthogonal to the cut (i.e. radially masked) mean. Within the cut dipole modes, the first radial mode is the Local Group mode.

The mean mode is used in determining the maximum likelihood normalization of the selection function, but is then discarded from the analysis, since it is impossible to measure the fluctuation of the mean mode.

Since all KL modes are orthogonal to the mean, there is no need to subtract the mean. Making the KL modes orthogonal to the mean immunizes their amplitudes against uncertainty in the mean density, and ensures that the measurement of power is unbiased by the so-called ‘pair-integral constraint’. This trick was first used by Fisher et al. (1993).

The fact that the selection function is measured in redshift space, not real space, leads to a correction to the linear redshift distortion operator for the monopole modes, described in §4.4 of Hamilton (1998). The correction leads to a noticeable reduction in measured power at the largest scales.

The motion of the LG through the CMB is known (Courteau & van den Bergh 1999; Lineweaver et al. 1996). We correct the amplitudes of the four LG modes (one cut monopole and three cut dipole) for this motion, and include them in the analysis, as did Tadros et al. (1999).

The amplitudes of the 4096 KL modes constitute the ‘data’ which now feed into the Gaussian likelihood function. Rather than attempting to solve the likelihood function directly for values of cosmological parameters, we instead quadratically compress (T97; TTH) the data into an intermediate number (here 147, plus one for the normalization of the selection function) of parameters that contain essentially all the information of interest. The advantage of quadratic compression is that it allows large numbers of parameters to be measured rapidly, without a time-consuming nonlinear search for an extremum in a high-dimensional space. Key to the quadratic method is that the intermediate parameters should be linearly related to the mean and covariance. In a linear redshift distortion analysis, the intermediate parameters that emerge naturally are the galaxy-galaxy, galaxy-velocity, and velocity-velocity power spectra (Kolatt & Dekel 1997; Tegmark 1998; Pen 1998; Tegmark & Peebles 1998; Dekel & Lahav 1999):

$$\begin{array}{ccccc}\hfill \text{galaxy-galaxy power}& :& P_{\mathrm{gg}}(k)& =& b(k)^2P(k)\hfill \\ \hfill \text{galaxy-velocity power}& :& P_{\mathrm{gv}}(k)& =& r(k)b(k)fP(k)\hfill \\ \hfill \text{velocity-velocity power}& :& P_{\mathrm{vv}}(k)& =& f^2P(k)\hfill \end{array}$$
(2)
where $`b(k)`$ is the (possibly scale-dependent) galaxy-to-mass bias factor, $`r(k)[1,1]`$ is a (possibly scale-dependent) galaxy-velocity correlation coefficient, $`f\mathrm{\Omega }_\mathrm{m}^{0.6}`$ is the dimensionless linear growth rate, and $`P(k)`$ is the mass power spectrum. More correctly, the ‘velocity’ here refers to minus the velocity divergence, which in linear theory is related to the mass (not galaxy) overdensity $`\delta `$ by $`f\delta +𝒗=0`$, where $``$ denotes the comoving gradient in velocity units.

Although the three power spectra are interesting in their own right, here we report results in terms of conventional parameters: the linear redshift distortion parameter $`\beta f/b`$, the galaxy-velocity correlation coefficient $`r`$, and the galaxy-galaxy power spectrum $`P_{\mathrm{gg}}(k)`$.

The estimate of power that emerges most directly from quadratic compression is the power spectrum convolved with the Fisher matrix (T97; TTH). Physically, this represents the power spectrum convolved with the Fourier transform of the optimally weighted survey window. We call this estimate, appropriately renormalized, the ‘correlated’ power, since the errors are correlated. Deconvolving the correlated power yields the ‘raw’ power spectrum. The raw power tends to be noisy and anti-correlated, and is not very useful. Half way between the correlated and raw power spectra is the ‘decorrelated’ power spectrum, which is partially deconvolved in such a way that estimates of power at different wavenumbers are uncorrelated with each other (Hamilton 1997; Hamilton & Tegmark 2000).

The quadratic method requires that prior powers $`P_{\mathrm{gg}}(k)`$, $`P_{\mathrm{gv}}(k)`$, and $`P_{\mathrm{vv}}(k)`$ be specified. The maximum likelihood (ML) solution is that power spectrum for which the estimated power equals the prior power. The ML solution can be obtained by folding the estimate back into the prior and iterating to convergence (Bond, Jaffe & Knox 1998).

However, in this paper the power is permitted to be an arbitrary function of wavenumber, with no penalty against violently varying power spectra. To express our Bayesian prejudice in favour of a smoothly varying power spectrum, we choose to fold back into the prior not the raw power spectrum, but rather the correlated power spectrum, which is smooth(ish). Moreover we iterate only once, since there is no point in attempting to overfit the bumps and wiggles in the power. As the initial prior power, we choose a COBE-normalized flat $`\mathrm{\Lambda }`$CDM model with $`\mathrm{\Omega }_\mathrm{m}=0.3`$, $`\mathrm{\Omega }_\mathrm{\Lambda }=0.7`$ from Eisenstein & Hu (1998).

Further, we set the redshift distortion parameter in the prior power to a constant, $`\beta =0.5`$, reflecting our Bayesian prejudice that this parameter should be constant in the linear regime. The data do not constrain the galaxy-velocity correlation coefficient $`r`$ tightly, but are consistent with $`r1`$ at linear scales, so we fix $`r=1`$ in the prior.

## 4 Results

Figure 3 shows the galaxy-galaxy power spectrum measured from the PSCz 0.6 Jy survey. This is the real space power spectrum, the redshift distortions having been isolated into the galaxy-velocity and velocity-velocity power spectra. At nonlinear scales, $`k0.3h\mathrm{Mpc}^1`$, the interpretation as a real space power spectrum becomes suspect, since the extent to which nonlinear redshift distortions can be approximated as ($`k`$-dependent) linear combinations of the linear galaxy-velocity and velocity-velocity powers is unknown.

Figure 3 also shows a linear model power spectrum from Eisenstein & Hu (1998), an observationally concordant COBE-normalized, untilted, flat $`\mathrm{\Lambda }`$CDM model with $`\mathrm{\Omega }_\mathrm{m}=0.3`$, $`\mathrm{\Omega }_\mathrm{\Lambda }=0.7`$, baryonic content $`\mathrm{\Omega }_bh^2=0.02`$, and Hubble constant $`h=0.65`$. Interestingly, the $`\mathrm{\Lambda }`$CDM power spectrum fits the PSCz data well at linear scales, $`k0.3h\mathrm{Mpc}^1`$, with no bias, $`b=1`$. The enhancement of the observed power over the model power at smaller scales, $`k0.5h\mathrm{Mpc}^1`$, can be attributed to nonlinearity (Peacock & Dodds 1996).

The real space power spectrum measured here appears consistent with the redshift space power spectrum measured from PSCz by Sutherland et al. (1999).

Figure 4 shows the redshift distortion parameter $`\beta f/b`$ measured as a function of wavenumber $`k`$. The left panel shows the values of $`\beta `$ measured separately from the ratio $`P_{\mathrm{gv}}(k)/P_{\mathrm{gg}}(k)`$ of the galaxy-velocity to galaxy-galaxy power on the assumption of unit galaxy-velocity correlation coefficient, $`r=1`$, and from the ratio $`\text{sign}(P_{\mathrm{vv}}/P_{\mathrm{gg}})\left|P_{\mathrm{vv}}/P_{\mathrm{gg}}\right|^{1/2}`$ of the square root of velocity-velocity to galaxy-galaxy power, equations (2). These are ratios of correlated power spectra, since the decorrelated galaxy-velocity and velocity-velocity powers are confusingly noisy.

The central panel of Figure 4 shows the best fit value of $`\beta `$ that comes from combining the two ratios of galaxy-velocity and velocity-velocity to galaxy-galaxy power at each wavenumber $`k`$. Also shown for comparison is the value of $`\beta `$ inferred from the ratio of smoothed quadrupole power to smoothed monopole power, measured in the fashion of Hamilton (1995, 1998). The two measures of $`\beta `$ agree at scales $`k0.04h\mathrm{Mpc}^1`$. The difference at $`k0.04h\mathrm{Mpc}^1`$ may plausibly be attributed to the inability of the quadupole-to-monopole ratio to probe the largest scales.

The right panel of Figure 4 shows the best fit value of $`\beta `$ that emerges from a least squares fit to the decorrelated powers including data cumulatively up to the wavenumber specified on the $`x`$-axis, again on the assumption of unit galaxy-velocity correlation coefficient, $`r=1`$. The shaded region bounds the 1-$`\sigma `$ uncertainty in $`\beta `$, where $`\chi ^2`$ exceeds its minimum value by 1.

In carrying out the fit we fix the shape of the galaxy-galaxy power spectrum to be that of the measured correlated power spectrum shown in Figure 3, but we allow the overall amplitude of the power spectrum to vary. The 1-$`\sigma `$ uncertainty in Figure 4 is the uncertainty marginalized over the amplitude of the power spectrum. We also tried allowing the shape of the power spectrum to vary, but the shape is already measured about as well as can be, and we could not thereby reduce $`\chi ^2`$ in a statistically significant way (reduction by at least 1 for each additional degree of freedom).

The fit is consistent with constant $`\beta `$ for $`k0.3h\mathrm{Mpc}^1`$, and as an overall best fit we quote the value $`\beta =0.41_{0.12}^{+0.13}`$ at $`k=0.274h\mathrm{Mpc}^1`$. $`\mathrm{\Lambda }`$CDM with $`\mathrm{\Omega }_\mathrm{m}=0.3`$ and $`\mathrm{\Omega }_\mathrm{\Lambda }=0.7`$ predicts $`f=0.513`$, consistent with our best fit result if IRAS galaxies are unbiased or mildly biased.

We also tried allowing the galaxy-velocity correlation coefficient $`r`$ to vary. At linear scales, $`k0.15h\mathrm{Mpc}^1`$, the best fit value is of order unity, $`r1`$, as suggested by the approximate agreement in the two curves in the left panel of Figure 4, but the uncertainty is large. At smaller scales, $`k0.15h\mathrm{Mpc}^1`$, the measured coefficient decreases below 1. This is almost certainly caused by nonlinearity, not by stochasticity in the linear galaxy-velocity correlation. In fact fingers-of-god, which appear extended along the line-of-sight in redshift space rather than compressed, are expected to cause the galaxy-velocity power to go negative at nonlinear scales, yielding a negative galaxy-velocity correlation coefficient $`r`$. The left panel of Figure 4 shows that the galaxy-velocity power indeed goes negative at $`k0.2h\mathrm{Mpc}^1`$.

## 5 Conclusions

This Letter reports a measurement of the linear power spectrum and its redshift distortions in the IRAS PSCz 0.6 Jy survey (Saunders et al. 2000) using state-of-the-art methods, including logarithmic wave expansion, Karhunen-Loève compression, quadratic compression, decorrelation, and isolation of the mean mode in order to eliminate large scale bias from the ‘pair-integral-constraint’.

The measurement yields three separate linear power spectra, the galaxy-galaxy, galaxy-velocity, and velocity-velocity power spectra.

The best fit overall value of the linear redshift distortion parameter, $`\beta f/b=0.41_{0.12}^{+0.13}`$, is consistent with, albeit slightly lower than, that reported by Tadros et al. (1999) for a more conservative 0.75 Jy cut of the survey.

The galaxy-velocity correlation coefficient $`r`$ is not tightly constrained by the data, but is consistent with being about unity at linear scales, $`r1`$.

The inferred real space power spectrum and its redshift distortions are consistent with a COBE-normalized, untilted, flat $`\mathrm{\Lambda }`$CDM model with $`\mathrm{\Omega }_\mathrm{m}=0.3`$ and $`\mathrm{\Omega }_\mathrm{\Lambda }=0.7`$, if IRAS galaxies are an unbiased tracer of mass at large scales in the Universe.

## Acknowledgements

This work was supported by NASA ATP grant NAG5-7128 and by NASA LTSA grant NAG5-6034. AJSH thanks Michael Culhane for many conversations.


text:
# Modern nuclear force predictions for the 𝛼-particle

## Abstract

We present new calculations of the $`\alpha `$-particle which are based on the most modern nucleon-nucleon interactions alone and combined with the Tucson-Melbourne or the Urbana IX three-nucleon interaction. Results for the binding energies and some properties of the wave function are given. On that phenomenological level little room is left for the action of a possible four-nucleon force.

Few nucleon bound states have received increasing attention in recent years. The possibility of solving increasingly complex systems allows to probe the underlying dynamics directly (for a detailed review, see ). Several methods have been developed and applied to the 4N system and realistic forces, the GFMC , the CHH , the SV , the CRCGV  and the FY  method. Very recently a new no core shell model calculation appeared . The 4N system is an important test ground for both, the NN and the 3N nuclear interactions, because of its strong binding. In this article we want to address the question whether at least the ground state energy of the $`\alpha `$-particle can be described by the most modern nuclear Hamiltonians.

In recent years NN forces have been tuned very well to the rich NN data base, which led to a new generation of so called realistic NN forces: Nijm I, II, AV18 and CD-Bonn. While Nijm II and AV18 are purely local, Nijm I has a weak nonlocality in form of a $`^2`$-dependence and CD-Bonn is quite nonlocal keeping the underlying Dirac structure of the one-boson-exchange without p/m expansions. All potentials describe the NN data base with a $`\chi ^2`$/data very close to 1. But they also have a large phenomenological character with a typical number of 40 fit parameters. The AV18 and CD-Bonn forces distinguish nn, pp and np interactions and thus include charge symmetry(CSB) and charge-independence(CIB) breaking. The interaction AV18 has in addition a whole set of electromagnetic corrections built in.

It is well known that the two 3N bound states, 3H and 3He, are theoretically underbound using only NN forces . We compare our latest theoretical results based on fully converged Faddeev calculations to the experimental binding energies in Table I. The two nonlocal potentials lead to less underbinding than the local ones. The 3N calculations take the full CSB and CIB of the NN force into account including the isospin $`T=3/2`$ admixtures. In case of 3He and 4He the Coulomb interaction is included and, in addition, we use also the electromagnetic corrections in case of AV18. Because of an implementation error the results for Nijm I, II changed with respect to . Table I includes also the strongly model dependent kinetic energies which are correlated to the NN correlations. The local potentials have a somewhat harder core  leading to a higher kinetic energy.

There are two additional dynamical ingredients, which should cure that underbinding, relativistic effects and three-nucleon forces (3NF). Though one can in principle always find an unitary (but very complex) transformation to a Hamiltonian without any 3NF , we keep the given 2N interactions for practical reasons and because of their generally accepted physical origin. In relation to those forces one or both of the two mentioned dynamical ingredients are needed. There are still controversies about the role of relativistic effects  and as we think also open conceptual questions. Thus we consider here only a strictly nonrelativistic framework, having however in mind that relativistic corrections should finally be added.

The topic of 3NF’s is as old as nuclear physics  and based on meson exchanges various processes have been proposed in the past (for a review see ). Among them the Fujita-Miyazawa force  with an intermediate $`\mathrm{\Delta }`$ generated by the exchange of two pions sticks out and is implemented in all modern 3NF models. Here we mention the rather popular $`2\pi `$-exchange Tucson-Melbourne model (TM) , the Brazilian version thereof  and the Urbana 3NF . There are also extensions to $`\pi `$-$`\rho `$ and $`\rho `$-$`\rho `$ exchanges . The 2$`\pi `$-exchange model has been critically reviewed recently and a modified version, TM’, has been proposed in  which satisfies at least chiral symmetry. In this article we use the TM, TM’ and the Urbana IX 3NF’s.

In the TM force enters the strong $`\pi `$NN vertex function parameterized in the form of a monopole form factor with a cut-off parameter $`\mathrm{\Lambda }`$. Choosing it around a generally accepted value one achieves the right order of the lacking binding energy in the 3N system . Because there is a strong correlation between the $`\alpha `$ binding energy, some 3N scattering observables  and the 3N binding energy, we fine tune $`\mathrm{\Lambda }`$ to the 3N binding energies of 3H and 3He and this separately for each of the four NN forces. In this step we have not yet included $`T=3/2`$ admixtures. We list the fit results in Table II. Thereby the original TM parameters for the constants a,b,c and d have been used . In case of the Nijmegen interactions we did not adjust $`\mathrm{\Lambda }`$ to the triton, because they do not include a modern specification of the 1S0 nn force. Since for 3H the lack of binding energy is a bit larger than for 3He the $`\mathrm{\Lambda }`$ values for 3H are slightly larger than for 3He.

Now having those 3N Hamiltonians at our disposal we can study the $`\alpha `$-particle. For results referring to older forces see for instance . We rewrite the Schrödinger equation into the Yakubovsky equations (YE)  and thus decompose the wave function $`\mathrm{\Psi }`$ into 18 Yakubovsky components (YC). Due to the identity of the nucleons the YC’s are not independent from each other and we can reduce their number to two: $`\psi _1`$ and $`\psi _2`$. Then the wave function reads

$$\mathrm{\Psi }=(1(1+P)P_{34})(1+P)\psi _1+(1+P)(1+\stackrel{~}{P})\psi _2$$
(1)
and is expressed with the help of the permutations $`P=P_{12}P_{23}+P_{13}P_{23}`$ and $`\stackrel{~}{P}=P_{13}P_{24}`$ where $`P_{ij}`$ are transpositions of particles $`i`$ and $`j`$. Thus $`P`$ acts on the 3-body subcluster (123) and $`\stackrel{~}{P}`$ interchanges the two two-body subclusters (12) and (34). This decomposition is highly advisable for scattering states since the boundary conditions can most easily be expressed for YC’s. In our case of a bound state it would in principle be possible to solve directly the Schrödinger equation, but the usage of two YC’s introduces in a natural manner two kinds of Jacobi coordinates which accelerates the convergence of a partial wave decomposition.

For 4 identical particles the YE’s reduce to two coupled integral equations

$`\psi _1`$ $`=`$ $`G_0t_{12}P\left[(1P_{34})\psi _1+\psi _2\right]`$ (2)
$`\psi _2`$ $`=`$ $`G_0t_{12}\stackrel{~}{P}\left[(1P_{34})\psi _1+\psi _2\right]`$ (3)
Here in addition to the permutations the free 4N propagator $`G_0`$ and the NN t-operator $`t_{12}`$ occur. $`t_{12}`$ is driven by the NN force $`V_{12}`$ through the Lippmann-Schwinger equation $`t_{12}=V_{12}+V_{12}G_0t_{12}`$. In case of 4N scattering one has to go one step further and solve the 3-body and 2+2 subcluster problems beforehand in order to define the correct cut-structure . This is not necessary for bound states and the form of Eqs. (2)-(3) is easier to handle numerically. One, presumably the most effective manner, to include 3NF’s has been given in . We stick to that. Then only the first of the two YE’s is changed into

$$\psi _1=G_0t_{12}P\left[(1P_{34})\psi _1+\psi _2\right]+(1+G_0t_{12})G_0V_{123}^{(3)}\mathrm{\Psi }$$
(4)
Here $`V_{123}^{(3)}`$ is that part of the 3NF which is symmetrical under exchange of particles 1 and 2. In case of the 2$`\pi `$-exchange TM or the Urbana 3NF such a separation into 3 parts is very natural. We solved the sets Eqs. (2)-(3) or Eqs. (4)-(3) in momentum space and in a partial wave representation. The first YC $`\psi _1`$ stands for the “3+1” partition and is naturally described by two Jacobi momenta for the 3-body subcluster and one for the relative motion of the 4th particle to the other 3. The second YC $`\psi _2`$ stands for the “2+2” partition and is naturally described by two relative momenta for the inner motion of the subclusters (12) and (34) and by the relative motion of the two subclusters. There are a lot of orbital and spin angular momenta as well as isospin quantum numbers to be coupled to $`J^\pi =0^+`$ and $`T=0`$. Both basis sets for the “3+1” and “2+2” partitions comprise about 1800 different combinations thereof in order to reach a converged description. The various momenta are discretized with roughly 35-45 grid points each. This leads to a huge absolutely full kernel matrix of dimensions $`10^8\times 10^8`$. We solve the eigenvalue problem by a Lanczos type algorithm  and make intensive use of a massively parallel supercomputer.

We would like to remark that the introduction of the two kinds of Jacobi momenta in Eqs. (2),(3) and (4) leads to additional coordinate transformations which are hidden in the operator form of the YE’s. They are equivalent to permutations. The crucial point in our calculation is the treatment of these permutations and coordinate transformations. The direct interchange of arbitrary particles is unfeasible because of the huge dimension of the problem. It is therefore necessary to interchange particles in two steps in such a way that at least one of the Jacobi momenta is not changed. This guarantees a block diagonal structure for the permutations. Only in this manner the calculation becomes feasible. For a detailed description see .

Since we allow for CIB and CSB in the NN forces in principle the dominant total isospin state $`T=0`$ has to be supplemented by $`T=1`$ and $`T=2`$ admixtures. Our estimations lead to the result that their admixtures will change the binding energy only very slightly ( $`<`$ 10 keV) and thus at this stage we neglected them. But CIB and CSB lead to the prescriptions that the NN t-operators occur in the form $`t=\frac{1}{3}t_{np}+\frac{1}{3}t_{pp}+\frac{1}{3}t_{nn}`$ in the NN isospin 1 channels. This is different from the 3N system, where the corresponding linear combination is $`t=\frac{1}{3}t_{np}+\frac{2}{3}t_{pp}`$. The $`pp`$ $`t`$-matrix also includes the effect of the Coulomb interaction. Since the bound nucleons are confined to a small space region, we can put to zero the Coulomb interaction outside a radius of $`1020`$ fm. Then the Fourier transformation of this interaction is nonsingular. The results are cut-off independent and numerically stable.

We show in Table I the $`\alpha `$-particle binding energies using NN forces only. As expected this theory underbinds the $`\alpha `$-particle. The CD-Bonn result compares well with . A correlation between the 3H and 4He binding energies found in  for simple forces remains valid also for the most modern ones. This is depicted in Fig. 1. The experimental value is close to that straight line correlation, which nourishes the hope that by curing $`B_t`$ one possibly will also cure $`B_\alpha `$.

In Table II we present our results adding the TM 3NF adjusted individually to the 3N binding energies. We find a slight overbinding of about 300 keV for AV18, Nijm I and II . In case of CD Bonn the overbinding reaches about 800 keV. Note that adjusting the TM 3NF to 3H leads to a somewhat larger overbinding. Nevertheless all these numbers indicate that with those Hamiltonians one reaches the $`\alpha `$-particle binding energy rather closely and there is little room left for the action of 4N forces. On such a phenomenological level, however, it is not possible to decide about the need of 4N forces. It is always possible to add another piece of a 3N force which in this case should be repulsive to reach the $`\alpha `$-particle binding energy more accurately. For instance the $`\pi `$-$`\rho `$ exchange 3N force would provide at least one more parameter and both nuclei, 3He and 4He, could be described. In that respect the Urbana 3NF is adjusted in such a manner  that in addition to the 3N binding energy also nuclear matter is taken into account. This then fixes a repulsive piece in that 3NF, with the result that the $`\alpha `$-particle binding comes out essentially right. We repeated that calculation for the Urbana 3NF firstly performed with the GFMC method . Our result based on YE’s is given in Table II together with the previous one. There we also show a triton result by the CHH-method . We estimate our numerical accuracy to be about $`\pm 3`$keV ($`\pm 50`$keV) for the 3N (4N) system. There appear to be small differences between our and the GFMC results, especially in the kinetic energies.

Like in 3H one can also separate the 4He wave function into $`S`$, $`S^{}`$, $`P`$ and $`D`$-state probabilities. Here $`S`$ is spatially symmetric, $`S^{}`$ has two-dimensional mixed symmetry and $`P`$ and $`D`$ are the total $`L=1`$ and 2 orbital angular momentum parts . We compare the two nuclei in Tables III. For 4He $`S^{}`$ is reduced, $`P`$ somewhat enhanced and probabilities shifted from $`S`$ to $`D`$. It is remarkable that the TM 3NF together with CD-Bonn reduces the $`D`$-state probability, in contrast to all other cases.

In  it is argued that chiral symmetry requires that the $`c`$-term in the TM 3NF should be dropped, leading to a TM’ 3NF (one keeps the remaining constants $`b`$,$`d`$ unchanged and $`a`$ is replace by $`a^{}=a2c`$). As seen in Table II the resulting $`\alpha `$-article binding energy after fitting the cut-off to 3He coincides essentially with the $`\alpha `$-particle binding energy. It is also interesting to see that the $`c`$-term has a significant effect on the wave function as demonstrated in Table IV. There we show the expectation values of the four different terms contributing to the TM 3NF evaluated with the wave function including TM and the $`a`$,$`b`$ and $`d`$ expectation values based on the wave function including TM’. In relation to that latter wave function one can also evaluate the $`c`$-term. Interestingly it turns out to be repulsive, whereas for the TM wave function it is attractive. Also it is interesting to see that for TM the $`a`$-term is negligible, whereas it has some importance for TM’. These examples demonstrate again that 3NF’s change wave functions and cannot be treated perturbatively, a fact known since long time .

Summarizing, after adjusting 3NF’s to 3N binding energies the $`\alpha `$-particle binding energies based on modern nuclear forces are rather close to the experimental value. This indicates that 4N forces (at least for $`T=0`$) will be unimportant. More details on wave function properties will be published elsewhere. The FY equations are perfectly well under control for realistic NN and 3NF’s and will be a perfect tool to study upcoming new force structures given in chiral perturbation theory .

This work was supported financially by the Deutsche Forschungsgemeinschaft (A.N. and H.K.). The numerical calculations were performed on the CRAY T3E of the NIC in Jülich.


text:
# Strain versus stress in a model granular material: a Devil’s staircase.

## Abstract

The series of equilibrium states reached by disordered packings of rigid, frictionless discs in two dimensions, under gradually varying stress, are studied by numerical simulations. Statistical properties of trajectories in configuration space are found to be independent of specific assumptions ruling granular dynamics, and determined by geometry only. A monotonic increase in some macroscopic loading parameter causes a discrete sequence of rearrangements. For a biaxial compression, we show that, due to the statistical importance of such events of large magnitudes, the dependence of the resulting strain on stress direction is a Lévy flight in the thermodynamic limit.

The mechanical properties of granular media are currently an active field of research, both in the condensed matter physics and in the mechanics and engineering communities BJ97 ; WG97 ; HHL98 .

Granular packings close to mechanical equilibrium are traditionnally modelled, in the framework of continuum mechanics, with elastoplastic constitutive laws Muirwood ; HHL98 , *i.e.,* incremental stress-strain relations. Such laws, despite their practical success, were never clearly related to grain-level mechanics. Moreover, cohesionless granular systems seem to be quite different from ordinary solids. Many experimental, theoretical and numerical studies WG97 ; CLMNW96 ; RJMR96 ; OR97b ; MJN98  have recently been devoted to the peculiar features of stress transmission in granular systems at equilibrium, with correlations over length scales significantly larger than the grain size.

Observations of displacement fields and strains, as the system moves from one equilibrium to another, are scarcer. Systems of rigid grains are expected to deform because of rearrangements of the packing, rather than contact elasticity. How such rearrangements average to produce a macroscopic strain, related to stress variations, remains rather mysterious. The rather singular, unilateral form of the local interaction in such systems led some authors CWBC98  to expect quite unusual macroscopic properties, for which the very concept of strain, so familiar in mechanics of solids, would be irrelevant.

Direct grain-level approaches are, in principle, possible by numerical simulations. However, one has then to define a complete mechanical model to enable a calculation of particle trajectories. In practice, dynamical parameters ruling energy dissipation are often chosen according to computational convenience as much as physical accuracy. It would be desirable to assess the influence of such choices on the results.

The present numerical study addresses those problems, as follows.

Disordered, dense assemblies of rigid, circular, frictionless discs, are prepared by isotropic compaction. The force law reduces to the condition that contacts transmit repulsive normal forces of unknown magnitude. Then, the direction of the load is gradually altered, thus simulating the biaxial test of fig. 1. Exploiting the isostaticity property MO98a ; TW99 ; JNR2000  of such systems, which is exactly satisfied provided impenetrability is enforced accurately enough, we designed a prescription JNR2000  for the computation of sequences of equilibrium states, that we call the geometric quasi-static method (GQSM), in which the only inputs are the geometric data. The strain versus stress evolution in biaxial compression tests is recorded, and a statistical analysis of fluctuations and system size dependence is carried out. Then the GQSM predictions are compared with those of a standard molecular dynamics (MD) method.

First, samples of different sizes are generated (51 samples of N=1024 discs, 23 with N=1936, 10 with N = 3025 and 15 with N=4900). Disc diameters are uniformly distributed between 0.5 and 1 (the largest diameter is chosen as unit length), and packed in a loose state within a rigid, square box. After some amount of random mixing (using some dynamical method with energy conservation) we proceed to the isotropic compaction: two of the walls, 1 and 2 on fig. 1, are now mobile, and submitted to compressive external forces $`F_1`$ (along the x axis on the figure) and $`F_2=F_1`$ (along the y axis). Stress components $`\sigma _{11}=F_1/L_2`$ and $`\sigma _{22}=F_2/L_1`$ are kept constant, equal to $`p`$, while the system lengths along directions 1 and 2 ($`L_1`$ and $`L_2`$) decrease.

To produce a dense, isostatic equilibrium state we use the ‘lubricated granular dynamics’ method of refs. OR97b ; JNR97b . Then, the initial, reference configuration of the biaxial experiment is ready: $`F_2=(p+q)L_1^0`$ is gradually increased while $`F_1=pL_2^0`$ stays constant, strain components are defined as the relative decrease of lengths $`(L_i)_{1i2}`$, with reference to their initial values $`(L_i^0)_{1i2}`$ as $`ϵ_{ii}=\mathrm{\Delta }L_i/L_i^0`$. One also defines the volumetric 3D  strain as $`ϵ_v=ϵ_{11}ϵ_{22}`$. We use units such that $`p=1`$.

Our essential result here is the obtention of the $`ϵ(q)`$ curves in the thermodynamic limit, as loading parameter $`q`$ increases monotonically, at constant $`p`$. Let us first describe the GQSM procedure. One starts from an equilibrium state in which the force-carrying structure is isostatic. This means that the equilibrium conditions are sufficient to compute all contact force values, on the one hand (the structure is devoid of hyperstaticity, it is not ‘overbraced’ or ‘overconstrained’ MO98a ), and that the force-carrying structure is rigid (devoid of mechanisms or ‘floppy modes’ JNR2000 ) on the other hand. The first of these two properties stems from the condition that two grains need be exactly in contact to transmit a force to one another TW99 ; JNR2000 , and cannot interpenetrate. It entails that force values, once equilibrium positions are known, are geometrically determined, all material properties being irrelevant in the limit of rigid grains. The second property is satisfied, for stability reasons, because the grains are circular and contacts do not withstand tension JNR2000 . It entails that an assembly of rigid discs will not deform at all until some initially active contact opens. This cannot occur as long as contact forces are compressive, since this would require the potential energy to increase from equilibrium. As soon as one contact force vanishes, this contact will open JNR2000 , because the resulting motion corresponds to an instability. Hence the following algorithm:

(1) At equilibrium, as $`q`$ increases from its initial value $`q_0`$, the contact forces depend linearly on $`q`$ (equilibrium equations are linear). When $`q`$ reaches some value $`q_0+\delta q`$, the force vanishes in one contact, say $`l_0`$.

(2) Open $`l_0`$, all other contacts being maintained. Due to isostaticity, this entirely determines the initial direction of motion for the whole structure. Keep moving the grains with the same prescription.

(3) When another contact, say $`l_1`$, closes, the new contact structure (the old one, minus $`l_0`$, plus $`l_1`$) is isostatic and may carry the load with geometrically determined contact forces. If there is no traction, a new equilibrium state has been found: go back to step 1. Otherwise, pick up the largest traction, call the corresponding contact $`l_0`$ and go back to step 2.

This procedure determines a series of equilibrium configurations, that are separated by rearrangements occurring for discrete values of $`q`$. The strain versus stress curve is a staircase (see fig. 2).

As long as the same contacts carry the load, the system does not deform; as soon as a rearranging event occurs, strain variables jump to the values corresponding to the next equilibrium configuration. This algorithm clearly involves, in steps 2 and 3, an arbitrary ingredient: the prescription that contacts open one by one. The main merit of GQSM, however, is that it does not introduce parameters other than geometric ones.

We now focus on the rise of $`ϵ_{22}`$ with $`q`$, close to the origin, and ask whether the staircase approaches a smooth curve in the thermodynamic limit. (Its initial slope, if finite, would be the effective compliance of the material). To do so, one studies the statistics of stress ($`\delta q`$) and strain ($`\delta ϵ_{22}`$, $`\delta ϵ_{11}`$) steps.

Successive $`\delta q`$ and $`\delta ϵ`$ values are found to be independent, and the width $`\delta q`$ of a stability interval is not correlated to the following strain steps. Throughout the investigated $`q`$ interval, the probability distributions of increments $`\delta q`$, $`\delta ϵ_{22}`$, and $`\delta W=p\delta ϵ_vq\delta ϵ_{22}`$, which is the variation in potential energy corresponding to the current load, do not appreciably change KS . No significant difference between samples is observed either.

Both $`q`$ and $`ϵ_{22}`$ values reached at a given stage can thus be regarded as sums of equidistributed independent random increments. The distribution of stress increments $`\delta q`$ is displayed on fig. 3.

It decays exponentially for large $`\delta q`$, and is shifted to smaller and smaller values as $`N`$ increases, so that the probability distribution of the rescaled increment $`\delta qN^\alpha `$ is size-independent. We denote as $`\delta q_0`$ its average. The exponent can be estimated as $`\alpha =1.16\pm 0.04`$. Stability intervals shrink to zero as $`N`$ increases: in the thermodynamic limit, any macroscopic load variation entails some motion of the grains. This property of packings of frictionless rigid grains, known as *fragility* CWBC98 ; JNR2000 , is unambiguously established by our simulations. The value of $`\alpha `$ should be related to the shape of the force distribution for small values, and to the varying sensitivity of the contact forces to macroscopic load increments. The classical central-limit theorem, applied to $`q`$ increments, will relate for large systems the number of steps $`M`$ from the beginning of the biaxial compression to the current value of $`q`$, as

$$M\frac{q}{\delta q_0}N^\alpha .$$
(1)
The distribution of strain increments $`\delta ϵ_{22}`$ is shown on figure 4.

A reduced variable $`N^\beta \delta ϵ_{22}/\delta ϵ_0`$ can also be defined, with an $`N`$-independent probability distribution. Remarkably, the (power-law distributed) number $`c`$ of contact losses (or gains) in one strain step, which tends to grow with $`\delta ϵ_{22}`$, does not show any significant size dependence, and contacts that simultaneously open or close are homogeneously scattered throughout the sample, whatever $`N`$. We estimated $`\beta `$ as $`\beta =2.18\pm 0.12`$. Unlike for stress increments, the distribution now decays as a power law:

$$p(\delta ϵ_{22})(\delta ϵ_{22})^{(1+\mu )},$$

with $`\mu =0.46\pm 0.03`$. As $`\mu <1`$, it should be remarked that, although the *typical* strain increment decreases as $`N^\beta `$ in the limit of large samples, the *average* strain increment does not exist. The standard central limit theorem does not apply, but, due to the power law decay of the distribution function, one may resort to a generalized central limit theorem BG90 . In our case, for large numbers $`M`$ of increments, the asymptotic form of their sum $`ϵ_{22}`$ is

$$ϵ_{22}\delta ϵ_0N^\beta M^{\frac{1}{\mu }}\xi ,$$
(2)
with an $`M`$-independent random variable $`\xi `$ abiding by an asymmetric Lévy distribution of index $`\mu `$.

The stress-strain relationship is obtained on combining eqns 1 and 2:

$$ϵ_{22}=\frac{\delta ϵ_0}{(\delta q_0)^{\frac{1}{\mu }}}q^{\frac{1}{\mu }}N^{\beta +\frac{\alpha }{\mu }}\xi .$$
(3)
The presence of $`\xi `$ in relation 3 implies that the strain versus stress curve will never express a deterministic dependence. Some size effect – a dependence on $`N`$ – remains in the thermodynamic limit, unless one has:

$$\alpha =\beta \mu $$
(4)
Our estimates of $`\alpha `$, $`\beta `$ and $`\mu `$ are compatible with this relation. If it is satisfied, then the distribution of axial strain increments $`\mathrm{\Delta }ϵ_{22}`$ corresponding to a given, fixed, $`q`$ increment, $`\mathrm{\Delta }q`$, should no longer depend on the system size as soon as $`N`$ is large enough for $`\mathrm{\Delta }q`$ to involve, typically, sufficiently many elementary rearranging steps. This was checked, confirming (see fig. 5), within statistical errors, relation 4 and the absence of size effects.

In order to compare the GQSM results to the predictions of more conventional methods, we also simulated biaxial compressions on samples of 1024 and 3025 discs using MD, successively imposing deviator increments $`\mathrm{\Delta }q=10^3`$. Contacts obey an elastic unilateral law, with a normal stiffness equal to $`10^5`$. After each stress step, one waits for a new equilibrium state, requesting forces on each grain to balance up to an accuracy of $`10^5`$. Successive strain increments $`\mathrm{\Delta }ϵ_{22}^{MD}`$ are uncorrelated and distributed according to the same probability law, which coincides (within statistical uncertainties) as shown on fig. 5, with that of increments $`\mathrm{\Delta }ϵ_{22}`$ obtained for the same value of $`\mathrm{\Delta }q`$ by GQSM. Similarly, the(power-law) distribution of the *fraction* of the total number of contacts that change within one fixed increment $`\mathrm{\Delta }q`$ does not depend on $`N`$ (fig. 5).

We therefore conclude that, in a biaxial test, the axial strain $`ϵ_{22}`$ dependence on deviator $`q`$, as expressed by relations 3 and 4, is a Lévy stochastic process. It does not become deterministic in the thermodynamic limit, but it is devoid of size effect: the strain versus stress curve approaches a Devil’s staircase with a dense set, on the $`q`$ axis, of discontinuities of random magnitudes. Moreover, it appears not to depend on dynamical parameters (at least within the accuracy of the present study) and to be determined by the sole system geometry. These results apply, without appreciable change in the statistics, throughout the interval $`0q/p0.2`$.

The evolution of $`ϵ_v`$ (which can be of either sign, typically one order of magnitude smaller than $`ϵ_{22}`$) is somewhat more complicated, since successive increments $`\delta ϵ_v`$ are not equidistributed (unlike $`\delta W`$). Leaving a discussion of volumetric strains to a future publication, let us further comment here the behavior of axial 3D  strains.

Our results preclude the existence of a constitutive law in the usual sense. As $`q`$ is increased, $`ϵ_{22}`$ is typically of order $`q^{1/\mu }`$, but its actual value is essentially impredictible, with the remarkable consequence that macroscopic models for granular mechanics should be of a stochastic, rather than deterministic, nature. Of course, such a behavior might be limited to frictionless grains, and is thus perhaps more relevant for some colloidal glasses than for sand, although important fluctuations and high noise levels were sometimes reported for soil materials or glass beads. One does obtain an extremely noisy curve on numerically compressing our system at constant strain rate instead of controlling the stress: understandably, between two equilibrium states that can be relatively far apart, the contact structure does not consistently oppose any given level of stress, and the observed response depends on dynamical parameters. Laboratory tests, at constant strain rate, on such systems, would be sensitively influenced by the apparatus itself. It will be interesting to study the dependence of parameter $`\mu `$ on grain shape and polydispersity. Preliminary MD results on dense random packings of monodisperse spheres in 3D yielded Lévy-distributed large strain steps with a value of $`\mu `$ in the 0.4-0.6 range.

Physically, large strain increments are due to rearrangements involving many contact changes (GQSM steps 2 and 3 have to be repeated). Had we resorted to the *approximation of small displacements* (ASD), in which all relevant quantities are dealt with to leading order in the displacements, then, as shown in ref. JNR2000 , no iteration of steps 2 and 3 would have been necessary. Within the ASD, contacts are replaced one by one, and we have checked that the resulting strain increment distribution does possess an average value. This geometric approximation would ensure JNR2000  uniqueness of the equilibrium state and a one-to-one correspondence between stress and strain. The ASD was used for slightly polydisperse discs on a triangular lattice, in which case it can be justified and constitutive laws are obtained JNR2000 ; JNR97b  (the second ref. TW99  proposes an independent implementation of the GQSM, with the ASD). The distinctive mechanical features of model granular assemblies that are reported here are thus due to the motion of the system in configuration space along a tortuous path between local minima in a complex potential energy lanscape. Cooperative rearranging events are also observed in glassy relaxation DDKPHG98 , which might open interesting perspectives.

We are grateful to Jean-Philippe Bouchaud for useful suggestions.


text:
# On the fate of gas in ultraluminous infrared galaxies at low and high redshift

## 1 Ultraluminous galaxies at low redshift

Ultraluminous infrared galaxies (ULIGs) are peculiar galaxies with infrared luminosities in excess of 1012 L (Sanders & Mirabel 1996). These galaxies have very dense gas cores, as indicated by molecular line measurements, both of CO (Sanders, Scoville & Soifer 1991), and also of high-density traces like HCN (Solomon, Downes & Radford 1992) and CS (Solomon, Radford & Downes 1990). A number of measurements at radio, infrared, and optical wavelengths (Condon et al. 1991, Surace et al. 1998, Surace & Sanders 1999, Soifer et al. 1999) suggest that most of the very substantial bolometric luminosities of these galaxies comes from these dense molecular cores.

The existence of these extreme objects leads to the question: what do they evolve into? One approach to answering this question is to look for stellar systems at redshift $`z=0`$ with central densities (see Fig. 1 of Binggeli 1994 for measrements of stellar densities in local galaxies) of about 100 M pc-3 or greater, the gas densities of the ULIG cores. Kormendy & Sanders (1992) made this comparison and showed that the only realistic candidates on these grounds at $`z=0`$ are giant elliptical (gE) cores. Galaxy disks (which have lower surface-brightnesses than gE cores; Freeman 1970) do not have high enough stellar densities. These gE cores are dense but have profiles that are flatter than an extrapolation of a de-Vaucouleurs $`r^{\frac{1}{4}}`$ law to small radii; they are rarely isothermal (Lauer 1985). Recent $`Hubble`$ $`Space`$ $`Telescope`$ imaging (Faber et al. 1997) has confirmed the existence of such cores in the most luminous ellipticals (they are, however, less common in lower-luminosity ellipticals). Often the cores are kinematically decoupled from the rest of the galaxy (e.g. Forbes et al. 1996), suggesting that their formation happens in a different way to that of the rest of the galaxy. These cores also tend to have supermassive black holes (BHs) in their centers, as inferred from stellar-kinematical measurements. The bigger the galaxy, the bigger this BH (Magorrian et al. 1998, van der Marel 1999, Ferrarese & Merritt 2000, Gerbhardt et al. 2000). Additionally, the bigger the galaxy, the bigger the core (Lauer 1985).

An additional connection between the dense molecular cores of the ULIGs and the gE stellar cores + BHs is that their masses and sizes are similar, in addition to their densities. Typical masses for both are several times 109 M and typical radii are several hundred parsecs (see Sakomoto et al. 1999 for details of Arp 220, the nearest ULIG, and Lauer et al. 1985 for the core properties of gE cores). The BH contribution to the gE core masses can be added in from the correlations of Magorrian et al. 1998 or van der Marel 1999 which relate the BH mass to the galaxy mass, in combination with Tables 2 and 3 of Lauer 1985, which relate the core size to the size of the whole galaxy. Masses of the stellar cores can be computed directly from the core sizes, along with the mass-to-light ratio of Fukugita, Hogan & Peebles (1998).

But locally, gE galaxies are much more common than ULIGs. The number density of ULIGs with far-infrared luminosities in excess of 1012 L is about $`1.1\times 10^7\mathrm{Mpc}^3`$ (Saunders et al. 1990). Typical ULIGs have total $`K`$-band (2.2 $`\mu `$m) absolute magnitudes of $`M_K=25.4`$ (Carico et al. 1988, 1990; here, as throughout this paper we assume an Einstein-de Sitter cosmology with $`H_0`$ = 50 km s-1 Mpc-1), most of which comes from pre-existing stars in the progenitor galaxies (the dense cores are optically thick at $`K`$-band; Goldader et al. 1995). The number density of normal galaxies with $`M_K<25.4`$ is $`1.6\times 10^4\mathrm{Mpc}^3`$ (Szokoly et al. 1998), a factor of about 1500 larger than the density of ULIGs. Most of these $`K`$-band-luminous normal galaxies are early-type galaxies, which have the cores + BHs discussed in the previous two paragraphs; the contribution of luminous disks, which tend to be blue, to the $`K`$-band luminosity function at the very bright end is small (see e.g. Binggeli, Sandage & Tammann 1988 for the optical luminosity function of the different Hubble types, and Huang et al. 1997 for the optical-near-infrared colors).

Therefore if gE cores were made by ULIGs, the comoving density of ULIGs must have been a factor 1500 higher in the past than it is today, approximately 10-4 Mpc-3. But this is very similar to the comoving density (e.g. Trentham, Blain & Goldader 1999) of the SCUBA ($`>1`$ mJy at 850 $`\mu `$m in the observer frame) sources (Smail, Ivison & Blain 1997, Barger et al. 1998, Hughes et al. 1998, Eales et al. 1999), which are presumably high-redshift objects, and which Barger et al. (1998) argue are ULIGs based on their similar bolometric luminosities and spectral energy distribution (SEDs). Therefore a picture in which SCUBA sources make gE cores and their associated black holes seems very likely not just based on similarities of their SEDs and luminosities to those of local ULIGs but also due to the comoving number density of local gE cores \+ BHs being similar to that of the high-redshift SCUBA galaxies. It is this possibility that we investigate in the rest of this paper.

## 2 Ultraluminous galaxies at high redshift

In the previous section we argued that the SCUBA sources are probably the high-redshift analogues of low-redshift ULIGs, both in terms of their SED and luminosity concordance, and also in terms of the comoving number density concordance between the SCUBA galaxies and the likeliest remnants of the local ULIGs (the gE cores + BHs).

This SED (Ivison et al. 1998) and luminosity (Barger et al. 1999) concordance is observed for individual sources, and is also implied on statistical grounds. By this it is meant that: if the SCUBA sources were significantly hotter than local ULIGs, then they would overproduce the far-infrared background shortward of about 500 $`\mu `$m measured by Fixsen et al. (1998) using $`COBE`$. If they were signifiantly colder than local ULIGs, they would severely underpredict the background at 450 $`\mu `$m: the SCUBA sources are observed to generate at least one-third of the background measured by Fixsen et al. at 450 $`\mu `$m (Blain et al. 1999a). These two statements assume a median source redshift of approximately 3 (see the radio measurements of Smail et al. 2000) and then follow directly from the result that almost the entire background at 850 $`\mu `$m is generated by the SCUBA sources (Blain et al. 1999b, Barger, Cowie & Sanders 1999). Both statements are independent of any specific model of SCUBA source luminosity and density evolution.

Now let us make the following additional assumption, which we will call A1: The density structure of the gas in the SCUBA sources is the same as in local ULIGs. This is an assumption, not an observation, since we have no direct probes of the gas densities in the SCUBA sources. We only know about CO in two of the SCUBA sorces (Frayer et al. 1998, 1999), and have no measurements at all of high-density tracers like HCN or CS in these objects. There are two indirect pieces of evidence supporting A1. Firstly, HCN has been detected in the Cloverleaf quasar (Barvainis et al. 1997), a strongly-lensed infrared quasar which has an SED at far-infrared wavelengths similar to those of the SCUBA sources. Secondly, if the burst producing the high far-infrared luminosity is shortlived (see e.g. Blain et al. 1999a), then it is unlikely that it can be maintained in a coherent way over a galactic scale (in this case, the large scale over which the luminosity emerges from requires that the ultimate power source is star formation). Therefore the gas fuelling the burst is unlikely to be distributed over the whole galaxy. Scenarios do exist in which positive feedback can happen where explosions in one part of a galaxy can trigger those in another part (e.g. Taniguchi, Trentham & Shioya 1998), but these mechanisms are not powerful enough to generate the kinds of luminosity required here over a whole galaxy.

The important point is that if A1 is correct, then the above suspicion of the connection between local ULIGs and the SCUBA sources based on the (remnant) number density and SED concordance now becomes a very strong assertion. This is because the only things that the SCUBA sources can evolve into are elliptical galaxy cores and their associated black holes.

We do not consider the possibility that stars are made in the nucleus of the galaxies (as would be required by A1) and then redistributed throughout the galaxy in stellar-kinematical mergers because (1) the dense stellar cores so produced will not be disrupted by mergers (or secular evolution), (2) we would require about 100 mergers per L elliptical and the observed merger rate (e.g. Carlberg et al. 2000) is not that high, and (3) the number density concordance described above is then lost due to the large number of mergers.

Van der Marel 1999 finds the black hole mass (in solar units) $`\mathrm{log}_{10}M_{\mathrm{BH}}1.83+\mathrm{log}_{10}L_V`$. Therefore for a Schechter (1976) $`L^{}`$ elliptical galaxy ($`M_V=21.5`$) at redshift $`z=0`$, the central black hole has mass $`5.1\times 10^8\mathrm{M}_{}`$. For the same elliptical galaxy, the stellar core mass is $`M_c=1.6\times 10^9\mathrm{M}_{}`$, assuming that the core is approximately isothermal at small radii ($`r<r_c`$ where $`r_c`$ is the core radius), that the galaxy has a de Vaucouleurs $`r^{1/4}`$ profile at large radii, that the ratio of $`r_c`$ to the de Vaucouleurs effective radius $`r_e`$ is 0.033 (derived from the sample of Lauer 1985), and a value of $`\left(M/L_V\right)=3.86`$ for elliptical galaxies (see the discussion in Section 2.1.3 of Fukugita et al. 1998 – this assumes the stellar initial mass function (IMF) of Gould, Bahcall & Flynn 1996, which turns over at about 0.5 M, similar to the IMF of Kroupa, Tout & Gilmore (1993) that we adopt below). Therefore, for typical luminous elliptical galaxies

$$\frac{M_{\mathrm{BH}}}{M_c}\frac{1}{3}$$
(1)
Given the scatter in the correlations presented by van der Marel and Lauer, and the fact that many cores deviate significantly from isothermality, the most extreme galaxies may deviate from this ratio by as much an order of magnitude.

For accretion onto a black hole,

$$L_{\mathrm{acc}}=\eta \dot{M}c^2$$
(2)
where $`\eta `$ is the efficiency of accretion and $`\dot{M}`$ is the mass accretion rate.

For dust-enshrouded star formation where almost all the flux is absorbed and reradiated at far-infrared wavelengths, the star formation rate (in M yr-1) is (Rowan-Robinson et al. 1997)

$$\dot{M_{}}=1.2\times 10^{10}L_{}.$$
(3)
where $`L_{}`$ (in L) is the total luminosity generated by the young stars. This relation adopts a value of $`\varphi =0.45`$ appropriate to the stellar IMF of Kroupa et al. (1993), and a value of $`ϵ=1`$ (see Section 3 of Rowan-Robinson et al. 1997 for a definition of those parameters and a discussion of the derivation of this equation and other approaches). The IMF chosen here has been shown to be consistent with most current observations of both star-forming regions and the solar neighbourhood (Gilmore & Howell 1998).

Combining the two previous equations yields

$$L_{\mathrm{acc}}/L_{}1700\eta .$$
(4)
Accretion is a lot more efficient than star-formation at generating luminosity per unit mass: a 1012 L galaxy has a star formation rate of about 80 M yr-1 but a 2 $`\times `$ 1012 L quasar has an accretion rate of about 1 M yr-1 for $`\eta =0.1`$ (e.g. Chokshi & Turner 1992). Therefore averaged over time in the SCUBA sources:

$$<L_{\mathrm{acc}}/L_{}>_t\mathrm{\hspace{0.17em}560}\eta $$
(5)
This ratio is very big, meaning that the total power generated in the SCUBA sources comes almost entirely from accretion if A1 is correct. Note that throughout this calculation, we have been working backwards, deriving the properties of the progenitor from the properties of the remnants under the assumption that we can match particular progenitors to particular remnants; this is quite a different calculation from one in which we start out with a cloud of gas and attempt to track its progress.

The fact that the power sources are so heavily accretion-dominated means that most SCUBA sources must be dust-enshrouded active galactic nuclei (AGNs) if A1 is correct. This is not surprising, since the SCUBA sources contribute at least 10 per cent of the bolometric luminosity density of the Universe, and this could not be generated from the stars in gE cores alone, which comprise less than 1 per cent by mass of the stars in the Universe, and these are the only stars that the SCUBA sources are allowed to make if A1 is correct.

## 3 Additional evidence for the SCUBA source - AGN connection

There exist other reasons for suggesting a connection between the SCUBA sources and dust-enshrouded AGNs.

1. The local supermassive BH density (Magorrian et al. 1998, van der Marel 1999) is high and its production generates a bolometric background of

$$5.1\left(\frac{\rho _{\mathrm{AGN},\mathrm{dusty}}}{2\times 10^5\mathrm{M}_{}\mathrm{Mpc}^3}\right)\mathrm{nWm}^2\mathrm{sr}^1$$

(Trentham & Blain, 2000) if they radiate at one-tenth of the Eddington luminosity and if the luminosity density of obscured AGNs follows that measured by Boyle & Terlevich (1998) for optical quasars. This equation follows from a consideration of the energy released by accretion processes per comoving volume element given a final MDO density some part $`\rho _{\mathrm{AGN},\mathrm{dusty}}`$ of which was generated by dust-enshrouded accretion. The fiducial $`\rho _{\mathrm{AGN},\mathrm{dusty}}`$ is the estimate of Salucci et al. (1999) for obscured quasars. The SCUBA sources in Fig. 1 generate about 7 $`\mathrm{nWm}^2\mathrm{sr}^1`$, which is close to this value.

2. There is growing evidence that the hard (30 keV) X-ray background originates from absorbed quasars which reradiate energy absorbed at optical, ultraviolet, and soft X-ray at far-infrared and submillimetre wavelengths. Were this the case, Fabian & Iwasawa (1999) find a total reradiated energy density of about $`3\mathrm{nWm}^2\mathrm{sr}^1`$. This is close to the numbers in the previous paragraph, and it is therefore plausible that these same AGNs that we are hypothesizing are the SCUBA sources are also the objects which contribute to the hard X-ray background (see also Almaini et al. 1999 and Bautz et al. 2000). Note that if Compton-thick sources are very common (see Maiolino et al. 1998), the reradiated energy will be somewhat higher than $`3\mathrm{nWm}^2\mathrm{sr}^1`$.

3. The Madau plot – that is, the cosmic star-formation rate of the Universe as a function of redshift – when determined by ultraviolet and optical measurements alone (e.g. Steidel et al. 1999) and integrated over time, produces a total stellar density in critical units of $`\mathrm{\Omega }_{}0.004`$ (Pettini 1999), equal to the observed local value (Fukugita et al. 1998) Therefore an additional population of objects which are star-forming would lead to $`\mathrm{\Omega }_{}`$ being overproduced. Since $`>10`$ per cent of the bolometric luminosity of the Universe is generated in the SCUBA sources, were they to be star-forming, they would be producing a substantial fraction of $`\mathrm{\Omega }_{}`$. We therefore require an alternative power source, such as accretion onto a black hole. The Madau plot described above includes a correction for star-formation obscured by dust in normal galaxies, which are not the same as the ULIGs or SCUBA sources described in the previous sections: their predicted 850-$`\mu `$m fluxes, even for the most luminous normal star-forming galaxies, are substantially lower (see e.g. Trentham & Blain 2000).

4. All supermassive BHs at low $`z`$ are found in spheroids (Magorrian et al. 1998), except those in ULIGs like Mrk 231 which are probably in the process of forming spheroids (the core in the dissipative collapse, the rest of the galaxy in the stellar-kinematical merger of pre-existing stars). No supermassive BHs are known in pure disks. Such a scenario follows from the physical processes described in Section 2.

5. In their spectroscopic survey, Barger et al. (1999) find some dust-enshrouded AGNs. But it is clear that powerful Seyfert 1 galaxies are not present in substantial number, so that if this sample is AGN dominated, they must be very highly obscured. Such objects are not very common in local ULIG samples (see Section 3.2 of Lutz, Veilleux & Genzel 1999).

## 4 Problems with the SCUBA source - AGN connection

The previous two sections have provided some evidence in support of a picture where A1 is true. But there are problems with such a scenario:

1. The ratio of the number of ULIGs powered by accretion onto a BH to the number of ULIGs powered by star formation must change from about 1/3 at $`z=0`$ (Genzel et al. 1998) to about 80 at the redshift of the SCUBA sources (for $`\eta =0.1`$) if A1 is correct. This is very large effect, and there exists no clear physical mechanism to account for such a big difference, particularly since the low-redshift ULIGs and high-redshift SCUBA galaxies are similar in every other respect if A1 is true.

2. The $`q`$ value (relating the far-infrared flux to the 4.85 GHz radio flux; Condon, Frayer & Broderick 1991) of SMM J02399$``$0136 is high, suggesting that about half of its power (at most 3/4) comes from star formation (Frayer et al. 1998). Recent X-ray measurements of this galaxy with Chandra (Bautz et al. 2000) also suggest that the fraction of the bolometric luminosity contributed by star-formation is at least 20 %, probably more. This is much higher than the value of $`1/(560\eta +1)`$ predicted in Section 2. For this not to be a concern, SMM J02399$``$0136 must be a very unrepresentative example of a SCUBA galaxy (which is not the case given the optical spectroscopy of Ivison et al. 1998 and Barger et al. 1999) or its AGN must emit an amount of radio luminosity given its bolometric luminosity that somehow mimics starbursts.

3. If the SCUBA sources are powered by AGNs, the infalling gas must lose at least 99 % of its initial angular momentum before reaching the BH event horizon. This is not achiveable by any known mechanism. Suggested mechanisms include global gravitational or magnetic processes (Begelman 1994), but the details are unclear. Note that the radius of the orbit of the two nuclear gas disks in Arp 220 is about 250 pc (Sakamoto et al. 1999), which is approaching the size scale on which the central black hole dominates the gravitational potential. Although this has been a long-standing problem, the increased angular resolution of the gas achieved by current inteferometry gives it an additional importance if A1 is to be regarded as a plausible hypothesis.

One possible solution that simultaneously addresses all three problems above could be that: the ULIGs undergo an intense star formation phase in the center that generates a bolometric luminosity of $`10^{12}\mathrm{L}_{}`$ followed immediately by an AGN phase that lasts for a shorter time, but generates a higher bolometric luminosity, say about $`10^{13}\mathrm{L}_{}`$. If one selects objects instantaneously with bolometric luminosities above $`10^{12}\mathrm{L}_{}`$, one would end up with a sample dominated by starbursts, perhaps in a ratio similar to what Genzel et al. (1998) found locally if the relative timescales and luminosities between the starburst and AGN phase are extremely large and small, respectively. This could remain true even if the time-averaged luminosity was generated predominantly by AGNs, as required by eq. (5). Such a scenario could help the angular momentum problem since the black hole could accrete the stars themselves (Begelman 1994); given eq. (1) about one-fourth by mass of the total stars produced would have to be accreted. It is also consistent with a picture in which the most luminous things that exist at any redshift are AGN-powered (see e.g. Sanders & Mirabel 1996).

Another possible solution is that the fraction of ULIGs at $`z=0`$ that are AGNs has been underestimated because the AGNs are so heavily obscured that they are invisible even at mid-infrared wavelengths. The detection of hard X-ray emission from NGC 6240 (Vignati et al. 1999) suggests that the bolometric luminosity of this infrared-luminous galaxy, previously thought to be starburst-powered, originates from an AGN. If this object is at all typical of ULIGs, then the possibility of a much higher ratio of local ULIGs that are AGN-powered than one-third seems possible (see however the discussion in the final section of Lutz et al. 1999).

## 5 Discussion and summary

We have shown that if the SCUBA sources, which appear on SED and luminosity grounds to be ULIGs at high $`z`$ have the same gas density structure as local ULIGs, then they must generate most of their power from embedded AGNs. Such a scenario is attractive for many reasons: it generates the correct local number density of elliptical galaxy cores of the correct mass, the correct local mass density of their associated black holes, and gives a formation mechanism for kinematically decoupled cores (which probably form in some other way from the rest of the galaxy). These SCUBA sources would then be very different from local ULIGs, which derive most of their power from star formation.

But there are problems (see Section 4). If these cannot be resolved, then A1 must be abandoned. In this event, the gas in the SCUBA sources must be far more diffuse than in local ULIGs because in order to generate such huge luminosities, a substantial fraction of the stars in the galaxies must be made; in $`z=0`$ galaxies, large fractions of the stars in galaxies do not reside at extremely high densities of 100 M pc-3 or greater so it is improbable that they were made out of gas at these densities. In local ULIGs all the bolometric luminosity seems to come out of regions associated with gas at this high density.

Either way, there is a big difference between local ULIGs and high-redshift SCUBA sources. This is the main result of this letter.

Testing the validity of A1 will be an important exercise in the future. HCN measurements (like for the Cloverleaf quasar) would be best, but until these and other millimetre wave spectral diagonstics (other than CO) are available, the best hope is to rely on indirect methods like those outlined in the letter.


text:
# Inverse meson mass ordering in color-flavor-locking phase of high density QCD: erratum

(April 2000)

## Abstract

We correct a mistake in the calculation of meson masses at large baryon chemical potential $`\mu `$ made in Ref. .

The values of the meson masses reported in Ref.  are not correct. The correct answer for the meson masses and mixings are given by Eqs. (45) and (46) in Ref.  with $`C=2c/f_\pi ^2`$, and

$$c=\frac{3\mathrm{\Delta }^2}{2\pi ^2},$$
(1)
where $`\mathrm{\Delta }`$ is the value of the superconducting gap at the Fermi surface. To leading order in perturbation theory, the constant denoted in Ref.  as $`c^{}`$ vanishes. Parametrically, the meson masses are of order $`m_q\mathrm{\Delta }/\mu `$ instead of $`m_q`$ as reported in Ref. . Moreover, since $`c^{}`$ vanishes to the leading order, a truly inverse mass hierarchy emerges instead of a partial one previously reported. This means the inverse hierarchy anticipated by the heuristic argument following Eq. (46) in Ref.  is indeed realized.

Technically, the mistake of our previous calculation of the masses came from the failure to take into account the fact that the quark mass insertion in a quark line changes the quark mass shell from $`E=|𝒑|\mu `$ to $`E=|𝒑|\mu `$ . As a result, an almost on-shell quark near the Fermi surface becomes an off-shell particle with virtuality $`2\mu `$.

The correct computation of the constants $`c`$ and $`c^{}`$ follows the same methodology of Ref.  of matching vacuum energies and involves the evaluation of two Feynman diagrams in Fig. 1, which both make leading-order contributions to the shift of the vacuum energy.

Note that the gap insertions, which reverse the direction of the quark lines in Fig. 1, occur only on the propagators which are near the particle mass shell. The lower part of the first diagram between $`M`$ insertions corresponds to the one-loop quark self-energy and is equal, by the gap equation, to a quantity commonly referred to as the “anti-particle gap”. The latter is, however, not gauge invariant even to leading order in perturbation theory. Consequently, the first diagram is not gauge-fixing independent on its own. However, the sum of the two diagrams is indeed gauge-fixing independent, which is a consequence of Ward-Takahashi identity and can be verified by a direct calculation.

The following comments concerning the result (1) are in order. The contribution of the magnetic gluon exchange cancels between the two diagrams. This happens because emission of a magnetic gluon by an on-shell quark is suppressed by parity conservation. As a result, $`c`$ is proportional to $`\mathrm{\Delta }^2`$ instead of $`\mathrm{\Delta }^2\mathrm{log}(\mu /\mathrm{\Delta })`$. The remaining electric gluon exchange is suppressed at small gluon momenta because of the angular momentum conservation and the requirement that the helicity flip must occur in the diagram. This removes a collinear divergence which would have given rise to an extra factor of $`1/\mathrm{log}g`$ .