text:
# A First Estimate of the Baryonic Mass Function of Galaxies

## 1 Introduction

The distribution of the mass in collapsed baryons (cold gas and stars) in galaxies is a fundamental prediction of galaxy formation models. Unfortunately, to date there is no robust estimate of the baryonic mass function (MF) of galaxies, leaving modelers with the non-trivial task of predicting stellar masses or, even worse, galaxy luminosities. Discrepancy between the model and data may indicate a problem with the predicted distribution of galaxy baryonic masses, or could represent poorly-constrained star formation (SF), stellar population or dust prescriptions. In this Letter, we present a first estimate of the baryonic MF of galaxies by assigning galaxies gas fractions statistically (based on an independent sample), under the assumption of a universally-applicable stellar initial mass function (IMF).111A time-varying IMF, as speculated on by Hernandez & Ferrara (2001) or Ferguson, Dickinson, & Papovich (2002), would invalidate this assumption.

The time is ripe to attempt this for the first time. With the advent of large, relatively complete surveys, the luminosity function (LF) is now well-constrained in the optical and near-infrared (NIR) (Gardner et al., 1997; Cole et al., 2001; Kochanek et al., 2001; Norberg et al., 2002; Liske et al., 2003; Blanton et al., 2003; Bell et al., 2003). Furthermore, under the assumption of a universally-applicable stellar IMF, the distribution of stellar masses is reasonably well-constrained, with an overall normalization uncertainty caused by our relatively poor knowledge of the faint end slope of the IMF (Cole et al., 2001; Bell et al., 2003). Crucially, there are also relatively large samples of galaxies with $`K`$-band data and gas masses, allowing a reasonably accurate characterization of the gas mass of galaxies as a function of their physical parameters (Bell & de Jong, 2000).

## 2 The Data and Methodology

Because of the lack of a large galaxy survey with both gas mass and $`K`$-band data, we take a sampling approach, analogous to that used by Loveday (2000) to estimate the $`K`$-band luminosity function from a $`B`$-band limited survey. Essentially, we estimate a stellar MF (§2.1) and then add representative gas masses to each galaxy (§2.2), allowing us to estimate the distribution of galaxy baryonic masses (§3). We assume $`\mathrm{\Omega }_\mathrm{m}=0.3`$, $`\mathrm{\Omega }_\mathrm{\Lambda }=0.7`$, and $`H_0=100h`$ km s-1 Mpc-1.

### 2.1 Estimating the Stellar Mass Function

We construct the baryonic MF using a combined sample of galaxies from the Two Micron All Sky Survey (2MASS; Skrutskie et al., 1997) and the Sloan Digital Sky Survey (SDSS; York et al., 2000). We use the SDSS early data release (EDR; Stoughton et al., 2002) to provide an 84% redshift complete $`r17.5`$ sample of galaxies with accurate $`ugriz`$ fluxes over 414 square degrees, which is $``$10% less than the whole EDR imaging area because some spectroscopic plates that were not attempted (Stoughton et al., 2002). The 84% redshift completeness within this area is our own direct estimate based on the fraction of galaxies fulfilling the Strauss et al. (2002) criteria that have spectra, in agreement with the EDR analysis of Nakamura et al. (2003). To account for light missed by the Petrosian magnitude estimator (Strauss et al., 2002; Blanton et al., 2003), we add 15% to the fluxes of galaxies morphologically classified as early-type using the SDSS $`r`$-band concentration parameter following Strateva et al. (2001). This correction produces only a $`5`$% effect on the LFs and MFs (Blanton et al., 2003; Bell et al., 2003). We also correct for an $`8`$% overdensity of galaxies in the EDR, as estimated by comparing the number of $`10K13.5`$ galaxies in the EDR spectroscopic area with that from the sky with $`|b|30\mathrm{deg}`$.

We use the now complete 2MASS extended source and point source catalogs to augment the SDSS $`ugriz`$ fluxes with $`K`$-band fluxes, and for extended sources $`K`$-band half-light radii. We correct 2MASS $`K`$-band fluxes to total following a comparison with deeper $`K`$-band data from Loveday (2000); for extended sources this amounts to a 0.1 mag correction (Bell et al., 2003). We do not use 2MASS $`J`$ or $`H`$-band data here because we cannot correct the magnitudes similarly. The optical and NIR magnitude zero points are accurate to $``$0.05 and $``$0.02 mag respectively, and the random errors are 0.05 mag (optical) and 0.2 mag (NIR).

To estimate $`k`$-corrections, evolution corrections, and the present day stellar mass-to-light ratios (M/Ls), we fit the $`ugrizK`$ observed fluxes222Not all galaxies have $`ugrizK`$ fluxes. We have checked that missing passbands do not significantly bias the estimated $`k`$-corrections, evolution corrections or stellar M/Ls (but do, of course, increase the random error somewhat). with model stellar populations. These populations have a range of metallicities and SF histories at a given redshift. We use the Pégase stellar population synthesis model (see Fioc & Rocca-Volmerange, 1997, for a description of an earlier version of their model) with a ‘diet’ Salpeter IMF (following Bell & de Jong, 2001) that has the same colors and luminosity as a normal Salpeter IMF, but with only 70% of the mass (due to a smaller number of low-mass stars). Corrections derived using this technique are consistent with those used by Blanton et al. (2003). The stellar M/Ls we derive are within 10% of those from the spectral modeling technique of Kauffmann et al. (2003), accounting for differences in IMF; the random and systematic uncertainties from dust and bursts of SF dominate, however, and are $`25`$% (Bell & de Jong, 2001). This IMF is ‘maximum disk’, inasmuch as IMFs richer in faint low-mass stars over-predict the rotation velocity of Ursa Major Cluster galaxies with $`K`$-band photometry and well-resolved Hi rotation curves. This prescription thus gives the maximum possible stellar M/L. Naturally, a different choice of IMF allows lower M/Ls. For example, the popular Kennicutt or Kroupa IMFs have $``$37% lower M/Ls than this IMF, and are thus ‘submaximal’ (see Bell & de Jong, 2001, for more discussion of this point).

We calculate LFs and stellar MFs using the $`V/V_{\mathrm{max}}`$ formalism (Felten, 1977), taking into account foreground Galactic extinction, $`k`$-corrections, and evolution corrections. In Bell et al. (2003), we match precisely published LFs; in particular, we reproduce the $`g`$-band and $`K`$-band LF and luminosity densities to within $`10`$% (Blanton et al., 2003; Cole et al., 2001; Kochanek et al., 2001). Furthermore, this method produces accurate stellar MFs that match the estimate of Cole et al. (2001) to $``$5% in total stellar mass density (accounting for IMF differences), but can do so using LFs limited by optical or NIR magnitude limits (Bell et al., 2003). For this Letter, we choose 11848 galaxies with $`13r17.5`$ and $`g17.74`$, which ensures that we have accurate $`gr`$ color estimates providing a stellar M/L accuracy of better than 25%, while avoiding potential biases against low surface brightness galaxies in 2MASS (Bell et al., 2003). The stellar MF estimated using this technique is shown in Fig. 2, along with the stellar MF of Cole et al. (2001) for comparison. A much more extensive description of the LF and stellar MF construction is given by Bell et al. (2003).

### 2.2 Estimating Gas Masses

Because there are no samples of galaxies with good number statistics, deep optical/NIR data and gas masses, we estimate the gas masses of SDSS$`+`$2MASS galaxies indirectly. We use galaxies from Bell & de Jong (2000) with $`K`$-band luminosities, half-light radii and gas masses to statistically assign a gas mass to every SDSS$`+`$2MASS galaxy, appropriate to its $`K`$-band luminosity and half-light radius.

Fig. 1 shows the $`K`$-band half-light radii and luminosities for the late-type subsample of the SDSS$`+`$2MASS galaxies (contours) and for the comparison sample of 156 galaxies with gas masses (filled circles) taken from Bell & de Jong (2000). We estimate gas masses by multiplying the Hi gas mass by 1.33 to account for helium, and by the morphological type-dependent H2 to Hi ratio presented by Young & Knezek (1989). For the 21% of the galaxies without $`K`$-band luminosities and sizes, we adopt the $`r`$-band half-light radii and estimate the $`K`$-band luminosity by dividing the $`g`$-band-derived stellar mass by the model $`K`$-band stellar M/L. Our assumptions are accurate to better than 40% in both cases, and make no difference to our results.

Each SDSS$`+`$2MASS galaxy is assigned a gas mass from a galaxy in the Bell & de Jong (2000) comparison sample with similar half-light radius and $`K`$-band luminosity. For galaxies morphologically classified as late-type we assign gas masses using a randomly chosen comparison galaxy that is within a factor of two in size, and within one magnitude in $`K`$-band luminosity. Galaxies are morphologically classified using the $`r`$-band concentration parameter following Strateva et al. (2001). We scale the gas mass by the difference in luminosity to conserve the galaxy gas fraction. As a consistency check, we also assign gas masses by choosing the nearest galaxy in half-light radius–luminosity space, and by assigning gas to galaxies of all morphological types (see Fig. 2). These changes make no appreciable difference to our results.

## 3 The Baryonic Galaxy Mass Function

Our baryonic galaxy MF is shown as the solid line with open circles in the left panel of Fig. 2. Shown for comparison is the stellar MF from Cole et al. (2001), and the $`g`$-selected stellar MF described above and in Bell et al. (2003). The baryonic galaxy MF follows the stellar MF at high masses, shows a modest offset at the ‘knee’ of the MF, and shows a reasonably steep faint-end slope. This behavior is expected, as low-mass field galaxies tend to have high cold gas fractions and more ongoing SF (Bell & de Jong, 2000; Kauffmann et al., 2003b), which steepens the baryonic MF compared to optical/NIR LFs that typically have $`\alpha 1.0`$. The Schechter Function fit parameters are $`\varphi ^{}h^3=0.0108(6)`$ Mpc-3 log$`{}_{10}{}^{}M_{}^{1}`$, $`M^{}h^2=5.3(3)\times 10^{10}M_{\mathrm{}}`$, and $`\alpha =1.21(5)`$, where the formal error estimates are in parentheses.

In the middle panel of Fig. 2, we show the effects of the systematic uncertainties. In particular, the dashed, dotted and dot-dashed lines show the effects of changing the gas mass assignment (see the figure caption for more details); these uncertainties have only small effects. The two bare solid lines show the effect of assuming different stellar IMFs. As stated earlier, we adopt as our default an IMF that has the largest stellar M/L permitted, without over-predicting the rotation velocities of spiral galaxies in the Ursa Major Cluster (solid line with open circles; Bell & de Jong, 2001). Yet, the stellar M/L may be lower than this maximal value; thus, we show two cases. First, we plot (upper solid line) the increasingly popular Kroupa, Tout, & Gilmore (1993) and Kennicutt (1983) IMFs, which both have M/Ls of $``$70% of maximal IMF ($``$50% Salpeter). The other case is an IMF that has M/Ls of only 40% of our maximal IMF (the lowest solid line), which corresponds roughly to the 63% disk velocity contribution to the rotation curve as argued by Bottema (1997). It is clear that, assuming a universal IMF, the factor-of-two uncertainty in the stellar M/L dominates the error budget in terms of total cold baryonic mass in the local Universe. Given the stellar IMF, the systematic uncertainties are $`25`$%. We estimate these uncertainties using different passbands to estimate stellar mass, using different gas mass assignment methods, and accounting for the effects of dust and bursts of SF on the M/Ls (Bell & de Jong, 2001; Bell et al., 2003). The Schechter Function fit parameters for the Kennicutt/Kroupa case and the Bottema case are: $`\varphi ^{}h^3=0.0116(5),0.0142(8)`$ Mpc-3 log$`{}_{10}{}^{}M_{}^{1}`$, $`M^{}h^2=3.78(11),2.24(8)\times 10^{10}M_{\mathrm{}}`$, and $`\alpha =1.22(3),1.20(3)`$.

It is important to make sure that our statistical procedure assigns gas masses consistent with the true galaxy population by comparing with the Hi or H2 MF of galaxies (the right hand panel of Fig. 2). The solid black line with open circles and error bars is our Hi MF derived in this way. The dashed and dotted lines show the effects of using the closest galaxy to estimate gas mass and allowing elliptical galaxies to have gas mass, respectively. For comparison, the Hi MF of the blind Hi Aricebo survey of Rosenberg & Schneider (2002) is plotted as the smooth solid curve. We also show our prediction of the H2 galaxy MF as the lower solid grey line with open circles and error bars. For comparison, we plot the Schechter Fit to the first observational estimate of the H2 MF (Keres, Yun, & Young, 2003). Given the uncertainties, it is clear that our method for estimating galaxy gas masses reproduces the Hi and H2 galaxy MFs with astonishing precision. Therefore, the sample of galaxies that we use to assign gas masses is indeed reasonably representative, and lends considerable credibility to our estimate of the baryonic galaxy MF.

## 4 Discussion

Even with the factor-of-two uncertainty from the contribution of low-mass stars to the overall stellar M/L, we can still draw some conclusions about the local Universe. It is clear that the overall efficiency of galaxy formation is very low. Firstly, the faint end slope of the baryonic MF is $`1.2`$, which is much shallower than the $`2`$ expected for the halo MF (e.g., White & Frenk, 1991). Secondly, integrating under the MF, we derive $`\mathrm{\Omega }_{\mathrm{cold}\mathrm{baryon}}h=2.4_{1.4}^{+0.7}\times 10^3`$, including the IMF and 25% systematic stellar M/L uncertainties. Our estimate agrees well with the value of $`2.9\pm 1.5\times 10^3`$ from Fukugita et al. (1998), and is preferred due to our better accounting for stellar M/Ls compared with Fukugita et al. (1998) who use (harder to convert into stellar mass) $`B`$-band luminosity densities assuming a similar IMF to the maximum-disk IMF we adopt here. Taking the value of the total baryon density from the Big Bang nucleosynthesis value of O’Meara et al. (2001), and assuming h=0.7$`\pm `$0.07, we find $`\mathrm{\Omega }_{\mathrm{cold}\mathrm{baryon}}/\mathrm{\Omega }_\mathrm{b}8_5^{+4}`$%, where the error estimates account for the uncertainties in IMF, H0, $`\mathrm{\Omega }_\mathrm{b}`$, our gas assignment method, and the $``$25% uncertainties in stellar M/Ls from dust and bursts of SF. Our value is quite consistent with the low galaxy formation efficiency characteristic of most current models, which have low efficiencies at the low and high-mass ends because of feedback from supernovae and inefficient gas cooling, respectively (e.g., Cole et al., 2000).

Accounting for the possible gaseous content of elliptical galaxies, for sub-maximal M/Ls, and for the effects of dust and bursts of SF on stellar M/Ls, the universal gas fraction, $`f_g=\mathrm{\Omega }_{\mathrm{cold}\mathrm{gas}}/\mathrm{\Omega }_{\mathrm{cold}\mathrm{baryons}}`$, should lie in the range $`0.2f_g0.5`$. For the ‘maximal’ IMF, we find $`f_g25`$%. Fukugita et al. (1998) and Keres, Yun, & Young (2003) find values of 15–20% when our IMF is adopted; their slightly lower determinations stem primarily from a lower estimate of Hi mass density. Nevertheless, all the studies agree that $`f_g0.5`$; therefore, the dynamically cold baryons (i.e., the gas and stars in disks and spheroids) are primarily in the form of stars, even for low stellar M/Ls.

It is well-known observationally that cluster optical/NIR LFs have steeper faint-end slopes than field LFs (e.g., Trentham & Tully, 2002). Furthermore, cluster galaxies tend to have little ongoing SF and little gas, so that most cluster galaxies are star-dominated with large stellar M/Ls (e.g., Kuntschner, 2000). Thus, the trend of increasing faint end slope with increasing cluster mass noted by e.g., Trentham & Tully (2002) may be more naturally interpreted as a constant baryonic MF, with a suppression of recent SF in massive clusters of galaxies. Obviously, a deeper investigation of this issue is warranted before speculating any further.

## 5 Conclusions

Together with the baryonic (stellar$`+`$cold gas) luminosity–linewidth relation (e.g., McGaugh et al., 2000; Bell & de Jong, 2001), the baryonic galaxy MF is an ideal test of models of galaxy formation and evolution. In this Letter, we have estimated the baryonic galaxy MF in the local Universe for the first time assuming a universally-applicable stellar IMF. We assign gas fractions statistically to a large sample of galaxies from 2MASS and SDSS, using a local sample with accurate $`K`$-band and gas fraction data. We cross-check this statistical procedure against independent Hi and H2 surveys, finding excellent agreement. The baryonic MF is similar to the stellar MF at the high mass end (with a slightly higher density normalization), and has a reasonably steep faint end slope, $`\alpha 1.2`$, due to the typically high cold gas fractions and low stellar M/Ls of low-mass galaxies. Integrating under the baryonic MF, we find that the baryonic (stellar$`+`$cold gas) mass density implied by this estimate is $`\mathrm{\Omega }_{\mathrm{cold}\mathrm{baryon}}=2.4_{1.4}^{+0.7}\times 10^3`$, or 8$`{}_{5}{}^{}{}_{}{}^{+4}`$% of the Big Bang nucleosynthesis expectation. This clearly implies a low overall efficiency of galaxy formation.

E. F. B. was supported by the European Research Training Network Spectroscopic and Imaging Surveys for Cosmology. D. H. M. and M. D. W. acknowledge support by JPL/NASA through the 2MASS core science projects. This publication makes use of data products from the Two Micron All Sky Survey which is a joint project of the University of Massachusetts and the Infrared Processing and Analysis Center/California Institute of Technology, funded by the National Aeronautics and Space Administration and the National Science Foundation. This publication also makes use of the Sloan Digital Sky Survey (http://www.sdss.org/). We thank the referee, Andrea Ferrara, for his helpful suggestions and feedback.


format: text

text:
# Probe Initial Parton Density and Formation Time via Jet Quenching

## 1 Introduction

In high-energy heavy-ion collisions, one of the pressing issues is the direct determination of the initial parton density of the dense medium and its formation time. Jet quenching or attenuation of leading particles from jet fragmentation can provide an effective tool to measure directly the parton density of the medium with which the jet interacts strongly during its propagation. Recent theoretical studies  have shown that parton energy loss via induced radiation is directly related to the gluon density of the medium. The attenuation will suppress the final leading hadron distribution giving rise to modified parton fragmentation functions . By measuring the medium modification of the fragmentation function one can thus extract the effective parton energy loss. Most importantly, one can compare the effective parton energy loss extracted from heavy-ion collisions to that of cold nuclei and extract the initial parton density of the hot dense medium relative to a cold nucleus.

## 2 Parton Energy Loss in Cold Nuclei

In deeply inelastic scattering (DIS) off a nuclear target, a quark suffers multiple scattering with other nucleons inside the nucleus after it was knocked out of its parent nucleon. The induced gluon bremsstrahlung before hadronization leads to suppression of leading hadrons and gives rise to a modified quark fragmentation function. Including the leading twist-4 contributions from double scattering processes, the nuclear medium induced correction to the fragmentation function is found to be similar to the radiative correction in the vacuum except that the normal splitting functions are replaced by the modified ones. These modified splitting functions are found  to depend on the quark-gluon correlation function $`T_{qg}^A(x,x_L)/q_A(x)`$. Here

$`T_{qg}^A(x,x_L)`$ $`=`$ $`{\displaystyle \frac{dy^{}}{2\pi }𝑑y_1^{}𝑑y_2^{}e^{i(x+x_L)p^+y^{}+ix_Tp^+(y_1^{}y_2^{})}(1e^{ix_Lp^+y_2^{}})(1e^{ix_Lp^+(y^{}y_1^{})})}`$ (1)
$`{\displaystyle \frac{1}{2}}`$ $`A|\overline{\psi }_q(0)\gamma ^+F_\sigma ^+(y_2^{})F^{+\sigma }(y_1^{})\psi _q(y^{})|A\theta (y_2^{})\theta (y^{}y_1^{})`$
and $`q_A(x)`$ is the twist-two quark distribution of the nucleus. Such a two-parton correlation function contains essentially four independent twist-4 parton matrix elements in a nucleus \[$`x_T=k_T^2/2p^+q^{}z(1z)`$\]. The dipole-like form-factor $`(1e^{ix_Lp^+y_2^{}})(1e^{ix_Lp^+(y^{}y_1^{})})`$ arises from the interference between the final state radiation of the $`\gamma ^{}q`$ scattering and the gluon bremsstrahlung induced by the secondary quark-gluon scattering. By generalizing the factorization assumption to these twist-four parton matrices, we have

$$T_{qg}^A(x,x_L)/q_A(x)\stackrel{~}{C}(Q^2)m_NR_A(1e^{x_L^2/x_A^2}),$$
(2)
with a Gaussian nuclear distribution $`\rho (r)\mathrm{exp}(r^2/2R_A^2)`$, $`R_A=1.12A^{1/3}`$ fm. Here, $`x_A=1/m_NR_A`$, and $`m_N`$ is the nucleon mass.

Since the two interference terms in the dipole-like form-factor involve transferring momentum $`x_Lp^+`$ between different nucleons inside a nucleus, they should be suppressed for large nuclear size or large momentum fraction $`x_L`$. Notice that $`\tau _f=1/x_Lp^+`$ is the gluon’s formation time. Thus, $`x_L/x_A=L_A/\tau _f`$, with $`L_A=R_Am_N/p^+`$ being the nuclear size in the infinite momentum frame. The effective parton correlation and the induced gluon emission vanishes when the formation time is much larger than the nuclear size, $`x_L/x_A1`$, because of the Landau-Pomeranchuck-Migdal (LPM) interference effect. Therefore, the LPM interference restricts the radiated gluon to have a minimum transverse momentum $`\mathrm{}_T^2Q^2/m_NR_AQ^2/A^{1/3}`$. The nuclear corrections to the fragmentation function due to double parton scattering will then be in the order of $`\alpha _sA^{1/3}/\mathrm{}_T^2\alpha _sA^{2/3}/Q^2`$, which depends quadratically on the nuclear size. For large values of $`A`$ and $`Q^2`$, these corrections are leading; yet the requirement $`\mathrm{}_T^2Q^2`$ for the logarithmic approximation in deriving the modified fragmentation function is still valid.

Shown in Fig. 2 are the calculated nuclear modification factor of the fragmentation functions for $`{}_{}{}^{14}N`$ and $`{}_{}{}^{84}Kr`$ targets as compared to the recent HERMES data . The predicted shape of the $`z`$\- and and energy dependence agrees well with the experimental data. A remarkable feature of the prediction is the quadratic $`A^{2/3}`$ nuclear size dependence, which is verified for the first time by an experiment. The only parameter in our calculation is found to be $`\stackrel{~}{C}(Q^2)=0.0060`$ GeV2 with $`\alpha _\mathrm{s}(Q^2)=0.33`$ at $`Q^23`$ GeV2.

If one defines theoretically the quark energy loss as that carried by the radiated gluons, then the averaged total energy loss is,

$$\mathrm{\Delta }E=\nu \mathrm{\Delta }z_g\stackrel{~}{C}\alpha _\mathrm{s}^2(Q^2)m_NR_A^2(C_A/N_c)3\mathrm{ln}(1/2x_B).$$
(3)
With the determined value of $`\stackrel{~}{C}`$, $`x_B0.124`$ in the HERMES experiment  and the average distance $`L_A=R_A\sqrt{2/\pi }`$ for the assumed Gaussian nuclear distribution,one gets the quark energy loss $`dE/dL0.5`$ GeV/fm inside a $`Au`$ nucleus.

One can also calculate the transverse momentum broadening of the quark jet which is also related to a similar twist-four parton matrix elements.

$$\mathrm{\Delta }q_{}^2=\frac{2\pi \alpha _\mathrm{s}}{N_c}\frac{\stackrel{~}{T}_{qg}^A(x_B)}{q_A(x_B)}\stackrel{~}{C}\frac{\pi \alpha _\mathrm{s}}{N_c}m_NR_A0.011A^{1/3}\mathrm{GeV}^2$$
(4)
## 3 Initial Parton Density in $`Au+Au`$ at RHIC

To extend our study to jets in heavy-ion collisions, we assume $`k_T^2\mu ^2`$ (the Debye screening mass) and a gluon density profile $`\rho (y)=(\tau _0/\tau )\theta (R_Ay)\rho _0`$ for a 1-dimensional expanding system. Since the initial jet production rate is independent of the final gluon density which can be related to the parton-gluon scattering cross section  \[ $`\alpha _sx_TG(x_T)\mu ^2\sigma _g`$\], one has then

$$\frac{\alpha _sT_{qg}^A(x_B,x_L)}{f_q^A(x_B)}\mu ^2𝑑y\sigma _g\rho (y)[1\mathrm{cos}(y/\tau _f)],$$
(5)
where $`\tau _f=2Ez(1z)/\mathrm{}_T^2`$ is the gluon formation time. One can recover the form of energy loss in a thin plasma obtained in the opacity expansion approach ,

$$\frac{dE}{dL}\frac{\pi C_aC_A\alpha _s^3}{R_A}_{\tau _0}^{R_A}𝑑\tau \rho (\tau )(\tau \tau _0)\mathrm{ln}\frac{2E}{\tau \mu ^2},$$
(6)
where we have assumed $`\sigma _gC_a2\pi \alpha _s^2/\mu ^2`$ ($`C_a`$=1 for $`qg`$ and 9/4 for $`gg`$ scattering).

Neglecting the logarithmic dependence on $`\tau `$, the averaged energy loss in a 1-dimensional expanding system can be expressed as

$$\frac{dE}{dL}_{1d}\frac{dE_0}{dL}\frac{2\tau _0}{R_A};\frac{dE_0}{dL}=\frac{\pi C_aC_A\alpha _s^3}{2}\rho _0R_A\mathrm{ln}\frac{2E}{\tau _0\mu ^2}$$
(7)
where $`dE_0/dL`$ is the energy loss in a static medium with the same gluon density $`\rho _0`$ as in a 1-d expanding system at an initial time $`\tau _0`$. Because of the expansion, the averaged energy loss $`dE/dL_{1d}`$ is suppressed as compared to the static case and does not depend linearly on the system size. Similarly, one have the total transverse momentum broadening in an expanding system,

$$\mathrm{\Delta }q_{}^2=2\pi C_a\alpha _s^2_{\tau _0}^R𝑑\tau \rho (\tau )=\mathrm{\Delta }q_{}^2_0\frac{\tau _0}{R_A}\mathrm{ln}\frac{R_A}{\tau _0};\mathrm{\Delta }q_{}^2_0=2\pi C_a\alpha _s^2\rho _0$$
(8)
Fitting the PHENIX data  on suppression of large $`p_T`$ $`\pi _0`$ spectra in Fig. 2 yields $`dE/dL_{1d}0.34\mathrm{ln}E/\mathrm{ln}5`$ GeV/fm. Taking into account the expansion, this would be equivalent to $`(dE/dL)_0=0.34(R_A/2\tau _0)\mathrm{ln}E/\mathrm{ln}5`$ in a static system with the same gluon density as the initial value of the expanding system at $`\tau _0`$. With $`R_A6`$ fm and assuming $`\tau _00.2`$ fm, this would give $`(dE/dL)_07.3`$ GeV/fm for a 10-GeV parton, Since the parton energy loss is directly proportional to gluon density of the medium, this implies that the gluon density in the initial stage of $`Au+Au`$ collisions at $`\tau _0=0.2`$ fm/$`c`$ is about 15 times higher than that inside a cold nucleus. To extract the parton density and the initial formation time, one has to measure the energy loss and $`q_T`$ broadening simultaneously.

In Fig. 2, we also show the predicted suppression factor at large $`p_T`$ with the same logarithmic energy-dependence of the energy loss as in the cold nuclear matter which gives a suppression factor that increases with $`p_T`$. However, new RHIC data  at 200 GeV show a constant suppression factor at larger $`p_T`$. This indicates the importance of the detailed balance  which gives much stronger energy dependence of the energy loss. We show as the lower dashed line in Fig. 2 the result for an effective parton energy loss

$$\mathrm{\Delta }E(E/\mu 1.6)^{1.20}/(7.5+E/\mu )$$
(9)
which is parameterized according to the result from Ref.  where both stimulated gluon emission and thermal absorption are included in the calculation of the total energy loss. The detailed balance between emission and absorption reduces the effective parton energy loss and on the other hand increases the energy dependence. The threshold is the consequence of gluon absorption that competes with radiation that effectively shuts off the energy loss. The parameter $`\mu `$ is set to be 1 GeV in the calculation shown in Fig. 2.

This work was supported by DOE under Contract No. DE-AC03-76SF00098.


text:
# Electromagnetic Meissner effect in spin-one color superconductors

## Abstract

It is shown that color-superconducting quark matter, where quarks of the same flavor form Cooper pairs with spin one, exhibits an electromagnetic Meissner effect. This is in contrast to spin-zero color superconductors where Cooper pairs consist of quarks with different flavors.

Quantum chromodynamics (QCD) is an asymptotically free theory asymp  and, thus, quark matter at large quark chemical potential $`\mu `$ is a weakly coupled system. In this case, the dominant interaction between two quarks is single-gluon exchange. Single-gluon exchange is attractive in the color-antitriplet channel. Consequently, at sufficiently low temperatures, the quark Fermi surface is unstable with respect to the formation of Cooper pairs. Since this is analogous to what happens in ordinary superconductors bcs , this phenomenon was termed color superconductivity bailin .

Color superconductivity was studied from first principles in the framework of QCD at weak coupling QCDgapeq ; rischke1  as well as in more phenomenological Nambu–Jona-Lasinio-type models alford1 . Both approaches indicate that the color-superconducting state is the true ground state of quark matter at any density beyond the quark-hadron phase transition and at sufficiently low temperature. They also agree in the magnitude of the color-superconducting gap parameter, $`\varphi `$, which they predict to be of the order of tens of MeV (for quark Cooper pairs with total spin zero) at densities of the order of ten times the nuclear matter ground state density. Gap parameters of this order of magnitude may have enormous phenomenological implications, since the transition temperature to the normal conducting phase is typically of the order of $`\varphi `$ rischke1 ; wang ; schmitt . For instance, during the evolution of a neutron star, the temperature ranges from a few tens of MeV down to a few keV pons . If its core is sufficiently dense to consist of quark matter, this quark matter core is then very likely a color superconductor.

Color superconductivity is more complex than ordinary superconductivity, because quarks do not only carry electric, but also color and flavor charge. Two different quark flavors may form Cooper pairs in the color-antitriplet, flavor-singlet, total spin-zero channel (the so-called 2SC phase) bailin . For three different quark flavors, the favored state is the so-called color-flavor locked (CFL) phase alford2 , with spin-zero Cooper pairs in the color-antitriplet, flavor-antitriplet representation.

A necessary condition for pairing of quarks of different flavor is that their respective Fermi surfaces are equal. In physical systems, however, this condition may be hard to achieve. For instance, compact stellar objects are neutral with respect to electric and color charge. This requires the introduction of separate chemical potentials for quarks which differ in color and flavor neutrality . Imposing the constraints of color and electric charge neutrality leads to different values for these chemical potentials. This effect and the mass difference between the light up and down quarks and the heavier strange quark berges  may then lead to different Fermi surfaces for each quark species. If this difference is of the order of, or larger than, the color-superconducting gap parameter $`\varphi `$, a color-superconducting state where quarks with different flavor form spin-zero pairs with zero momentum is no longer favored. Then, besides a transition to the normal conducting state bedaque , there are at least four other possibilities. The first two possibilities are either a displacement bowers  or a deformation muether  of the Fermi spheres of the two quark species forming the Cooper pair, breaking translational or rotational invariance, respectively. The third possibility is an interior gap structure for the quark species with the larger Fermi momentum liu .

The fourth possibility is that quarks of the same flavor form Cooper pairs with total spin one bailin ; iwasaki ; rischke1 ; schaefer ; alford3 . Quarks with the same flavor have the same mass and the same electric charge, and thus pairing is neither affected by a mass difference nor a nonzero electric chemical potential which may be required to fulfill the constraint of electric neutrality. Moreover, a potentially nonzero color chemical potential does not destroy the Cooper pairs, because we expect it to be parametrically much smaller than the gap, $`\mu _{\mathrm{color}}\varphi ^2/(g\mu )\varphi `$, where $`g`$ is the strong coupling constant colorneutrality . (A color chemical potential is necessary in models which treat color as a global charge. In QCD, the role of the color chemical potential is assumed by a constant gluon background field $`A_a^0`$ colorneutrality .)

According to the results of Refs. schmitt ; schaefer , the gap in spin-one color superconductors is of the order of 20 - 400 keV (assuming that the gap in the 2SC phase is of the order of 10 - 100 MeV). The critical temperature is therefore of the order of 10 - 400 keV schmitt . Consequently, when the core of a neutron star cools below this temperature, it could very well consist of quark matter in a spin-one color-superconducting state.

The interesting question is how such a state affects the properties of the star and whether this leads to observable consequences. The best known properties of a neutron star are its radius and mass, which are determined by the equation of state. A recent study alfordreddy  shows that radii and masses of compact stellar objects with a color-superconducting quark matter core do not change appreciably from the values expected for ordinary neutron stars. Another property of a neutron star with potentially observable consequences is its magnetic field. Due to an admixture of protons in neutron matter, the core of an ordinary neutron star is a superconductor, and magnetic fields experience the Meissner effect. However, if the core of a neutron star is a spin-zero color superconductor (for instance in the 2SC or CFL phase), there is no electromagnetic Meissner effect emprops . Since charge neutrality may favor a spin-one over a spin-zero color-superconducting state, a natural question is whether the electromagnetic Meissner effect is also absent in a spin-one color superconductor. In order to answer this question, in this letter we study the pattern of symmetry breaking of the local gauge symmetries, and then also present the results of an explicit calculation of the Meissner masses.

In ordinary superconductors, the electromagnetic gauge group $`U(1)_{em}`$ is broken due to the fact that the electrons in a Cooper pair carry electric charge. This leads to the electromagnetic Meissner effect, i.e., magnetic fields only penetrate a finite distance into the superconductor. The inverse distance can be associated with a nonzero photon mass, the so-called Meissner mass. Since quarks are not only electrically charged but also carry color, besides the electromagnetic $`U(1)_{em}`$ symmetry also the $`SU(3)_c`$ gauge symmetry of the strong interaction is broken in a color superconductor. This leads to the color Meissner effect, i.e., the gluons obtain Meissner masses and color-magnetic fields are expelled. The question is whether all eight gluons and the photon become massive. This depends on the particular pattern of how the local symmetries are broken in the color superconductor. If there is a residual local symmetry, the corresponding gauge bosons remain massless. This residual symmetry group of $`SU(3)_c\times U(1)_{em}`$ leaves the gap matrix $`\mathrm{\Delta }`$ invariant,

$$(g_c\times g_{em})\mathrm{\Delta }(g_c^T\times g_{em}^T)\stackrel{!}{=}\mathrm{\Delta },$$
(1)
where $`g_cSU(3)_c`$, $`g_{em}U(1)_{em}`$, and $`T`$ denotes the transpose. In general, the gap matrix $`\mathrm{\Delta }`$ is a matrix in color, flavor, and Dirac (spin) space bailin ; QCDgapeq ; rischke1 ; alford1 . Since pairing occurs in the attractive color-antitriplet channel, the color structure of the gap matrix corresponds to the color-antitriplet $`\overline{\mathrm{𝟑}}_c`$ representation of the $`SU(3)_c`$ gauge group. In a spin-one color superconductor, the spin structure of the gap matrix corresponds to the symmetric spin triplet $`\mathrm{𝟑}_J`$ representation of the $`SU(2)_J`$ spin group, which is also a representation of $`SO(3)_J`$. The gap matrix is diagonal in flavor space, since quarks in a Cooper pair carry the same flavor. For the moment, let us consider quark matter with a single flavor only, $`N_f=1`$. The case of several quark flavors (where each flavor pairs at its respective Fermi surface) will be discussed below. The gap matrix can be written as

$$\mathrm{\Delta }=\mathrm{\Phi }_a^iJ_av^i,$$
(2)
where $`J_a`$ and $`v_i`$ $`(a,i=1,2,3)`$ are bases of $`\overline{\mathrm{𝟑}}_c`$ and $`\mathrm{𝟑}_J`$, respectively, and $`\mathrm{\Phi }_a^i`$ is the order parameter. The form of the order parameter defines the phase of the condensate. As in helium-3, there is a multitude of possible phases for spin-one condensates vollhardt . In this letter we only consider the polar phase and the color-spin locked (CSL) phase schmitt ; schaefer . The order parameters in these phases are

$$\left(\mathrm{\Phi }_a^i\right)_{polar}\delta _{a3}\delta ^{i3},\left(\mathrm{\Phi }_a^i\right)_{CSL}\delta _a^i.$$
(3)
In the polar phase, the condensate points in a fixed direction in real space, which breaks the global spatial symmetry group $`SO(3)_J`$ to $`SO(2)_J`$. The condensate also points in a fixed direction in color space, which spontaneously breaks the local $`SU(3)_c`$ gauge symmetry to a residual $`SU(2)_c`$ gauge group. From Eq. (1) we deduce that the residual subgroup which leaves the order parameter invariant is generated by the three generators of $`SU(2)_c`$ (corresponding to $`T_1,T_2,T_3`$ of the original $`SU(3)_c`$) and the generator $`\stackrel{~}{Q}_{polar}=Q2\sqrt{3}qT_8`$, where $`T_8`$ is one of the generators of $`SU(3)_c`$ and $`Qq\mathrm{𝟏}_J`$ is the generator of $`U(1)_{em}`$. The constant $`q`$ is the electric charge of the single quark flavor considered here (2/3 for $`u`$ quarks and $`1/3`$ for $`d`$ or $`s`$ quarks). The generator $`Q`$ is proportional to the unit matrix in spin triplet space, $`\mathrm{𝟏}_J`$, since all states of the spin triplet have the same electric charge. Consequently, the symmetry breaking pattern is $`SU(3)_c\times U(1)_{em}SU(2)_c\times \stackrel{~}{U}(1)`$, where $`\stackrel{~}{U}(1)`$ is generated by $`\stackrel{~}{Q}_{polar}`$. The existence of a nontrivial residual gauge symmetry is equivalent to the fact that there are charges with respect to which the Cooper pairs are neutral. These are the two color charges corresponding to the $`SU(2)_c`$ gauge symmetry and the $`\stackrel{~}{Q}_{polar}`$ charge corresponding to the $`\stackrel{~}{U}(1)`$ gauge symmetry. The gauge boson of the latter is a superposition of the photon and the eighth gluon of $`SU(3)_c`$. This superposition is mathematically given by an orthogonal rotation of the original gauge fields by an angle $`\theta `$. In general, a generator $`\stackrel{~}{Q}=Q+\eta T_8`$ results in a mixing angle given by $`\mathrm{cos}^2\theta =g^2/(g^2+\eta ^2e^2)`$, where $`e`$ is the electromagnetic coupling constant emprops . In our case the mixing angle, $`\theta _{polar}`$, is determined by this expression with $`\eta =2\sqrt{3}q`$. Since $`eg`$, the mixing angle is small, $`\theta _{polar}2\sqrt{3}qe/gq/3`$. Consequently, the main contribution to the gauge boson of the local $`\stackrel{~}{U}(1)`$ symmetry comes from the original photon, with a small admixture of the eighth gluon. This justifies to call this gauge boson the “new” photon. There is no electromagnetic Meissner effect, since the new photon can penetrate the color-superconducting phase. This is similar to other color-superconducting phases, for instance the 2SC phase or the CFL phase emprops . In both cases, there is a $`\stackrel{~}{U}(1)`$ gauge symmetry and thus no electromagnetic Meissner effect.

In the CSL phase, the order parameter breaks $`SU(3)_c\times SO(3)_J`$ to the diagonal subgroup $`SO(3)_{c+J}`$ schaefer . This is analogous to the breaking of color-flavor symmetries in the CFL phase, where $`SU(3)_c\times SU(3)_fSU(3)_{c+f}`$. The residual $`SO(3)_{c+J}`$ and $`SU(3)_{c+f}`$ groups are global symmetries, and thus not gauged. This similarity between CSL and CFL phases does, however, not extend to the behavior concerning electromagnetism. Unlike the CFL phase, it turns out that for the CSL phase there is no nontrivial subgroup of $`SU(3)_c\times U(1)_{em}`$ that leaves the gap matrix invariant, i.e., the only possible solution to Eq. (1) is $`g_c=g_{em}=\mathrm{𝟏}`$. Consequently, the symmetry breaking pattern is $`SU(3)_c\times U(1)_{em}\mathrm{𝟏}`$. This is equivalent to the fact that a Cooper pair in the CSL phase is not neutral, neither with respect to ordinary electric charge nor with respect to any possible new combination of color and electric charge. This fact has the physical consequence that there is an electromagnetic Meissner effect for the CSL phase schaefer .

In the following, we confirm the above qualitative arguments by an explicit calculation of the Meissner masses. The Meissner mass is defined as the zero-energy, zero-momentum limit of the spatial $`(ii)`$ components of the polarization tensor meissner

$$m_{ab}^2\text{lim}_{p0}\mathrm{\Pi }_{ab}^{ii}(0,p),a,b=1,\mathrm{},8,9,$$
(4)
where the first eight indices correspond to the gluons, and the ninth index to the photon, $`9\gamma `$. If there is mixing of gluons and photons, the $`9\times 9`$ Meissner mass matrix $`m_{ab}^2`$ is not diagonal. The Meissner masses for the physically relevant modes are obtained by diagonalizing this matrix. In weak coupling, the polarization tensor may be computed in one-loop approximation. At zero temperature, only the quark loop contributes. Our calculation follows the method of Ref. meissner ; details are deferred to a future paper future . It turns out that the gluon part of the mass matrix is diagonal, $`m_{ab}^2\delta _{ab}m_{aa}^2`$ for $`a,b=1,\mathrm{},8`$. In Table 1 we collect all results.

For the polar phase, the vanishing masses for gluons 1, 2, and 3 indicate the unbroken $`SU(2)_c`$ subgroup. The nonzero Meissner mass $`m_{8\gamma }^2`$ reflects the mixing of eighth gluon field, $`A_8`$, and the photon, $`A_9A_\gamma `$. In terms of the physically relevant modes, $`\stackrel{~}{A}_a`$, the Meissner mass matrix is diagonal,

$$\underset{ab}{}A_am_{ab}^2A_b\underset{a}{}\stackrel{~}{A}_a\stackrel{~}{m}_a^2\stackrel{~}{A}_a.$$
(5)
For $`a=1,\mathrm{},7`$, $`\stackrel{~}{A}_aA_a`$ (and, correspondingly, $`\stackrel{~}{m}_a^2m_{aa}^2`$), whereas the new gluon, $`\stackrel{~}{A}_8`$, and the new photon, $`\stackrel{~}{A}_\gamma `$, are obtained by an orthogonal rotation of $`A_8`$ and $`A_\gamma `$,

$$\left(\begin{array}{c}\stackrel{~}{A}_8\\ \stackrel{~}{A}_\gamma \end{array}\right)=\left(\begin{array}{cc}\mathrm{cos}\vartheta & \mathrm{sin}\vartheta \\ \mathrm{sin}\vartheta & \mathrm{cos}\vartheta \end{array}\right)\left(\begin{array}{c}A_8\\ A_\gamma \end{array}\right),$$
(6)
where $`\mathrm{tan}(2\vartheta )2m_{8\gamma }^2/(m_{88}^2m_{\gamma \gamma }^2)`$. With the numbers of Table 1 we find that the rotation angle $`\vartheta `$ is identical to the mixing angle $`\theta _{polar}`$ found in the above group-theoretical argument. The Meissner masses for the rotated gauge fields are

$$\stackrel{~}{m}_8^2=\left(\frac{1}{3}g^2+4q^2e^2\right)\frac{\mu ^2}{6\pi ^2},\stackrel{~}{m}_\gamma ^2=0.$$
(7)
The massless new photon confirms that there is no electromagnetic Meissner effect in the polar phase, whereas there is a color Meissner effect for the new gluon.

For the CSL phase, the particular pattern for the gluon masses reflects the residual global $`SO(3)_{c+J}`$ symmetry which is generated by a combination of the antisymmetric Gell-Mann matrices $`T_2`$, $`T_5`$, and $`T_7`$ and the generators of $`SO(3)_J`$. In contrast to the polar phase, now all mixed masses are zero, indicating that there is no mixing between gluons and the photon. All gauge bosons have a nonzero Meissner mass. This coincides with the above group-theoretical argument.

Let us now discuss the situation with more than one quark flavor, where the quarks of each flavor separately form spin-one Cooper pairs. In general, each quark flavor has a different chemical potential, $`\mu _i`$, $`i=1,\mathrm{},N_f`$. In this case, the results of Table 1 have to be modified. All gluon Meissner masses $`m_{aa}^2`$ have to be multiplied by a factor $`_i(\mu _i/\mu )^2`$. In the mixed masses $`m_{a\gamma }^2`$ we have to replace the single quark electric charge $`q`$ by a factor $`_iq_i(\mu _i/\mu )^2`$. Finally, one has to substitute the square of the charge $`q^2`$ in the photon masses $`m_{\gamma \gamma }^2`$ by a factor $`_i(q_i\mu _i/\mu )^2`$.

For the polar phase, these modifications have the following effect. Only if $`_{i,j}q_i(q_iq_j)(\mu _i\mu _j)^2=0`$, there is a massless combination of photon and gluon. For a system with two quark flavors, for instance $`u`$ and $`d`$, or $`u`$ and $`s`$, this condition is equivalent to $`\mu _u^2\mu _d^2=0`$, or $`\mu _u^2\mu _s^2=0`$, respectively. This means that for these two combinations of flavors, there is always a Meissner effect, unless the chemical potential of one of the two quark flavors is zero, i.e., this flavor is absent. For a two-flavor system consisting of $`d`$ and $`s`$ quarks, however, the above condition is always trivially fulfilled, i.e., there is no Meissner effect in this case. This is to be expected, since $`d`$ and $`s`$ quarks carry the same electric charge and thus appear as one single flavor with respect to electromagnetic interactions. For three flavors, the above condition is equivalent to $`\mu _u^2(\mu _d^2+\mu _s^2)=0`$, and can only be fulfilled if either $`\mu _u=0`$, or $`\mu _d=\mu _s=0`$. Both cases have been discussed above. We therefore conclude that, for more than one quark flavor, there is an electromagnetic Meissner effect in the polar phase, unless all quarks carry the same electric charge.

In the CSL phase, the mixed masses always vanish, while the results for the gluon masses and the photon mass are modified by the same factors as for the polar phase. Thus, the conclusion that there is no photon-gluon mixing and an electromagnetic Meissner effect in the CSL phase remains valid also for more than one quark flavor.

While spin-zero color superconductors could be of type II at small $`\mu `$ blaschke ; baym , a spin-one color superconductor is most likely always of type I, because the ratio of the penetration depth to the coherence length is of order $`\varphi /(e\mu )10^3[100\mathrm{MeV}/(e\mu )]1`$. Consequently, the magnetic field is completely expelled from the core of a compact stellar object, if it is a spin-one color-superconductor. This is true unless the magnetic field exceeds the critical field strength for the transition to the normal conducting state. The magnetic field in neutron stars is typically of the order of $`10^{12}`$ Gauss. This is much smaller than the critical magnetic field which, from the results of Ref. baym , we estimate to be of the order of $`10^{16}`$ Gauss.

In conclusion, a compact stellar object with a core consisting of quark matter in the spin-one color-superconducting state is, with respect to its electromagnetic properties, different from an ordinary neutron star: ordinary neutron star matter is commonly believed to be an electromagnetic superconductor of type II, while a spin-one color superconductor is of type I. The question whether neutron star matter is a type-I or type-II superconductor has recently stirred a lot of attention zhit , because it was shown that type-II superconducting matter is incompatible with the observation of pulsars with precession periods of order 1 year link . The presence of spin-one color-superconducting quark matter in the pulsars’ core could explain this observation.

Acknowledgment. We would like to thank M. Alford, M. Lang, K. Rajagopal, T. Schäfer, and I.A. Shovkovy for comments and discussions. Q.W. acknowledges financial support from the Alexander von Humboldt-Foundation.


text:
# Are Scale Free Networks Better?

## I Introduction

Random graphs were studied since the middle of the 20th century with the initiative works of mathematicians Pál Erdős and Alfréd Rényi. In 1960 they published their paper On the evolution of random graphs, with the first thorough study of the graph theory er  . However, these results were strictly theoretical, since tools to measure real-world random graphs at that time were not available. With the evolution of the personal computers, we have now the possibility to study this topic in practice, and, interestingly, the computers themselves provide one of the most exciting real random graphs: the computer networks.

As these networks evolved, their properties were analyzed, and soon it was clear, that the Erdős-Rényi models (ER) were not appropriate for these graphs www . The most important difference was the *scale free* property, that the computer networks have, but the classical models lack. It turned out, that this feature is shared by other types of networks, like social social  and metabolic meta  ones. It means that the distribution of the degrees of the nodes follows a power-law distribution, while the classic graphs has a Poisson degree distribution, with an exponential tail. It shows that, in the real-world graphs, nodes with a high number of connections are much more likely than expected. Therefore new models were necessary to describe these kind of random graphs.

One of these models was developed by László Barabási and Réka Albert (BA model) alb ; abr . Their method has two key features: *incremental growth* and *preferential attachment*. Incremental growth means that the graph is constructed by adding nodes to the existing graph, and connecting them according to a construction rule, contrary to the original Erdős-Rényi picture with a static graph. Preferential attachment means the the likelihood of a connection depends on the degree of a node, again lifting the classical assumption of equal probabilities. This model has more variations, and they are able to describe a wide class of random graphs. There is another class of models study such networks without incremental growing, purely on their statistical properties zck .

In this paper, we investigated network properties of BA-like models with different parameters and raise the question whether non-classical graphs may perform better transmitting data over them. Certain value of parameters yield a graph with an exponential tail degree distribution, hence allowing us to compare these classical type models to the ones with power-law distribution. In Section II we analyze the properties of the extended BA models, estimating their node distribution. Next, in Section III we construct networks based on these models, and simulate a network traffic on such a graph in a simplified model estimating the traffic load on the nodes. The scaling of the load with nodal degree is presented, allowing to estimate the proper bandwidth allocation when upgrading a node. In Section IV we study a more realistic setup with TCP dynamics and compare the theoretical result of the previous section to the simulated ones, while in Section V we discuss the overall performance of the different simulations. Finally, we conclude out analysis.

## II The model

In the original BA-model the graph is constructed as follows. Starting from a small initial graph, we extend it by adding a new node in each step and connecting it to $`m`$ randomly selected existing nodes. The probability of choosing a particular node is proportional to its degree,

$`p_i={\displaystyle \frac{d_i}{_{j=1}^{n1}d_i}}.`$ (1)
This model constructs a scale free graph, where the cumulative degree distribution has a power-law decay with an exponent 2 alb .

In the following, we modify the construction law, and use the more general form

$`p_i={\displaystyle \frac{d_i^\alpha }{_{j=1}^{n1}d_i^\alpha }},`$ (2)
weighting the probability with a power $`\alpha `$ of the nodal degree, similarly to Ref. krapivsky . The BA model corresponds to $`\alpha `$=1, while for $`\alpha =0`$ the preferential connectivity is cancelled, and we are back to a classical ER-like graph model with a uniform distribution, leading to exponential node distribution. While this model is similar to the original ER model, it differs in some aspects, such as it has a minimal guaranteed degree, and the ordering of the nodes presents nonzero correlations in the degrees callaway . With $`\alpha (0,1)`$, the models provide a smooth transition between the classical and the scale free models.

The degree distribution for these models can be derived following the method introduced in alb . Here we give a fast estimate on the distribution, for the exact result see krapivsky . First, we estimate the rate of growth of the degree at each node, assuming that the growth of the degree is continuous in time. At time $`t`$ there are exactly $`t`$ nodes and $`mt`$ links between them. Hence, the expected degree value $`k_i=E(d_i)`$ of node $`i`$ is growing as

$`_tk_i(t)={\displaystyle \frac{mk_i(t)^\alpha }{_{j=0}^tk_j(t)^\alpha }}.`$ (3)
For $`\alpha =0`$, the denominator is simply conunting the number of nodes, and is equal to $`t`$, while for $`\alpha =1`$, it is (double) counting the number of links, and thus is $`2mt`$, both being a linear function of the time $`t`$. Numerical simulations showed that for $`\alpha (0,1)`$ the denominator is well approximated by a linear function $`ct`$, where $`c`$ is the $`\alpha `$ dependent measured slope. Hence, generally our differential equation can be written as

$`_tk_i(t)={\displaystyle \frac{k_i(t)^\alpha }{ct}}.`$ (4)
Fortunately, these equations are solvable for each $`\alpha `$ in the choosen $`(0,1)`$ range. Specifically, for $`\alpha =0`$ we get

$`k_i(t)={\displaystyle \frac{\mathrm{log}t}{c}},`$ (5)
and more generally for $`\alpha (0,1)`$

$`k_i(t)=\left((1\alpha ){\displaystyle \frac{\mathrm{log}t}{c}}\right)^{\frac{1}{1\alpha }}.`$ (6)
Finally, for the scale free case $`\alpha =1`$ the structure changes, and we arrive at

$`k_i(t)=t^{\frac{1}{c}}.`$ (7)
In order to get a formula for the degree distribution, note that for each $`i`$, $`k_i(i)=m`$, since at the time the node is added it has exactly $`m`$ connections. The growth process for node $`i`$ has a similarity to the growth of the previous nodes, and this suggests that

$`k_i(t)=k_1\left({\displaystyle \frac{t}{i}}\right).`$ (8)
The cumulative distribution $`P(k_i>x)`$ counts the number of nodes with $`k_i(t)>x`$, i.e. $`k_1\left(\frac{t}{i}\right)>x`$, which leads to

$`{\displaystyle \frac{t}{k_1^1(x)}}>i.`$ (9)
Normalizing back with the total number of nodes at time $`t`$ (being also $`t`$) we get the approximate cumulative probability

$`P(k_i>x){\displaystyle \frac{1}{k_1^1(x)}}.`$ (10)
Inverting Eq. (6) we find

$`P(k_i>x)e^{\frac{c}{1\alpha }x^{1\alpha }},`$ (11)
the Weibull-distribution for $`\alpha (0,1)`$ with scale parameter $`c`$, and shape parameter $`1\alpha `$. We note, that the exact result is of the form krapivsky

$`P(k_i>x)x^\alpha e^{\frac{c}{1\alpha }x^{1\alpha }\mathrm{}},`$ (12)
where $`\mathrm{}`$ denotes higher order correction terms.

For $`\alpha =0`$ ($`c=1`$) one recovers the classical limit $`P(x)e^x`$. In the region $`\alpha (0,1)`$ the distribution still vanishes exponentially, so strictly speaking no heavy tails are present, however, the probability of finding nodes with large degree increases dramatically. Numerical studies revealed, that the distribution indeed follows the form (11), however, the parameters are slightly different from the predicted ones, due to the approximation used in (4). In the limit $`\alpha =1`$ ($`c=2`$) we arrive at the genuine scale free heavy tail result

$`P(x)x^c.`$ (13)
We note, that the analysis can be extended to $`\alpha >1`$, however, such a distribution actually tends to the degenerate case with one node being connected with all others in a star-like topology. For $`\alpha =1`$ the original BA model results in $`c=2`$, however, different growth stategies are able to vary the value of $`c`$, hence the decay exponent. In this work we restrict ourselves to the original BA model.

## III Estimation of the network load in a simplified model

In the following, first we generate random graphs on a computer according to the extended BA model (3) and study the load in a simplified network model. We assume, that each node generates a constant traffic (data stream) to all the other nodes. The amount of traffic is the same for all connections, and furthermore, we assume that the link capacity can handle the accumulated traffic. The properties of the network are analysed tuning the geometrical parameter $`\alpha `$.

For each value of $`\alpha `$ we performed measurements on 8 different graphs (generated with the same statistical properties), and then averaged. Each topology was made out of $`N=100000`$ nodes and $`m=3`$ links per node. The traffic was routed using the standard shortest-path strategy, and the number of data streams, passing a node was counted for each node. Since by construction at least $`m(N1)`$ data channel are open per node, we define the load $`l_i`$, to be the number of data channels divided by $`N`$, the number of nodes, at node $`i`$. In the large $`N`$ limit this gives the number of links $`m`$, if no “foreign” connection is going through the selected node. One can regard this quantity as “weighted degree”, with weights equal to the link loads.

To study the effect of congestion, we assume, that for $`N`$ nodes a data stream occupies $`1/N`$ part of the link capacity. With $`N(N1)`$ connections there will be certainly “overused” links, where the link bandwidth capacity constrains the througput of the node. Thus we define the throughput of a link to be the minimum of the load and the link capacity, while the throughput of the node to be the sum of the link troughputs over the links connected to the given node.

Figure 1 shows the tails of the cumulative distribution functions ($`1.0CDF`$) for degrees, load and (node) throughput on graphs with $`\alpha `$=0, 0.5 and 1. The distribution of the *degrees* show the expected tails, exponential ($`\alpha =0`$, ER graph), Weibull ($`\alpha =0.5`$), and the power-law ($`\alpha =1`$, BA graph), respectively. However, the distribution of load shows interesting deviations from the one of the nodal degree. For $`\alpha =0`$ the weighted nodal degree (load) distribution can be approximated much better by a Weibull distribution than an exponential one, describing the nodal degree distribution. This is a clear indication, that the distribution of the network traffic even in a simulated “uniform” situaton does not follow one-by-one the underlying network topology, rather developes a heavier tail. For $`\alpha =0.5`$ the load remains Weibull, however, the shape parameter (theoretically $`1\alpha `$), changes from the numerical value $`0.62`$ to $`0.32`$. In case of $`\alpha =1`$ the power-law decay survives, but with changing CDF tail exponent, descreasing from the numerical value $`1.97`$ to $`1.25`$. We may conclude, that the load distribution has considerably fatter tails than the underlying nodal distribution.

The per-node *throughput* shows a transition between the nodal degree and the load distribution. For low throughput values it follows the load distribution (the bandwidth is enough to hold the traffic), however, at higher values it approaches the nodal degree distribution, simply counting the number of links. The transition is governed by the link capacities, in a real network environment with large available capacity one expects the througput to follow the load distribution, however, in an underdesigned network it follows the degree distribution, as we show in the next section.

The change of the distribution from exponential to Weibull, and the parameter changes in the Weibull and in the power-law case in the load distribution suggests to examine the correlation between the degree and the load. This is shown in Figure 2 (left) on a scatter-plot. Since it is log-log scaled, the linear clusters indicate a power-law correlations, $`l_i=d_i^\beta `$. The dependence of the scaling exponent $`\beta `$, on the network geometry parameter $`\alpha `$, is shown in Figure 2 (right). The load is pushed to have fatter tails than the degree distribution, and the more “classical” is the network the larger the deviation. For the BA geometry the the load distribution decays $``$ 1.6 times faster than the corresponding degree distribution.

This means that not only the distribution of the load is more heavy-tailed than the distribution of the degrees, but also means that this dependence is quite strictly a power-law function. It also explains the deviations from the degree distributions, since a power of an exponentially distributed ($`\alpha =0`$) random variable is Weibull, while in the case of Weibull or power-law degree distributon the transformation results only in a parameter change for the distribution. It seems that the traffic pushes the distribution into heavier tails, from the exponential distribution, and converging to a power-law through Weibulls. The initial push at $`\alpha =0`$ is extremly high, with an exponent 1/2 in the Weibull distribution. The load distribution is much less sensitive to the underlying network as the degree distribution.

Since the nodal throughput is bounded from above by the degree, it limits the throughput on the high-degree nodes, where the links are already fully utilised. Therefore the throughput-degree joint distribution is different from the load-degree distribution, the power-law correlation is only valid for low and medium-degree nodes. But in real life the congestion at the overloaded nodes also affects the other parts of the network, since every data flow through these nodes is jammed. Furthermore, the TCP dynamics is also known to be chaotic pisti  which may change the scalings observed in a simpler model. To simulate the real life situation, and compare them with the results of the simplified model, next, we simulate a realistic traffic on a computer network, too.

## IV Estimation of the network load in a simulated traffic

In order to compare the theoretical results from the previous section to a more realistic setup, we simulated a TCP/IP network with $`N=1024`$ nodes, $`m=3`$ links per node, and a uniform link bandwidth of $`1`$ Mb/s. Every node communicated with a randomly selected target node. To study the effect of congestion we modelled three scenarios: a low, a medium, and a high traffic one, using constant bitrate data flows of $`16`$ Kb/s, $`64`$ Kb/s, and $`256`$ Kb/s, respectively. The simulation was ran using the Berkeley Network Simulator package ns . The link throughput was calculated as the number of packets sent through that link, and the nodal values as the sum of the throughput of the incident links.

To our surprise, in each scenario the throughput showed a good scaling through the whole degree range. Since the throughput at the high-degree nodes is obviously limited, it must mean that the congestion limits the throughput of the other nodes in such way that the power-law throughput-degree correlation remains valid. It also means that the congestion affects each other node proportionally to its degree.

The correlation exponent $`\beta `$ was, however, different in the three simulation. In the low-traffic scenario the measured distribution is exactly the same as obtained from the numerical simulations of the load in the previous section. The congestion still has not set up, and the throughput is identical to the load. As the traffic intensity grows, exponent $`\beta `$ decreases, and flattens. It shows that the heavy traffic is more evenly spreads through the network, but its dependence on the degrees remains a power-law function.

## V Overall performance

The next quantity we studied is the total number of transferred packets in the simulated network. Each node transmits constantly TCP packects to its randomly choosen partner, and if a packet arrives, an acknowledgement is sent back to the originating node. We counted the number of packets for which the acknowledgement was received and the difference betweeen the number of transmitted packets and acknowledged packets is the loss.

The source of the loss of packets in the network is the congestion: whether the link capacity cannot handle the amount of traffic, or the nodes in between cannot cope with the routing of the packages. One would expect, that in a network with smaller shortest pathes between two randomly choosen nodes the load on the links and routers is higher, hence using the same devices it would drop the packets more often as a network with larger shortest path.

It is also known, that scale free network have smaller shortest pathes connecting two arbitrary nodes abr ; cluster , i.e. a scale free network uses less routing devices, however, the load on them is higher. Indeed, numerical simulation, performed in the previous section also showed, that with increasing traffic the performance of the scale free (BA) network downgrades, for example, with a drop rate of 24% already at medium traffic, while the classical (ER) network show a downgrade only of 6% for the same traffic.

## VI Conclusion

In this paper we showed, that the distribution of the network traffic even in a simulated “uniform” situaton does not follow one-by-one the underlying network, rather it developes fatter tails than the nodal degree distribution of the network. A scaling between the nodal degree distribution and the load on a node was established, showing a power like patter, $`ld^{\beta (\alpha )}`$, where the scaling exponent $`\beta `$, is a decreasing function of the network parameter, $`\alpha `$. For networks with fatter nodal distribution the exponent is smaller. As a consequence, attaching a new connection to the node requires less upgrade in the bandwidth for a scale free network to keep the performance, as for a classical (ER) one.

The above theoretical result was confirmed by simulation, the per-node throughput is still scaling with nodal degree. The scaling dependes on the amount of the traffic, for a completely congested situation the throughput distribution by definition agrees with the nodal degree distribution, hence the scaling exponent is 1, however, for partly congested or congestion free networks this scaling approaches the load distribution with scaling exponents in the range 1.4 (partially congested BA network) to 2 (congestion free ER). A scale free network requires less resource upgrade when a new node is added.

The overall performance of a scale free network is decreasing rapidly with the traffic, where the classical network still has almost no loss. It is due to the feature, that the ER network uses more routers along the shortest connection, and hence the traffic is distributed more evenly.


text:
# High-redshift clusters from NVSS: The TexOx Cluster (TOC) Survey

## 1 Introduction

Galaxy clusters have traditionally been used as cosmological probes bf98 , and despite the much heralded recent dawning of the age of precision cosmology, they will remain useful as such in determining the normalisation of the matter power spectrum, the cosmic equation of state haiman01 , and refining measurements of other cosmological parameters. The current concordance of cosmological constraints, however, now makes it feasible to rely on the redundancy between measurement methods tz02  to use the predictions of the standard $`\mathrm{\Lambda }`$CDM model to test cluster and AGN physics rather than the reverse.

Clusters, especially at high redshift, are dynamic environments, rife with mergers bh96 , galaxy harassment m98 , enhanced star formation z96 , tidal stripping go72  and other processes which transform morphology and stellar populations. Many of these processes hinder the identification of clusters, as all methods make implicit assumptions about what a cluster should look like, and then try to find objects fulfilling such criteria, making each method subject to (often ill-understood) selection effects.

Whilst radio selection of clusters does not offer a “clean” method by which to construct an unbiased cluster sample, the selection effects are at least different, and the lack of dust obscuration in the radio enables one to probe high redshifts more easily v02 . Radiogalaxies tend to be located in clusters hl91 ; best00  and links between cluster formation and radiogalaxy triggering sr02  are expected in the dynamic environments noted above. This makes the comparison of radio-selected clusters with those detected by other methods an attractive means by which AGN triggering, fuelling, evolution, duty cycle, obscuration and environment may be studied.

## 2 The TexOx Cluster (TOC) Survey

The cosmic evolution of the radiosource population w01  suggests that some high-redshift clusters will contain several radio-loud AGN. By targetting imaging observations towards fields containing overdensities of radiosources, we ought to be able to detect such clusters. The TexOx (Texas–Oxford) Cluster Survey utilises this technique by searching for $`7\sigma `$ over-densities at $``$ mJy levels in $`7^{}\times 7^{}`$ boxes within the NVSS radio catalogue NVSS , finding 331 within an $`1128\mathrm{deg}^2`$ sky area TOC .

The TOC algorithm allows each source to be a member of no more than one overdensity. This prevents large overdensities from being detected several times, artificially increasing the number of “real” clusters detected. Running the TOC algorithm multiple times on random simulations of the same sky area (with the same source density) as one of the regions of the real TOC survey resulted in a Gaussian distribution of number of clusters detected. The mean of this Gaussian, 65.5, was very significantly ($`6\sigma `$) less than the number of clusters, 103, found in the same region of the NVSS.

We conclude, therefore, that at most $`60\pm 7\%`$ of TOC clusters are chance superpositions of unrelated radiosources. However, many of these fields are expected to contain at least one real cluster due to the high probability that even single radiosources may be found in clusters (and even if all sources are at different redshifts, TOC is $`5`$ times more efficient than looking around isolated sources). Inspection of the NVSS radio maps suggests that $`20\%`$ of TOC fields are extended structures resolved by NVSS into several components, leaving $`20\%`$ of fields to contain at least two (and potentially many more) physically-related radiosources. This means that we would expect TOC to discover $`65`$ fields containing two or more radiosources at the same redshift.

By targetting those fields where none of the NVSS sources are identified on the POSS-II POSS  images, we are able to preferentially target high-redshift clusters. Around $`30\%`$ of NVSS sources have POSS-II counterparts s02 , so assuming we are able to remove the extended structures mentioned above from the sample, there is only a $`17\%`$ chance that a given field of 5 radiosources would have no POSS-II identifications if the radiosources were randomly distributed in redshift. This POSS-II constraint (POSS-II can see $`L_{}`$ galaxies at $`z0.3`$), the $`z1`$ median redshift of NVSS sources, and the fact that the TOC $`7^{}\times 7^{}`$ box size is comparable to the size of a cluster core at $`z0.2`$, means that TOC is highly efficient at finding high-redshift clusters.

Deep optical imaging (30 – 60 minutes per band in $`R`$ and $`I`$ at the McDonald 2.7-m in around 130 fields) was used to select around 20 promising high-redshift candidates for near-IR followup in $`J`$ and $`H`$ at the Calar Alto 3.5-m. Additionally some $`K`$ data was obtained at UKIRT, and deeper, higher-resolution radio maps with the VLA.

## 3 TOC J0233.3+3021 - a high redshift cluster?

One particularly promising high-redshift cluster candidate was TOC J0233.3+3021, a field in which the radiosources were identified with galaxies of $`K18`$, $`J20`$, $`R24`$. Our collaborators in Cambridge observed this field with the Ryle Telescope and measured a 15 GHz SZ decrement sz  of $`695\mu `$Jy, at $`7\sigma `$ significance, consistent with the presence of a cluster with a minimum gas mass of $`5\times 10^{13}M_{}`$ cotter02 .

The magnitudes of the radiogalaxies strongly suggest the presence of a cluster at $`z1`$; attempts to obtain spectroscopic redshifts of the NVSS identifications using ISIS on the WHT resulted in the detection of faint red continuum in the four brightest cases, but no evidence of emission lines. We used the publicly available software hyperz hyperz  to obtain photometric redshifts using our $`BRIJHK`$ photometry and determined best-fit redshifts of $`z1.4`$ in four cases, and $`z1.0`$ for the remaining radiogalaxy. The best solutions were $`36L_{}`$ elliptical templates in the former cases, and a $`2L_{}`$ S0 in the latter. One of the ellipticals has a poorly constrained redshift and is consistent with being at $`z1.0`$.

The consistency of the photometric redshifts suggests that at least three of the radiogalaxies may be members of a structure at $`z1.4`$; this is supported by their position on the radiogalaxy $`Kz`$ diagram (Fig. 1). From the integrated 1.4 GHz radio fluxes, we obtain rest-frame 1.4 GHz luminosities111Assuming a spectral index $`\alpha =0.8`$, where $`L_\nu \nu ^\alpha `$ of $`10^{25}`$ W Hz-1 sr-1, close to the FR-I/II break lo96 . One of the ellipticals hosts a radiosource which is clearly an FR-II from our 8 GHz VLA maps (with a projected size of 140 kpc); the rest are compact doubles or are unresolved. Two of the sources appear to have distorted FR-I type morphologies, as might be expected if they were in clusters, or perhaps groups falling into a cluster.

Although it seems probable that the SZ signal is associated with the inferred overdensity of radiogalaxies at $`z1.4`$, it is of course possible that it is associated with a cluster at a different redshift. In an attempt to detect the signature of the presence of a cluster we are substantially aided by the availability of deep radio, optical and IR data, but hampered by the lack of spectroscopic redshifts. In an attempt to identify a cluster red sequence, we plotted a colour-magnitude diagram for those objects in the field with both $`R`$ and $`H`$ detections and magnitude errors $`<0.3`$. The resulting 297 objects are shown in Fig. 3. There is perhaps a hint of two red sequences; the redder contains 3 of the radiosources and is similar in colour and magnitude to that seen by Rosati et al. rosati99  in a $`z=1.26`$ cluster. The lower redshift red sequence contains the other two radiosources (consistent with the photometric redshift analysis).

Number counts of objects from our Gemini GMOS image of this field show a slight excess over galaxy number counts obtained by other authors postman98  at $`i^{}20`$, possibly associated with the inferred $`z1.0`$ structure, but the presence of a cluster at $`z1.4`$ is almost impossible to ascertain from these data due to the faintness of the objects, the large surface densities of background galaxies, and a drop in completeness in our data at these magnitudes. Radio number counts, however, show an excess of sources, both at the several mJy level at which the overdensity was selected (hardly surprising!) and, very tentatively, at fainter flux densities, before the contrast becomes lost in the rising background source counts (Fig 3).

## 4 Summary

The TOC method provides an efficient way to find clusters over a wide range of redshift, out to $`z>1`$. The contrast provided against the background radiosource counts at mJy levels makes the selection relatively easy, but distinguishing the cluster population from background and foreground sources will require a good deal of effort, ideally using multi-object spectroscopy. We note that the difficulties of clearly identifying clusters at high-redshift, especially in the “redshift desert” ($`1.2z2`$) compounds the problem; indeed other authors haynes02  have also had difficulties identifying the systems responsible for SZ decrements in some fields.

Our technique suffers from the problem that it tends to find superpositions of systems at different redshifts, complicating matters further. This may prove to be an advantage, however, should one wish to study the population of a more distant cluster using a foreground cluster as a gravitational lens.

### Acknowledgments

We thank N. Drory, our Calar Alto collaborator; G. Cotter, H. Buttery, R. Das, K. Grainge, W. Grainger, M. Jones, G. Pooley, and R. Saunders, our Cambridge SZ collaborators; and J. Tufts and P. Gay for help with observing. The work of SC was performed with the support of a PPARC PhD studentship, and under the auspices of the U.S. Department of Energy, National Nuclear Security Administration by the University of California, Lawrence Livermore National Laboratory under contract No. W-7405-Eng-48. This material is based in part upon work supported by the Texas Advanced Research Program under Grant No. 009658-0710-1999.


text:
# A TIDALLY-DISRUPTED ASTEROID AROUND THE WHITE DWARF G29-38

## 1 INTRODUCTION

The formation and evolution of large rocky solids such as asteroids and terrestrial planets are of great interest. Even though ever-increasing numbers of extra Solar System Jovian-mass gas-giant planets are known (see Marcy & Butler 1998), no rocky planet orbiting a main sequence star outside of our own Solar System has been directly discovered. There is good indirect evidence for an asteroid belt around the nearby (distance = 22 pc) main sequence A-type star, $`\zeta `$ Lep (Chen & Jura 2001) and there are rocky planets around the pulsar PSR 1257+12 (Wolszczan & Frail 1992). Here, we propose that the infrared excess around the white dwarf G29-38 is the consequence of the tidal disruption of an asteroid.

Most white dwarfs have infrared colors characteristic of their photospheres, but G29-38 was discovered by Zuckerman & Becklin (1987) to have an infrared excess which can be explained by circumstellar dust at a temperature of approximately 1000 K. The origin of these warm grains which must be close to the white dwarf is an unsolved mystery (Zuckerman 2001).

Graham et al. (1990) proposed that the circumstellar dust around G29-38 arose from an unspecified catastrophic destruction of an asteroid or comet near the white dwarf. Here, we present a model where the dust around G29-38 is so close to the star that it lies inside the Roche radius where an asteroid parent-body would have been tidally disrupted, a process suggested by Debes & Sigurdsson (2002) to destroy comets and that would produce the observed circumstellar dust. Our model, where a flat, opaque dust ring ultimately evolves from a tidally-disrupted asteroid, may solve another mystery: the presence of calcium and other metals in the atmosphere of G29-38 (see Koester, Provencal & Shipman 1997, Zuckerman & Reid 1998) where the calculated dwell time for this material is only $``$ 20 years (Paquette et al. 1986). Alcock, Fristrom & Siegelman (1986) suggested that the atmospheric composition of a white dwarf could be detectably altered by the direct impact of a comet. In the specific case of G29-38, we suggest a somewhat similar scenario where enough accretion to pollute the atmosphere arises from the tidal disruption of an asteroid. The model that accretion comes from the disruption of asteroids or comets naturally explains the absence of much hydrogen in the infalling material – a major difficulty with the commonly suggested hypothesis (Dupuis, Fontaine & Wesemael 1993) that white dwarf atmospheres are contaminated by accretion of interstellar matter.

## 2 MODEL FOR THE DUST DISTRIBUTION AND EMISSION

In order to model the circumstellar dust emission, we need to characterize G 29-38, a ZZ Cet variable (see Kleinman et al. 1998). We adopt a distance of 14.1 pc (Tokunaga, Becklin & Zuckerman 1990), a stellar radius, $`R_{}`$, of 8.2 $`\times `$ 108 cm, an effective temperature, $`T_{}`$, of 11820 K, and a mass, $`M_{}`$, of 0.69 M (Bergeron et al. 1995). These parameters yield a mean density, $`\rho _{}`$, of 6.0 $`\times `$ 106 g cm-3, a stellar luminosity of 2.4 $`\times `$ 10-3 L and imply a cooling age of 4 $`\times `$ 108 yr (Winget et al. 1987).

We assume that this star illuminates a passive, opaque, flat circumstellar ring where the incident optical energy is re-radiated in the infrared. When variations in the vertical ring structure are ignored, the standard solution for the temperature, $`T_{ring}`$, as a function of distance from the star, $`R`$, is given by the expression (Chiang & Goldreich 1997):

$$T_{ring}\left(\frac{2}{3\pi }\right)^{1/4}\left(\frac{R_{}}{R}\right)^{3/4}T_{}$$
(1)
The predicted temperature of the ring from equation (1) is lower than that expected for unshielded, isolated black bodies, and the very innermost grains are likely to be warmer than predicted by equation (1). As with Saturn’s rings, the vertical thickness of the ring around G29-38 may not be much greater than the size of the particles (Nicholson 2000), which typically is less than 10 cm (Graham et al. 1990). Thus the amount of unshielded material compared to the total mass in the ring may be tiny.

The predicted flux from the ring, $`F_{ring}`$ is:

$$F_{ring}=\frac{2\pi cosi}{D_{}^2}_{R_{in}}^{R_{out}}B_\nu (T_{ring})R𝑑R$$
(2)
where $`D_{}`$ is the distance to the star from the Sun, $`R_{in}`$ and $`R_{out}`$ represent the inner and outer radii of the dust ring, $`i`$ is the inclination angle of the ring and $`B_\nu `$ is the Planck function. With the substitution that $`x`$ = $`h\nu /kT_{ring}`$, then:

$$F_{ring}\mathrm{\hspace{0.33em}12}\pi ^{1/3}\left(\frac{R_{}^2cosi}{D_{}^2}\right)\left(\frac{2k_BT_{}}{3h\nu }\right)^{8/3}\left(\frac{h\nu ^3}{c^2}\right)_{x_{in}}^{x_{out}}\frac{x^{5/3}}{e^x\mathrm{\hspace{0.33em}1}}𝑑x$$
(3)
The comparison of the predictions from equation (3) with published observations for the ring emission is shown in Figure 1. To fit the near infrared data, we assume that the innermost temperature of the ring is 1200 K, appropriate for the sublimation of refractories. From equation (1) and the assumed stellar parameters, this implies that $`R_{in}`$ = 1.0 $`\times `$ 1010 cm or 0.14 R. Also, because the ring is relatively bright compared to the stellar photosphere, we assume that $`i`$ = 0. Finally, in view of the scatter in the data, a range of values of $`R_{out}`$ fit the data about equally well with 0.4 $`R_{}`$ $`<`$ $`R_{out}`$ $`<`$ 0.9 $`R_{}`$, near the value of the tidal radius described below and corresponding to an outer temperature between 600 K and 300 K. The agreement between the model and the data shown in Figure 1 is notable, especially since the only parameter in equation (3) that is unconstrained by the physical model is the ring inclination.

Since the ring is opaque, its total mass cannot be determined from the infrared excess alone. A minimum mass can be derived from the assumption that the ring is optically thick even at 15 $`\mu `$m, the longest observed wavelength. The total surface area of the ring, $`A`$, is $``$1022 cm2. Although the size distribution and composition of the particles in the ring are unknown, the maximum likely grain opacity, $`\chi `$, at 15 $`\mu `$m, is 1000 cm2 g-1 (Ossenkopf, Henning & Mathis 1992). In this case, the minimum ring mass, given by $`A/\chi `$, is 1018 g. If the particles are 10 cm in diameter, then $`\chi `$ $``$ 0.05 cm2 g-1, and the minimum mass is 2 $`\times `$ 1023 g. For reference, Saturn’s rings have a total mass of 2 $`\times `$ 1022 g and an inner radius of 0.13 R (Esposito 1993).

## 3 RING ORIGIN AND EVOLUTION

We now consider the origin and evolution of the ring. Following Debes & Sigurdsson (2002), we picture that an asteroid’s initial relatively circular orbit was drastically altered after the Asymptotic Giant Branch evolution of the progenitor to G29-38 with the consequence that it strayed within the Roche radius of the star and was tidally disrupted. The debris from this disruption-event evolved into the opaque ring that we infer to be present, a model which is similar to that proposed by Dones (1991) for the formation of Saturn’s rings.

While a complete description of the tidal interaction between an asteroid and a white dwarf is complex, in the simplest picture, an asteroid of density, $`\rho _a`$ near a star of density, $`\rho _{}`$ and radius $`R_{}`$ is disrupted within a region at separation, $`R_{tide}`$, given by the expression (Davidsson 1999):

$$R_{tide}=C_{tide}\left(\frac{\rho _{}}{\rho _a}\right)^{1/3}R_{}$$
(4)
where $`C_{tide}`$ is a numerical constant of order unity which depends upon the orbital parameters of the asteroid, its rotation and composition. For an asteroid with $`\rho _a`$ = 3 g cm-3 around G29-38, $`R_{tide}`$ $``$ 1.5 $`C_{tide}`$ R, and therefore, the model dust ring described above lies within the region where an asteroid would be tidally destroyed. Once the asteroid is disrupted, the debris of dust particles and rocks undergoes a complex evolution. Mutual collisions dissipate energy but not angular momentum so the system evolves towards a flat ring composed of matter in circular orbits, reminiscent of Saturn’s rings (Dones 1991). Another important process is the inward drift of material towards the central star under the action of the Poynting Robertson effect. The innermost radius of the dust ring is established by the location where the particles sublimate.

We can indirectly infer the minimum mass of the tidally-disrupted asteroid if the ring material ultimately accretes onto the star. To account for the observed photospheric metallicity, Koester et al. (1997) estimated an infall rate of interstellar matter onto G29-38 of 5 10-15 M yr-1, a rate which which is probably not strongly affected by the star’s pulsations (see Corsico et al. 2002). Since interstellar matter is hydrogen-rich, the metal-infall rate is $``$10-16 M yr-1. If accretion occurs from a rocky ring for 20 years, the total mass of refractory material in the ring must be at least 4 1018 g, implying a radius greater than 10 km. The total mass of accreted matter onto the star is probably much greater than this value since it is likely that G29-28 has had calcium in its photosphere for much longer than 20 years. Since there are metals in the atmospheres of over 10% of all white dwarfs (Zuckerman & Reid 1998), it is plausible that G29-38 has been actively accreting for 10% of its cooling age or 4$`\times `$ 107 years. In this case, the total mass of accreted material is $``$4 $`\times `$ 1024 g, comparable to the 2 $`\times `$ 1024 g in the asteroid belt of the Solar System (Binzel et al. 2000) or, equivalently, about 1% of the mass of the asteroids around $`\zeta `$ Lep (Chen & Jura 2001). While G29-38 may have been bombarded by many asteroids with radii $``$ 10 km, it is also possible that the ring has arisen from the disruption of a single object with a radius of 1000 km.

Why did the asteroid stray near the white dwarf? One possibility is that the parent-body to the dust ring was initially in a nearly circular orbit where it was always substantially further from G29-38 than 1 AU so it could survive the second red giant branch evolution of the star. Additionally, there was at least one Jovian-mass planet orbiting the main sequence progenitor. As described by Debes & Sigurdsson (2002), when the Asymptotic Giant Branch progenitor to G29-38 lost over half of its initial main sequence mass, a previously stable configuration of planetary orbits could become unstable. In this new environment, an asteroid’s orbit would evolve dramatically and it could acquire a high eccentricity by interactions with the Jovian-mass planet(s) (Ford, Havlickova & Rasio 2001). The ultimate fate of an asteroid in its new orbit would be (1) to collide with a perturbing planet or another asteroid, (2) to collide with the central star (or, as we propose above, be tidally disrupted by the white dwarf) or (3) to be ejected from the system (see Tremaine 1993). The relative probability of each of these fates depends upon the particular set of planet and asteroid masses and orbits; according to Debes & Sigurdsson (2002), there are a number of configurations where the asteroid ultimately collides with the central star. The duration of this unstable phase might last over 108 years since it is sensitive to the initial conditions of the system.

## 4 RELEVANCE TO OTHER WHITE DWARFS

Both helium-rich (Dupuis et al. 1993) and hydrogen-rich (Zuckerman & Reid 1998) white dwarfs sometimes show metals in their atmospheres, yet only G29-38 is known to have an infrared excess. Do the stars with metals also happen to possess dusty rings which have not yet been detected? Consider the following evidence. (1) Most of the hydrogen-rich white dwarfs with metals identified by Zuckerman & Reid (1998) have effective temperatures less than 9,000 K. It can be shown from equation (3) that at 2.2 $`\mu `$m, the expected ring emission is only 10% of the photospheric flux and may not be easily recognized. Since most infrared surveys of white dwarfs have been restricted to the J, H and K bands, observations at longer wavelengths may be required to detect the rings (if they exist) around these white dwarfs. (2) Dupuis et al. (1993) summarize the detections of metals in helium-rich white dwarfs including 12 stars with $`T_{eff}`$ $`>`$ 10,000 K. If these systems had rings similar to G29-38’s, they would have infrared excesses greater than 20% of the photospheric emission at 2.2 $`\mu `$m, unless the ring inclination was nearly edge-on. Most of these stars are so faint (mK $`>`$ 15 mag) that the 2MASS data are not precise enough to determine with confidence whether there is a 20% excess. (See also Wachter et al. 2003 who describe a search for late-type companions to white dwarfs.) However, one of these 12 stars, G111-54, has been observed with sufficient precision and the stellar parameters are well enough known (Bergeron et al. 2001) that if $`i`$ = 0, the estimated upper limit to the ring flux at 2.2$`\mu `$m is only about 20% of that predicted by a ring model similar to G29-38’s. If, however, the innermost grains have a maximum temperature of 1000 K instead of 1200 K, then there could be an opaque ring around the star since the cooler grains would not emit very much at 2.2 $`\mu `$m. (3) Chary, Becklin & Zuckerman (1999) observed 12 mainly hydrogen-rich white dwarfs with ISO at 15 $`\mu `$m, and according to equation (3), the ring emission should dominate the total flux. Except for G29-38, the sources either were undetected or the emission is dominated by the photosphere of the star with the typical upper limits to the ring emission being 0.4 mJy. Three of the stars (WD 0208+396, WD 0322-019 and WD 1633+433) are cool enough and also far enough from the Sun that the predicted fluxes from ring emission at 15 $`\mu `$m for $`i`$ = 0 are $``$1 mJy, and, within the errors and uncertainties, consistent with the available upper limits. However, for G1-27 (van Maanen’s star), the predicted flux from the ring for $`i`$ = 0 is near 10 mJy, depending upon the exact inner and outer radii of the ring, which is well in excess of the observed upper limit. Finally, there is no detected infrared excess around the hydrogen-rich white dwarf G238-044 which has an effective temperature of 20,000 K, so the dwell time for metals is only a few days (Holberg, Barstow & Green 1997). With a stellar radius of 8.4 $`\times `$ 108 cm (Provencal et al. 1998) and a distance from the Sun of 25 pc (Vauclair et al. 1997), the predicted flux from the ring at 15 $`\mu `$m for $`i`$ = 0 is $``$5 mJy, depending upon the inner and outer ring radii, while the upper limit given by Chary et al. (1999) is 0.2 mJy. In summary, there are probably some white dwarfs with atmospheric metals that do not currently possess opaque, dust rings similar to G29-38’s.

Planetary rings in the Solar System are diverse; rings around white dwarfs may be as well. There are several possible explanations for the lack of detected infrared excesses around most white dwarfs even if they accrete asteroids. (1) Portions of Saturn’s rings and those of other planets are optically thin (Esposito 1993). If particles in a ring are far enough apart and the system is on average optically thin, equation (3) overestimates the infrared flux. (2) If the tidally-disrupted asteroid or comet was composed largely of volatile material, then at the Roche radius, solids would sublimate. In this case, the ring may be mainly gaseous and emit little infrared radiation. (3) The accretion may be time dependent. Asteroids might stray into the Roche region episodically, and the ensuing dust ring may be dissipated in a time short compared to the gravitational settling time out of the atmosphere of the white dwarf which in some cases can be over 105 years.

G29-38 may have an unusually large infrared excess because its disrupted asteroid was unusually massive. For example, in the Solar System, the most massive asteroid, Ceres, carries about half of the total mass in the asteroid belt (see Binzel et al. 2000; Michalak 2000). There is such a large range in asteroid masses that the rings that result from their destruction may vary substantially.

## 5 CONCLUSIONS

1. We have found that the infrared excess around the white dwarf G29-38 can be explained as an opaque ring of refractory dust with inner and outer radii of 0.14 R and $``$1 R, respectively.

2. We propose that the circumstellar ring, which lies within the Roche radius of the white dwarf, was produced by a tidally-disrupted asteroid.

3. Accretion from the dust ring can plausibly explain the abundance of calcium in the atmosphere of G29-38. Either as a bombardment by a series of asteroids or because of one large disruption, the total amount of matter accreted onto the white dwarf may have been $``$ 4 $`\times `$ 1024 g, about 1% of the mass of asteroids around $`\zeta `$ Lep.

I have had very useful conversations with E. Becklin, B. Hansen, and B. Zuckerman. I thank NASA for support.


text:
# A new transiting extrasolar giant planet

## Acknowledgements

We wish to thank A. Udalski and the OGLE team for generous contributions to this project. We are very grateful to S. Kulkarni for invaluable support, to R. Noyes and D. Latham for helpful comments, to T. Barlow for assistance with the spectroscopic reductions using MAKEE\[Barlow¡22¿\], and to K. Stanek for his continuous encouragement. The data presented herein were obtained at the W. M. Keck Observatory (operated by Caltech, Univ. of California, and NASA), which was made possible by the generous financial support of the W. M. Keck Foundation. M.K. gratefully acknowledges the support of NASA through the Michelson Fellowship programme. G.T. acknowledges support from NASA’s Kepler Mission. S.J. thanks the Miller Institute for Basic Research in Science at UC Berkeley for support via a research fellowship.

Correspondence and requests for materials should be addressed to M.K. (e-mail: maciej@gps.caltech.edu)