text:
# ISOCAMBased on observations with ISO, an ESA project with instruments funded by ESA Member States (especially the PI countries: France, Germany, the Netherlands and the United Kingdom) with the participation of ISAS and NASA. observations of Galactic Globular Clusters: mass loss along the Red Giant Branch.

## 1 Introduction

A complete, quantitative understanding of the physics of mass loss processes and the precise knowledge of the gas and dust content in Galactic globular clusters is crucial in the study of Population II stellar systems and their impact on the Galaxy evolution. It also has major astrophysical implications on related problems such as the ultraviolet excess found in elliptical galaxies (Greggio & Renzini 1990; Dorman et al. 1995) and the interaction between the intracluster medium and the hot halo gas (cf. e.g. Faulkner & Smith, 1991).

Indirect evidence for mass loss in Population II stars includes the observed morphology of the horizontal branch (HB) in the cluster color–magnitude diagrams, the pulsational properties of the RR Lyrae stars, and the absence of asymptotic giant branch (AGB) stars brighter than the red giant branch (RGB) tip (Rood 1973, Fusi Pecci & Renzini 1975; 1976, Renzini 1977, Fusi Pecci et al. 1993, D’Cruz et al. 1996). The expected mass loss is about $`0.2M_{}`$ along the RGB and about $`0.1M_{}`$ along the AGB (e.g. Fusi Pecci & Renzini 1976). As a consequence of such mass loss processes, dust and gas should be present in the intracluster medium, as diffuse clouds (cf. e.g. Angeletti et al. 1982) or concentrated in circumstellar envelopes. If no cleaning mechanism is at work between two Galactic plane crossings, a few tens solar masses of intracluster matter should be accumulated in the central regions of the most massive clusters (i.e. those with large central escape velocity).

Early searches (see Roberts 1988, Smith et al. 1990 and references therein) for this intracluster gas (HI, HII, CO) and more recently detection of ionized gas in 47 Tuc (Freire et al. 2001), yielded upper limits or marginal detections well below 1 $`M_{}`$.

Mid–IR excesses and scattered polarized light have also been observed in the central region of massive globular clusters (Frogel & Elias, 1988; Gillet et al. 1988; Forte & Mendez 1989; Minniti et al., 1992; Origlia, Ferraro & Fusi Pecci 1995, Origlia et al. 1997b). They are mainly associated with long period variables evolving along the AGB. These preliminary results seem to indicate that in globular clusters even the brightest IR point sources are fainter than 1 Jy in the 10 $`\mu `$m spectral region and most of the emission comes from warm circumstellar dust with a minor photospheric contribution. Any diffuse component is much fainter than 10 mJy$`\mathrm{arcsec}^2`$.

Other evidence of dust in the intracluster medium of globular clusters comes from the presence of IRAS sources (Lynch & Rossano 1990, Knapp & Gunn 1995, Origlia et al. 1996) in their central regions and from more recent surveys using ISOPHOT (Hopwood et al. 1999) and ground based sub-millimeter antennas (cf. e.g. Origlia et al. 1997a, Penny, Evans & Odenkirchen 1997, Hopwood et al. 1998).

However, all these studies show that globular clusters are deficient in diffuse intracluster matter and that some mechanism(s) must be at work to remove the gas and dust which should have accumulated between each passage through the Galactic plane (cf. e.g. Faulkner & Smith, 1991).

Given this observational scenario the desirability of a deep mid-IR survey of the central regions of globular clusters in the continuum and in selected dusty features is obvious. Such a deep, spatially resolved photometric survey became possible with the spectrophotometric capabilities of ISOCAM on board of the Infrared Space Observatory (ISO, Kessler et al. 1996). ISOCAM (Cesarsky et al. 1996) provided relatively fine spatial resolution, large field coverage, and high sensitivity in the 10 and 20 $`\mu `$m spectral regions.

Ramdani & Jorissen (2001) performed a deep survey at 12 $`\mu `$m in the external regions of 47 Tuc at different distances from the cluster center, with the specific goal of studying the mass loss during the AGB evolutionary stage. They cross-correlated their ISOCAM survey with the DENIS survey and found dust excess only in two well known bright variables, confirming the strong link between stellar pulsation activity and mass loss modulation.

With the goal of studying the mass loss during the RGB evolutionary stage, a deep survey of the very central regions of six, massive globular clusters has been made using ISOCAM with two different filters in the 10 $`\mu `$m spectral region. Mid-IR observations have the advantage of sampling an outflowing gas fairly far from the star (typically, tens/few hundreds stellar radii). At 10 $`\mu `$m the sampled material typically left the star a few decades previously. Thus mass loss rates inferred from the mid-IR will smooth the mass loss rate over a few decades. For astrophysical purposes it is the longterm average mass loss which is important. Conceptually one could sample different distances and smoothing times by observing at different wavelengths, however at this time that is not practical. At longer wavelengths detectors lack the requisite sensitivity; at shorter wavelengths the photosphere dominates.

## 2 The data set

Five massive Galactic globular clusters, NGC 104 (=47 Tuc), NGC 362, NGC 5139 (=$`\omega `$ Cen), NGC 6388, NGC 7078 (=M15) and NGC 6715 (=M54) in the Sagittarius Dwarf Spheroidal were observed with ISOCAM between February and August 1997, (proposals ISO\_GGCS and DUST\_GGC, P.I.: F. Fusi Pecci). Each of the clusters has an IRAS point source in its core (Origlia, Ferraro & Fusi Pecci 1996). The 12 $`\mu `$m IRAS flux is 1.7 Jy in 47 Tuc and 1.2 Jy in NGC 6388 and ranges between 0.4–0.5 Jy in the others. 47 Tuc, NGC 6388, and NGC 5139 show extended IRAS emission at 12 $`\mu `$m as well, with typical sizes of 2–4′ and total flux densities between 3 and 11 Jy.

Beam-switching observations were performed using CAM03 and 6″$`\mathrm{pixel}^1`$ scale in two filters, LW7\[9.6\]111For clarity we give the central wavelength (in $`\mu `$m) of the ISOCAM filters in brackets. (8.3–10.9 $`\mu `$m, centered on a silicate dust feature) and LW10 (8.6–15.4 $`\mu `$m, the IRAS 12 $`\mu `$m band). The total integration time for each cluster was about 2500s and the field of view was about 3′$`\times `$3′. 47 Tuc was also observed a second time in four filters, LW6\[7.7\] (6.9–8.5 $`\mu `$m), LW7\[9.6\], LW8\[11.3\] (10.7-11.9 $`\mu `$m) and LW10, for a total integration time of about 6700s and a field coverage of about 5′$`\times `$3′.

Near-IR photometry of the clusters observed with ISOCAM was obtained at ESO, La Silla (Chile), using the ESO-MPI 2.2m telescope and the near IR camera IRAC-2 (Moorwood et al. 1992) equipped with a NICMOS-3 256$`\times `$256 array detector, during different runs between 1992 and 1994. The frames were taken through standard $`J`$ and $`K`$ broad band filters. The overall spatial coverage is about 4′$`\times `$4′ with a 0$`\stackrel{}{\mathrm{.}}`$5 pixel. More details on the reduction procedure and the photometric calibration are given by Ferraro et al. (1994; 2000) and Montegriffo et al. (1995). $`\omega `$ Cen was also observed (May 1999) with the ESO-NTT telescope using SOFI, the near-IR imager/spectrometer equipped with a 1024$`\times `$1024 array detector. The images were acquired through the $`J`$ and $`K_s`$ filters using a 0$`\stackrel{}{\mathrm{.}}`$3 pixel scale. The $`K_s`$ photometry has been converted into the standard $`K`$ system using the stars in common with the IRAC-2 survey, and the calibration of Origlia, Ferraro & Fusi Pecci (1995). Accurate Point Spread Function (PSF) fitting reductions of the stellar fields have been carried out using the ROMAFOT package (Buonanno et al. 1983).

Fig. 1 shows the $`K`$ band images of the six clusters with the ISOCAM 12 $`\mu `$m contour plots superimposed. The latter were obtained using the background subtracted, mosaiced images (CMOS) generated by the ISOCAM auto-analysis pipeline (7.0 version of the off-line processing system).

## 3 Stellar Point Sources

The ISOCAM point source catalogs were also generated using the auto-analysis pipeline. For each beam–switching pointing the pipeline generates a calibrated flux image (CMAP) and the corresponding catalog of detected point sources. The point source searching and photometry was performed by maximum-likelihood estimates and PSF fitting. The nominal flux uncertainty is $`20`$%.

For each filter configuration we cross-correlated the various ISOCAM point source catalogs from the different beam-switching pointings and generated a final catalog with average fluxes and positions. This procedure has also allowed us to remove spurious detections due to possible transient or memory effects. The fluxes were first color–corrected, according to a blackbody energy distribution in the range of temperatures between 3000 and 5000 K, typical of the globular cluster red giants. The applied color–correction factors for the LW6\[7.7\], LW7\[9.6\], LW8\[11.3\] and LW10 filters are 1.06, 1.01, 1.00 and 1.26, respectively.

In order to convert the observed fluxes into magnitudes, 86.5, 66.4, 43.6 and 28.3 Jy were adopted as zero–magnitude fluxes at the reference wavelengths of the LW6\[7.7\], LW7\[9.6\], LW8\[11.3\] and LW10 filters, respectively. These zero point values were obtained using the IRAS flux density calibration at 12 $`\mu `$m (cf. Beichman et al. 1985) and a 10,000 K blackbody to obtain the zero-magnitude fluxes at the other reference wavelengths.

The near-IR stellar counterpart of each ISOCAM point source was identified by roto-translating the ISOCAM coordinate system into the near-IR system and cross-correlating the two catalogs with an overall accuracy $``$2 arcsec. Given the large pixel size of the ISOCAM images, several stars are often associated to each mid-IR source on the basis of the spatial coordinates alone. The brightest one within a square box of 30″ $`\times `$ 30″ (5 $`\times `$ 5 ISOCAM pixels) was selected. In all cases the brightest star turns to be also the reddest in (J-K).

The near-IR $`J`$ and $`K`$ magnitudes were de-reddened and transformed into bolometric magnitudes using the reddening and distance moduli given in Table 1, and the bolometric corrections as a function of the de-reddened $`(JK)_0`$ colors listed by Montegriffo et al. (1998). The distance moduli in Table 1 are those from Harris (1996) with a $`+0.2`$ mag correction, based on the new calibration by Ferraro et al. (1999).

The average accuracy of the relative photometry is $`0.2`$ mag in the mid-IR and $`0.1`$ mag in the near-IR. The uncertainty of the zero point calibrations for absolute photometry is somewhat larger, a conservative estimate being $`0.5`$ mag.

In the two closest clusters, 47 Tuc and $`\omega `$ Cen, the limiting bolometric magnitude is $`M_{\mathrm{bol}}1.0`$; in M15 and NGC 6388 (the most distant ones in our Galactic sample) the detection threshold is $`M_{\mathrm{bol}}2.5`$. Our statistical analysis is based on this second limit. M54 in the Sagittarius Dwarf Spheroidal is more distant and even the brightest stars are close to the detection limit of our ISOCAM survey. The number of stellar counterparts identified in each cluster is given in Table 2. The $`K`$ band stellar maps of the six clusters down to a bolometric magnitude $`M_{\mathrm{bol}}=2.5`$ and their associated ISOCAM point sources are plotted in Fig. 2. The only known variables in this sample are V42 in $`\omega `$ Cen (a long period variable) and six other stars with small amplitude variability and unknown period present in the observed field of view of 47 Tuc (cf. Montegriffo et al. 1995 and reference therein).

Using the ISOCAM mosaiced images (CMOS, cf. Fig. 1) we also performed aperture photometry of the central core region of each cluster within an $`1\stackrel{}{\mathrm{.}}5`$ radius. We then compared these integrated flux densities with those obtained by summing the contributions of each ISOCAM point source over the same field of view and determined an upper limit to the percentage concentrated in circumstellar (CS) envelopes with identified stellar counterparts. Both the total fluxes and these percentages are given in Table 2. These percentages scale approximately with the geometrical dilution, and very similar values are also inferred using the $`K`$ band flux, which is fully dominated by the stellar photosphere. Such a behavior suggests that the 12 $`\mu `$m emission not associated with bright circumstellar envelopes is most likely due to the unresolved stellar background, although some contribution by diffuse warm dust heated by the cluster radiation field can be also possible.

## 4 Color excess and dust parameters

Suitable color–magnitude and color–color diagrams of the ISOCAM point sources were constructed in order to quantify the relative photospheric and dust excess contribution. Fig. 3 shows the $`M_{\mathrm{bol}},(JK)_0`$ and $`M_{\mathrm{bol}},(K[12])_0`$ color–magnitude diagrams. Fig. 4 shows the corresponding $`(K[12])_0,(K[9.6])_0`$ color–color diagram.

The $`(JK)_0`$ color mainly traces the photospheric temperature, thus defining the RGB ridge line of each cluster (Ferraro et al. 2000). Possible circumstellar dust excess is best shown using a combination of near and mid-IR colors. Overall a good correlation between the $`(K[12])_0`$ and $`(K[9.6])_0`$ colors is found, indicating that both the broad band continuum filter LW10 and the narrow band filter LW7\[9.6\] centered on a silicate dust feature are good tracers of dust excess. However, since the broad LW10 filter is somewhat more sensitive in detecting mid-IR excess and taking into account the actual color scatter along the observed RGB, $`(K[12])_00.65`$ was adopted as a conservative, border value between pure photospheric emission and significant dust excess. Mass loss at a lower rate could be traced by bluer $`(K[12])_0`$ colors as well, but the observed photometric scatter and the spatial resolution of our mid-IR survey is not adequate to further constrain mass loss in this regime.

20 out of the 52 ISOCAM sources detected in the upper 1.5 bolometric magnitudes of the RGB show direct evidence of mid-IR circumstellar dust excess: 2 variables (V8 and V26, cf. Montegriffo et al. 1995) and 5 normal giants in 47 Tuc, 6 in NGC 6388, 3 in NGC 362, 2 in M15, the only candidate cluster member in M54 and the long period variable V42 in $`\omega `$ Cen (also see Table 2). The de-reddened photometry of these sources with dust excess is reported in Table 3. From pure evolutionary timescale considerations (e.g., Renzini & Fusi Pecci 1988), 10% or less of bright giants are expected to be in the AGB phase. Since we observed the very central region of each cluster where most of the luminous giants should be located, we expect less than 5 AGB stars in our global sample, hence when lacking further information we neglect them since they have a minor impact on the overall statistics. One exception is V42 in $`\omega `$ Cen: this star well fits the Period-Luminosity relation of classical Mira variables (especially if our slightly larger distance modulus is adopted, cf. Feast, Whitelock & Menzies, 2001) and it should be more likely on the AGB evolutionary stage, so we do not consider it in our further statistical analysis on the mass loss along the RGB.

In NGC 6388 and in M15, the most distant clusters in our Galactic sample, also very centrally concentrated, both the total number of ISOCAM sources and the number of those with mid-IR excess must be regarded as lower limits, since stellar counterparts to the 12 $`\mu `$m emission cannot be identified in the very central core where the crowding is too severe given the ISOCAM spatial resolution.

Somewhat surprisingly, no significant mid-IR excess has been found around any giant star in the observed central region of $`\omega `$ Cen. However, this cluster is known to be anomalous in many respects, with a multi-population RGB (cf. Pancino et al. 2000). More crucially in the current context, while the cluster is large in absolute size, it is very low density, and our sample size is small. Taking advantage of the wide field photometry recently published by Pancino et al. (2000), it has been possible to estimate that only $`<30`$% of the brightest giants are located in the central region covered by our mid-IR survey, while in the other clusters we sample almost all of the brightest giants. In order to obtain a similar statistical significance as in the other clusters of our survey, a much larger mid-IR spatial mapping is definitely needed.

In order to estimate the dust parameters a simple model using an optically and geometrically thin shell at constant temperature and with a $`\nu B_\nu `$ grain emissivity (cf. e.g. Natta & Panagia 1976; Origlia, Ferraro & Fusi Pecci 1996; Origlia et al. 1997b) has been computed. Standard grain radius and density of 0.1 $`\mu `$m and 3 g cm-3, respectively, are adopted.

The mid-IR fluxes of those sources showing dust excess have been re-calibrated according to this modeling, using an average dust temperature of 300 K. The applied color–correction factors for the LW6\[7.7\], LW7\[9.6\], LW8\[11.3\] and LW10 filters are 1.01, 1.01, 1.00 and 0.95, respectively. The fraction of flux due to the photospheric emission has been neglected since it is small and always within the photometric errors.

The dust temperature has been inferred from the $`([9.6][12])`$ color in all clusters, for homogeneity. Since 47 Tuc has been also observed in two additional narrow band filters (LW6\[7.7\] and LW8\[11.3\], cf. Sect. 2), the $`([7.7][11.3])`$ and $`([7.7][9.6])`$ colors can be also used to derive dust temperatures and to check the consistency with those derived from the $`([9.6][12])`$ color. Fig. 5 shows that only minor ($`<`$50 K) systematic differences exist between the three methods. The inferred values of $`T_{\mathrm{dust}}`$ are between $`200`$ and $`600`$ K and the dust equilibrium radii $`R_{\mathrm{dust}}`$ between 300 and 40 stellar radii, respectively. The corresponding dust masses range between a few 10-9 and a few 10-7 $`M_{}`$.

Fig. 6 shows that for lower dust temperature there is a somewhat larger dust mass, while there is not any evident correlation between the latter and the stellar luminosity.
The amount of material available to make dust in RGB stars must be proportional to the cluster metallicity. Only with helium burning stages could stars produce their own dust raw material; this should reach the surface only on the luminous AGB if then. Indeed, Fig. 6 shows that a large amount of dust has been only detected in the most metal rich clusters of our sample (47 Tuc and NGC 6388) while the two metal poor stars of M15 have low dust content.

## 5 Circumstellar mass loss

Empirical mass loss rates can be derived from dust masses by assuming a gas-to-dust ratio and a typical timescale for the outflow.

Since more dust should be present in more metal rich objects one also expects that the gas-to-dust ratio increases when metallicity decreases (cf. e.g. van Loon 2000). Assuming a gas-to-dust mass ratio of 200 in 47 Tuc, the value for the other clusters is scaled accordingly to their metallicity:

$$\rho =M_{\mathrm{gas}}/M_{\mathrm{dust}}=200\times 10^{([\mathrm{Fe}/\mathrm{H}]+0.76)}$$
(1)
An empirical estimate of the epoch of dust ejection can be inferred from the dust equilibrium radius and assuming an outflow velocity:

$$\tau _{\mathrm{outflow}}=R_{\mathrm{dust}}/v_{\mathrm{outflow}}$$
(2)
The outflow velocity should depend on the stellar luminosity and the gas-to-dust ratio (cf. e.g. Habing, Tignon & Tielens 1994). Larger luminosity leads to a higher momentum flow by the photons, and thus the kinetic energy and velocity in the flow could be higher. For lower gas-to-dust ratio the number of grains per molecule of gas is larger, and could lead to a higher outflow velocity. Adopting the Habing et al. (1994) relation:

$$v_{\mathrm{outflow}}L_{}^{0.3}\rho ^{0.5}$$
(3)
The dependence of the outflow velocity on the mass loss rate rapidly weakens as the latter increases (Habing et al. 1994), and it has been neglected in our computations.

Empirical mass loss rates averaged over the time required by dust to reach its equilibrium radius have been inferred using the relation:

$$dM/dt=\rho M_{\mathrm{dust}}/\tau _{\mathrm{outflow}}$$
(4)
In order to check the impact of our simplified modeling on the inferred mass loss rates, several simulations using the DUSTY code (Ivezić, Nenkova & Elitzur 1999) have been performed. A large grid of values for the overall optical depth at 12 $`\mu `$m, stellar temperatures, dust grain sizes, compositions and inner bound temperatures has been explored, using the analytic approximation to compute the density distribution for radiatively driven winds.

For homogeneity with our computations an average grain radius of 0.1 $`\mu `$m has been used, but very similar results (within 10%) are obtained with a standard grain size distribution between 0.005 and 0.25 $`\mu `$m. An average stellar temperature of 3500 K has been adopted. Temperature changes of $`\pm 500`$ K yield overall infrared excesses, mass loss rates and terminal outflow velocities which are only $`50`$%, $`30`$% and $`10`$% different, respectively. Lowering the inner boundary dust temperature from 1000 K (our reference value) to 500 K, gives a larger (but always within a factor of 2) variation in the inferred infrared excess, mass loss rates and outflow velocities. Varying the shell thickness from 10 to 10,000 times the shell inner radius produces negligible variations.

Fig. 7 shows comparisons between the observed infrared excesses, mass loss rates and outflow velocities inferred by our simplified modeling and some suitable DUSTY simulations, scaled to the stellar luminosity and gas to dust ratio of the giants in 47 Tuc. A good, simultaneous consistency (within a factor of two) between the empirical and simulated values have been obtained for an optically thin dust with $`\tau 0.05`$ and by assuming an outflow velocity of 14 $`\mathrm{km}\mathrm{s}^1`$ for a stellar luminosity $`L_{}=1000L_{}`$ and $`M_{\mathrm{gas}}/M_{\mathrm{dust}}=200`$, at the 47 Tuc metallicity.

For the stars with mid-IR excess in our sample the epochs of dust ejection are typically decades and total mass loss rates between 10-7 and 10-6 $`M_{}\mathrm{yr}^1`$ (cf. Table 3). The systematic uncertainty of these estimates, mainly due to the various model assumptions, is a factor of about 3.

Fig. 8 shows that the inferred mass loss rates are almost independent of the amount of circumstellar dust and stellar luminosity. Also there is not any obvious trend with metallicity within the scatter. Nevertheless, it is worth noticing that while the stars in the intermediate and metal rich clusters span the whole range of measured mass loss values, the two metal poor stars of M15 have rather high rates, despite the relative low content of dust (cf. Fig. 6). This is not surprising since one expects that among the most metal poor objects only those with large mass loss rates can have detectable dust. The blue tail HB morphology of M15 also shows that some of its stars (if genuine post-RGB stars) have lost almost their entire envelope on the RGB.

These inferred mass loss rates are finally compared with different empirical laws by Reimers (1975a,b), Mullan (1978), Goldberg (1979), and Judge & Stencel (1991) in Fig. 9. The reference formulae are taken from Catelan (2000, cf. Appendix and Figure 10), who revised the original ones taking advantage of the most extensive data set made available by Judge & Stencel (1991) and more recent data on intermediate and low mass giants. Our rates are about one order of magnitude larger than those predicted by these relations and do not show any clear dependence on stellar parameters. However, as pointed out by Catelan (2000) these relations are calibrated on Population I giants of relatively low luminosity and extrapolation to the short-lived evolutionary phases near the tip of the RGB might be fraught with uncertainty. These lower mass loss rates could be indeed consistent with the lower color excess (if any) possibly measured in 47 Tuc and $`\omega `$ Cen at $`M_{bol}2.5`$ (cf. Fig. 3 and Sect. 4).

## 6 Discussion

Dust grains are very efficient in absorbing and scattering starlight and could contribute to driving or accelerating stellar winds. However, as widely discussed in the recent review by Willson (2000), pure stationary flow models seem inadequate to describe the mass loss phenomena since a strongly dust-driven wind in an oxygen-rich star requires large light-to-mass ratios ($`L/M>10^4L_{}/M_{}`$). Only the coolest and most luminous AGB stars with direct signatures of dusty winds possibly satisfy this criterion. If winds are dust driven, non-steady phenomena like pulsational effects must be invoked to drive or enhance a stellar wind and to promote dust formation. The role of pulsating atmospheres in AGB mass loss has been considered for many years both on theoretical and observational basis (cf. Willson 2000 and references therein).

For stars on the RGB basic mass loss mechanisms are still quite unknown (cf. e.g. Judge & Stencel 1991) and matter of debate. There is observational evidence of small-amplitude and short period variability among giants close to the RGB-tip (cf. e.g. Welty 1985; Edmonds & Gilliland 1996), consistent with some oscillation activity (low-overtone radial or nonradial pulsations). Cacciari & Freeman (1983) and Gratton, Pilachowski & Sneden (1984), observed H$`\alpha `$ emission lines in hundreds of giants in a large sample of globular clusters. They initially argued that the H$`\alpha `$ emission was direct evidence for mass loss. However, interpreting H$`\alpha `$ emission can be complicated (Dupree 1986). For example, Dupree et al. (1994) suggest from its time variability that H$`\alpha `$ emission is due to atmospheric shocks which might drive the mass loss. In either interpretation the presence of H$`\alpha `$ emission in RGB spectra would generally indicate ongoing mass loss.

In the star in which mass loss is best understood, the Sun, the mass loss mechanism is ultimately tied to the magnetic field generated by the convective zone (e.g., Dupree 1986). Since cluster RGB stars have convective zones they could well have magnetic fields. There are no direct measurements of such fields, and interpreting signatures of chromospheric activity which might arise from magnetic fields is quite complex (e.g., Ayres et al. 1997). There is some recent evidence for RGB starspots which would imply magnetic activity (Stefanik et al. 2002). Soker (2000) and García-Segura et al. (2001) argue that magnetic activity might play a role in AGB mass loss on the basis of planetary nebula morphology.

Some non-canonical deep mixing in the upper RGB has been also invoked to produce excess luminosities and enhancing the mass loss close to the RGB-tip (cf. e.g. Sweigart 1997; Cavallo & Nagar 2000, Weiss, Denissenkov & Charbonnel 2000), but a solid physical picture of the overall impact on the red giant evolution is still lacking.

Can our observations give further information on the RGB mass loss process?

According to our mid-IR survey, strong mass loss processes seem to occur only close to the RGB-tip, during the last few million years of the RGB lifetime at $`M_{\mathrm{bol}}2.5`$ (Fusi Pecci & Renzini 1975; Salaris & Cassisi 1997). Moreover, only a fraction of giants detected by ISOCAM show strong mid-IR excess. This is especially clear in Fig. 3. Except near the very RGB-tip, stars with mass loss are scattered pretty uniformly along the RGB along with stars with no mass loss. Further, there is no smooth variation with the stellar luminosity. This suggests episodic mass loss processes. We will refer to the fraction of the time that mass loss is “on” as $`f_{\mathrm{ON}}`$. The fraction of time that the mass loss is “off,” or has dropped below our detection threshold, is $`f_{\mathrm{OFF}}`$. Since of the 52 ISOCAM sources 20 show evidence for mass loss, a first global estimate is $`f_{\mathrm{ON}}0.4`$. Considering the small numbers, there is no evidence that $`f_{\mathrm{ON}}`$ for the individual clusters varies from this global value.

Things are not so simple. In principle, all the stars brighter than $`M_{\mathrm{bol}}2.5`$ in each cluster should have been detected. In practice, due to the broad ISOCAM PSF, some stars are missed. First, while there is sometimes more than one star within an ISOCAM point source PSF, we selected only the brightest. Second, in the dense cores of M15 & NGC 6388 it is not possible to associate the ISOCAM objects with any near-IR stellar counterpart. The number of lost stars basically depends on the crowding, and thus on both the cluster distance and density—the larger their values the larger the number of stars in each pixel field of view, other things being constant.

To estimate the number of missed stars, we have computed the expected number of giants brighter than $`M_{\mathrm{bol}}2.5`$ within the ISOCAM field of view in two ways. First we used a theoretical approach: the fraction of light within the ISOCAM beam, was computed using a King model and the Specific Evolutionary Flux by Renzini & Buzzoni (1986), assuming an average lifetime $`6`$ Myr at $`M_{\mathrm{bol}}2.5`$ (Fusi Pecci & Renzini 1975; Salaris & Cassisi 1997). As an alternative we counted the number of stars detected in our $`J,K`$ survey within the same beam. The two methods provide similar results (within 20%). The numbers are given in Table 2. These values are significantly larger than the corresponding number of ISOCAM point sources detected in each cluster, as expected given the wide ISOCAM PSF.

How can the missed stars affect our estimate of $`f_{\mathrm{ON}}`$? The missed stars may or may not be losing mass. Only if they differ in some systematic way from the detected stars will they affect $`f_{\mathrm{ON}}`$. Since the brightest and reddest star was selected and since strong mass loss occurs only close to the tip of the RGB, the missed (fainter) stars may be less likely to have mass loss. A lower limit can be obtained by assuming that none of the missed stars has detectable mass loss. In that case for each cluster $`f_{\mathrm{ON}}`$ (number of stars with dust excess)/(total number of stars) (cf. last line of Table 2).

In 47 Tuc and NGC 362, where crowding does not severely affect statistics, $`f_{\mathrm{ON}}`$15%, while lower values have been inferred in M15 and NGC 6388 but in this cluster the very central core region was not sampled. In $`\omega `$ Cen 3–4 stars with dust excess should have been detected if $`f_{\mathrm{ON}}`$15% as in 47 Tuc and NGC 362, but taking into account that only $`<30`$% of the brightest giants are sampled (see Sect. 3), the expected number of detections drops down to $`1`$, fully consistent with the inferred zero value. The result of M54 is not statistical significant, since it is too distant (see also Sect. 3).

Broad time limits can be placed on the episodic variations. At average rates between 10-7 and 10-6 $`M_{}\mathrm{yr}^1`$ (cf. Table 3) each major mass loss event should occur on time scales somewhat shorter than the evolutionary times (5–6 Myr on the upper RGB) and longer than the time it takes the outflowing gas to reach the radii where we observe the dust, i.e., decades, whatever be the actual mass loss mechanism (several bursting episodes with multiple shells or more continuous flows). In principle such a lower limit to the episodic variation timescale can be as short as an instantaneous event, possibly at much higher mass loss rates, but our observations do not allow to constrain timescales with such a temporal resolution.

The total mass lost is crudely

$$\mathrm{\Delta }M\mathrm{rate}\times \mathrm{\Delta }t_{\mathrm{ON}}$$
(5)
where $`\mathrm{\Delta }t_{\mathrm{ON}}`$ is the total time interval in which a star experiences mass loss (either through many short-living episodes or a more continuous flowing).

In order to account for a total mass loss of about $`0.2M_{}`$ required by the models to explain the HB morphology, and by assuming an average $`\mathrm{rate}3\times 10^7M_{}\mathrm{yr}^1`$, $`\mathrm{\Delta }t`$ should be $`7\times 10^5\mathrm{yr}`$, well within the inferred limits and almost one order of magnitude smaller than the total evolutionary time on the upper RGB.

Another interesting result seen from Figs. 8, 9 is that there is not any clear dependence of mass loss on metallicity, although more data, in particular for the metal poor clusters, are needed to draw firm conclusions. Again this result is not unexpected. The earlier H$`\alpha `$ studies also suggested little metallicity dependence in mass loss rates and are consistent with $`f_{\mathrm{ON}}0.5`$. However our observed time scales are much longer than the day to day variation observed in H$`\alpha `$ emission. Both the time scale and $`f_{\mathrm{ON}}`$ could be consistent with RGB magnetic cycles. Further, studies of horizontal branch morphology (Catelan & de Freitas Pacheco 1993; Lee, Demarque & Zinn 1994) suggest at most a mild metallicity dependence on mass loss efficiency except possibly at the lowest metallicities. In a recent paper van Loon (2000) analyzing the mass loss rates of a sample of obscured AGB stars also finds very modest dependence on metallicity, but possibly for different reasons related to the pulsation activity (e.g. van Loon 2001).

The observed metallicity independence of mass loss coupled with the fact that mass loss for RGB stars can be far below the criterion for radiatively driven winds suggests that dust may not be related to the driving mechanism for RGB mass loss. Dust can simply be a tracer of the out flowing gas and not responsible for the outflow.

## 7 Conclusions

Our ISOCAM survey of the central region of six massive globular clusters provided mid-IR photometry for 78 sources with identified stellar counterparts. 52 sources are associated with bright giants close to the RGB-tip. Of these about 40% have strong mid-IR excess ascribed to the presence of dusty circumstellar envelopes. Correcting for stars not detected because of the low spatial resolution of ISOCAM, dusty envelopes are inferred around about 15% of the brightest giants ($`M_{\mathrm{bol}}2.5`$).

The inferred mass loss rates are in the range $`10^7<dM/dt<10^6M_{}\mathrm{yr}^1`$, assuming as reference values an outflow velocity of 14 $`\mathrm{km}\mathrm{s}^1`$ and a gas-to-dust mass ratio of 200 at the metallicity of 47 Tuc and for a stellar luminosity of 1000 $`L_{}`$.

The major astrophysical implications are:

* The mass loss occurs very near the RGB-tip and is episodic.
* The mass loss episodes must last longer than a few decades and less than a million years.
* There is no indication for a strict metallicity dependence of the frequency of mass loss occurrence.

The research has been supported by the Agenzia Spaziale Italiana through the grants 9920015AS and I/R/27/00. RTR is partially supported by NASA LTSA grant NAG 5-6403 and STScI grant GO-8709. We thank Jacco van Loon, the Referee of our paper, for the many suggestions and comments which improved the overall presentation and analysis of our results. We thank M. Bellazzini, G. Bono, C. Cacciari, M. Catelan, Andea Dupree, Jeff Linsky, and P. Persi for helpful comments and discussions.


format: text

text:
# PROBING GALAXY FORMATION WITH HIGH ENERGY GAMMA-RAYS

## 1 Introduction

The extragalactic background light (EBL) represents all the light that has been emitted by galaxies over the entire history of the universe. The EBL that we observe today is an admixture of light from different epochs, its spectral energy distribution (SED) distorted by the redshifting of photons as they travel to us from sources at different distances. It is therefore a constraint on both the intrinsic SEDs of the sources and their distribution in redshift. At present, there is more than a factor of two uncertainty in the amplitude of the EBL in the UV, optical, and near-infrared . The EBL in the mid-IR is even more uncertain. The far-IR background measured at $`\stackrel{>}{}100\mu `$m  represents at least half of the total energy in the EBL, yet the sources that produced it remain uncertain.

High energy $`\gamma `$-ray astronomy promises to help resolve these uncertainties by providing independent constraints on the EBL, in the mid-IR with $`E_\gamma `$ in the $`10`$ TeV energy range, and in the 0.1-3 $`\mu `$m range with $`E_\gamma 100`$ GeV via the new low-threshold instruments that will soon be available. High energy $`\gamma `$-rays from sources at cosmological distances are absorbed via electron-positron pair production on the diffuse background of photons that comprises the EBL. Thus, $`\gamma `$-ray observations of objects with known redshift and intrinsic spectral shape will constrain the EBL in these crucial wavelength regimes by measuring the optical depth of the Universe to photons of various energies. This in turn will help to constrain some of the most fundamental uncertainties in physical models of galaxy formation.

In order to illustrate this, in this paper we use a “forward evolution” approach, which attempts to model the essential features of galaxy formation using simple recipes. These semi-analytic models are set within the modern Cold Dark Matter (CDM) paradigm of hierarchical structure formation, and trace the gravitational collapse and merging of dark matter halos, the cooling and shock heating of gas, star formation, supernovae feedback, metal production, the evolution of stellar populations and the absorption and re-emission of starlight by dust. This machinery has been used extensively to predict optical properties of low-redshift galaxies, with good results (e.g., ; reviewed and extended in , hereafter SP and SPF). A semi-analytic approach was also used by Devriendt and Guiderdoni  to make predictions of counts and backgrounds in the mid-to-far-IR, with more detailed modelling of dust extinction and emission, but less detailed modelling of merging and star formation. We have now combined the strengths of these two approaches, by integrating the stellar SEDs and dust modelling of  into the galaxy formation SAM code of the Santa Cruz group.

Some parts of the “standard paradigm” of galaxy formation represented by our SAMs are relatively solid. For example, once a cosmological model and power spectrum are specified, it is straightforward to compute the gravitational collapse of dark matter into bound halos using $`N`$-body techniques, and analytic formalisms such as those used in our modelling  have been checked against these results . Within the range of values for the cosmological parameters allowed by existing observational constraints (i.e., $`\mathrm{\Omega }_{\mathrm{matter}}0.30.5`$, $`\mathrm{\Omega }_{\mathrm{matter}}+\mathrm{\Omega }_\mathrm{\Lambda }1`$, $`H_06080`$ km/s/Mpc; see e.g.  for a summary), these results do not change significantly. Similarly, modelling of gas cooling appears to be fairly robust and agrees well with hydrodynamic simulations . However, other aspects, notably the efficiency of conversion of cold gas into stars, the effect of subsequent feedback due to supernovae winds or ionizing photons, the stellar initial mass function (IMF), and the effects of dust, remain highly uncertain, and some predictions are quite sensitive to their details.

For example, SPF showed that the star formation history of the Universe and the number density of high redshift $`z\stackrel{>}{}2`$ “Lyman-break” galaxies (LBGs; e.g. ) may be quite different depending on whether star formation is primarily regulated by internal properties, such as gas surface density in a quiescent disk, or triggered by an external event such as an interaction. Because the largest samples of LBGs are primarily identified in the rest UV, model predictions are also quite sensitive to the high-stellar-mass slope of the IMF, and to dust extinction. At the other end of the spectrum is the sub-mm population detected by SCUBA, believed to be predominantly high redshift ($`z\stackrel{>}{}2`$) luminous and ultraluminous infrared galaxies (LIRGs and ULIRGs) powered by star formation rates of hundreds to thousands of solar masses per year (e.g., ). Theoretical predictions of the numbers and nature of these objects are highly sensitive to the same issues (the dominant mode of star formation, dust, the IMF), but provide a crucial counter-balance to the optical observations. However, the current mismatch between the sensitivity and spatial resolution of optical and sub-mm instrumentation has made it difficult to establish the connection between the two populations observationally.

The Milky Way, like most nearby galaxies, emits the majority of its light in optical and near-IR wavelengths; only about 30% of the bolometric luminosity locally is released in the far-infrared . This was generally believed to be typical of most of the starlight at all redshifts until the discovery of the far-IR part of the EBL by the DIRBE and FIRAS instruments on the COBE satellite, at a level ten times higher than the no-evolution predictions based on the local luminosity function of IRAS galaxies, and representing twice as much energy as the optical background obtained from counts of resolved galaxies . This result suggests that either the dust extinction properties of “normal” galaxies change dramatically with redshift, or a population of heavily extinguished galaxies (perhaps analogous to local LIRGs and ULIRGs) is much more common at high redshift than locally, or both. Some of these galaxies may have already been observed, at 15 $`\mu `$m by ISO , and at 850 $`\mu `$m by SCUBA .

Guiderdoni et al.  showed that their simplified semi-analytic model could reproduce the multi-wavelength data only if they introduced a population of heavily extinguished galaxies with high star formation rates, and with strong evolution of number density with redshift. This population was introduced ad-hoc by , but as discussed by these authors, by  (based on ), and also by SPF, the increasing importance of starbursts at high redshift, due to the increasing merger rate and higher gas fractions, is a natural mechanism to produce this population. The models of SPF contain a detailed treatment of mergers and the ensuing collisional starbursts, which has been calibrated against the merger rate in cosmological $`N`$-body simulations  and the starburst efficiency in hydrodynamical simulations . Moreover, they produced good agreement with observations of LBGs (e.g. ) and damped Lyman-$`\alpha `$ systems (SPF and ) as well as low redshift galaxies (SP). Therefore, it will be extremely interesting to see if these same models, when combined with the more sophisticated treatment of dust extinction and emission developed by Devriendt, Guiderdoni, and collaborators, will be able to simultaneously reproduce observations over the broad range of wavelengths and redshifts discussed above.

In the next section we briefly describe the ingredients of our models, and then present the results of the predicted EBL. Section 4 presents the implications for $`\gamma `$-ray attenuation, and §5 briefly discusses some alternative treatments and our own conclusions. The work summarized here is a brief, preliminary sample of the results which will soon be presented in a series of papers, now in preparation, on the EBL and its breakdown into various kinds of sources and on the implications for $`\gamma `$-ray astronomy.

## 2 Semi-analytic modelling

In this section we briefly describe the ingredients of our models. Readers can refer to SP and SPF for more details, and to  for a brief introduction.

Using the method described in , we create Monte-Carlo realizations of the masses of progenitor halos and the redshifts at which they merge to form a larger halo. These “merger trees” reflect the collapse and merger of dark matter halos within a specific cosmology (each branching in the tree represents a halo merging event — for examples, see e.g. ). We truncate the trees at halos with a minimum circular velocity of 40 km/s, below which we assume that the gas is prevented from collapsing and cooling by photoionization. Each halo at the top level of the hierarchy is assumed to be filled with hot gas, which cools radiatively and collapses to form a gaseous disk. The cooling rate is calculated from the density, metallicity, and temperature of the gas. Cold gas is turned into stars using several simple recipes, depending on the mass of cold gas present and the dynamical time of the disk. Supernovae inject energy into the cold gas and may expell it from the disk and/or halo if this energy is larger than the escape velocity of the system. Chemical evolution is traced assuming a constant yield of metals per unit mass of new stars formed. Metals are initially deposited into the cold gas, and may later be redistributed by supernovae feedback, and mixed with the hot gas or the diffuse (extra-halo) inter-galactic medium.

When halos merge, the galaxies contained in each progenitor halo retain their seperate identities either until they spiral to the center of the halo due to dynamical friction and merge with the central galaxy, or until they experience a binding merger with another satellite galaxy orbiting within the same halo. We take into account subhalo truncation due to tidal effects in the larger halo. All newly cooled gas is assumed to initally collapse to form a disk, and major (nearly equal mass) mergers result in the formation of a spheroid. New gas accretion and star formation may later form a new disk, resulting in a variety of bulge-to-disk ratios at late times.

For an assumed IMF, the stellar SED of each galaxy is then obtained using stellar population models. Here we use the multi-metallicity stellar SEDs of  for the Salpeter and Kennicutt IMF cases, and the solar metallicity GISSEL models  for the Scalo IMF. (We have found that using evolving metallicity rather than solar metallicity SEDs has a relatively small impact on the resulting EBL.) Dust extinction is modelled using an approach similar to that of . The optical depth of the disk is assumed to be proportional to the column density of metals. We then use a simple slab geometry where stars and gas are homogenously mixed, and assign a random inclination to each galaxy to compute the absorption. We use a metallicity dependent extinction curve, following .

All absorbed light is re-radiated at longer wavelength. The galactic dust emission spectrum is represented by a combination of three components: 1) hot dust (as in $`H_{\mathrm{II}}`$ regions), 2) warm dust (as in the diffuse $`H_\mathrm{I}`$), and 3) cold dust (as in molecular clouds). In the models of Devriendt et al. , these components are modelled as a mixture of polycyclic aromatic hydrocarbon molecules (PAH), very small grains, and big grains. Big grains may be either cold ($``$ 17 K), or heated by radiation from star-forming regions (as suggested by observations of typical local starburst galaxies like M82). A set of template spectra is then constructed for galaxies of varying IR luminosity, with admixtures of the various components selected in order to reproduce the observed relations between IR/sub-mm color and IR luminosity. A similar approach was used by , using a mixture of a typical Orion-like $`H_{\mathrm{II}}`$ spectrum and an $`H_\mathrm{I}`$ spectrum constructed to fit DIRBE observations of the diffuse ISM . Here, we use the more empirical emission templates of  (kindly provided in electronic form by E. Dwek), but we obtain very similar results with the models of .

The recipes for star formation, feedback, chemical evolution, and dust optical depth contain free parameters, which we set for each model (see SP) by requiring an average fiducial “Milky Way” galaxy to have a K-band magnitude, cold gas mass, metallicity, and average B-band extinction as dictated by observations of nearby galaxies.

Figure 1a shows the global star formation rate density for the two star formation recipes that we consider here. The “fiducial” model is the collisional starburst (CSB) model favored by SPF, in which bursts of star formation may be triggered by galaxy collisions. The “Late Star Formation” model is the Constant Efficiency Quiescent (CEQ) model of SPF, in which cold gas is converted to stars only in a quiescent mode with constant efficiency. This produces a star formation history similar to the models of the Durham group, in which the peak in the star formation history occurs at a more recent epoch ($`z1.5`$) than in the CSB model. For the CSB model, we consider three different choices of IMF: Scalo , Salpeter , and Kennicutt . These IMFs are graphed in Figure 2. For the CEQ model we show only the Kennicutt case. There is a noticable difference in the far-UV and the mid- to far-IR. The Scalo IMF produces less UV light relative to optical and near-IR light, compared to the Kennicutt and Salpeter IMFs, which produce more high mass stars than the Scalo IMF, and thus more ultraviolet light to be absorbed and re-radiated by dust in the far IR. In Fig. 1b we show the redshift evolution of the far-UV ($`2000\AA `$) luminosity density for these four models, compared with observations.222These models were also compared with the observed luminosity density from nearby galaxies, obtained by integrating the luminosity functions of galaxies resolved in recent redshift surveys at wavelengths ranging from 0.2 to 2.2 $`\mu `$m, in Fig. 3 of . However, the SAM outputs graphed there were inadvertently multiplied by a factor of $`h^40.25`$. In , we renormalized all the models by requiring that they all agreed with the K-band point at 2.2 $`\mu `$m. Here we do not do this since our current SAMs  use a corrected version  of the Press-Schechter formalism. The Scalo model falls short at all redshifts, and the CEQ model, which agrees at $`z=0`$, falls short at higher redshifts.333The Durham type models  of SPF also predict that LBGs have higher stellar mass than observations indicate , while predicted stellar masses from the CSB model of SPF are in good agreement with the observations. It is encouraging that our very simple model for dust extinction, which we normalized in the B-band at $`z=0`$, appears to yield the appropriate level of dust extinction in the UV at higher redshifts (SPF).

Recently, improved luminosity functions (LFs) have become available in optical bands from the Sloan Digital Sky Survey and in the K-band from 2MASS . We have found that our CSB model with Kennicutt IMF agrees well with all of these LFs when the average baryon fraction is $`f_b=0.1.`$ This model is also consistent with the number counts in the mid-IR (15 $`\mu `$m from ISOCAM ) and far IR (60 $`\mu `$m from IRAS, 175 $`\mu `$m from ISOPHOT) but not the sub-mm (850 $`\mu `$m from SCUBA ). The resulting EBL is similar to that from the Salpeter model discussed below.

## 3 The Integrated Extragalactic Background Light

Figure 3 shows the EBL produced by our four models, obtained by integrating the light over redshift (out to $`z=4`$) with the appropriate K-corrections due to cosmological redshifting. We compare this with a compilation of observational limits and measurements of the EBL. Fig. 4 presents the same four models and the same data, but with a linear rather than logarithmic vertical axis so that one can integrate the total energy in the EBL by eye. It is apparent that there is at least as much energy in the far-IR part of the EBL as in the entire optical and near-IR bands. For example, Puget and collaborators  estimated that the total energy in the EBL is between 60 and 93 nW m-2 sr-1, with between 20 and 41 nW m-2 sr-1 contributed by the optical and near-IR, and between 40 and 52 nW m-2 sr-1 coming from the far-IR. If the possible detection of the EBL at 60 $`\mu `$m by Finkbeiner et al.  were correct, that would further increase the far-IR EBL; however, it is very difficult to determine the EBL at 60 $`\mu `$m since the zodiacal light is so much brighter at that wavelength, and it was probably partly confused with the EBL .

The total energy in the EBL in units of critical density $`\rho _c`$ is $`\mathrm{\Omega }_{\mathrm{EBL}}=(4\pi /c)(I_{\mathrm{EBL}}/\rho _cc^2)=2.5\times 10^8I_{\mathrm{EBL}}h^2`$, where $`I_{\mathrm{EBL}}`$ is in units of nW m-2 sr-1. The total energy density in the EBL corresponding to the lower and upper estimates of  is $`\mathrm{\Omega }_{\mathrm{EBL}}=(3.65.5)\times 10^6(h/0.65)^2`$. Although the EBL includes energy radiated by active galactic nuclei (AGNs) as well as stars, it is unlikely that AGNs contributed more than a few percent of the total. This is because the total energy radiated by AGNs is $`E_{\mathrm{EBL}}^{\mathrm{AGN}}=\eta \rho _{\mathrm{BH}}c^2`$, where the efficiency of conversion of mass to radiated energy in AGNs is $`\eta 0.05`$. Correspondingly, $`\mathrm{\Omega }_{\mathrm{EBL}}^{\mathrm{AGN}}=\eta \mathrm{\Omega }_{\mathrm{BH}}(1+z_{\mathrm{BH}})^14.5\times 10^8h^1(\eta /0.05)[3/(1+z_{\mathrm{BH}})]\stackrel{<}{}0.02\mathrm{\Omega }_{\mathrm{EBL}}`$<sup>4</sup><sup>4</sup>4Updating , we have estimated $`\mathrm{\Omega }_{\mathrm{BH}}=(M_{\mathrm{BH}}/M_{\mathrm{spheroid}})\mathrm{\Omega }_{\mathrm{spheroid}}(1.5\times 10^3)(1.8\times 10^3h^1)`$, using the observed (loose) correlation  between a black hole mass and that of the galactic spheroid in which it is found, and the estimated cosmological density of spheroids . Note that the factor $`(1+z)^1`$ arises because of the dilution of the contribution of high-redshift sources due to the expansion of the universe. So for simplicity, in this paper we will neglect the contribution of AGNs to the EBL.

Several interesting features emerge from the comparison of our SAM models with the EBL data. In the UV to near-IR, the models are closer to the direct measures of the EBL obtained by  than to the lower limits from the Hubble Deep Field , although the Scalo IMF produces less light in the UV because it has fewer high-mass stars. As noted, our Kennicutt CSB SAM with $`f_b=0.1`$, which agrees well with the latest observed local luminosity density at $`z=0`$, produces an EBL close to the Salpeter one in Fig. 4. Of our four new EBL curves, the Late SF model and the fiducial Kennicutt model are also consistent with the DIRBE measurements at 140 $`\mu `$m. The Salpeter EBL lies a little more than 2$`\sigma `$ below the DIRBE measurement at 140 $`\mu `$m. The LateSF far IR is higher than the other models because its later star formation suffers less dilution due to the expansion of the universe. The models differ significantly in the mid-IR, $`1060\mu `$m, where the EBL can be probed by TeV $`\gamma `$-rays. The lower dotted curve in Fig. 3, representing our previous attempt  to model the EBL, is well below the 15 $`\mu `$m lower limit as well as the DIRBE measurements at longer wavelengths. As we stated in , we expected our EBL results to change as we improved our dust emission modelling. In addition to inclusion of the PAH features, the new dust emission model has more warm dust than the one used in .

We now discuss constraints from the TeV $`\gamma `$-ray observations.

## 4 Attenuation of high-energy $`\gamma `$-rays

Figure 5 shows the $`\gamma `$-ray attenuation predicted by the four $`\mathrm{\Lambda }`$CDM models considered here, for sources at redshifts $`z_s=0.03`$ and 0.10. All of the models predict rather little absorption at $`E_\gamma \stackrel{<}{}5`$ TeV for sources at $`z_s=0.03`$, but fairly sharp cutoffs above $`5`$ TeV, especially for the Late SF model. That model may be in conflict with the data from Mrk 501. For the blazars Mrk 421 and 501, both at $`z0.03`$, the synchrotron self-Compton (SSC) model, in which $``$keV synchrotron X-radiation from a very energetic electron beam is Compton up-scattered by the same electrons to produce the observed $``$TeV $`\gamma `$-rays, appears to explain both the keV-TeV spectra and their time variation (see, e.g.,  and references therein). Using a simplified SSC model and keV X-ray data to predict the unattenuated TeV spectrum of Mrk 501, Guy et al.  used CAT and HEGRA data to estimate the amount of $`\gamma `$-ray attenuation. They find that there is a rather good fit to the observed attenuation for the $`\mathrm{\Lambda }`$CDM-Salpeter EBL from our earlier work  when it is scaled upward by a factor of up to about 2.5 across the wavelength range 1-80 $`\mu `$m; this is the upper Salpeter curve on Fig. 3. Our new Salpeter curve appears to be rather consistent with this rescaling of our old Salpeter one, the Kennicutt curve may be a little high, and the Late SF curve appears to be definitely too high. As noted earlier, the new Salpeter EBL curve in Figs. 3,4 is similar to our latest $`f_b=0.1`$ Kennicutt SAM, which is in good agreement with the latest local luminosity functions in the optical and K bands from SDSS and 2MASS, and also in good agreement with IR number counts from IRAS and ISO satellites.

The large flares in Mrk 501 in spring 1999 allowed an accurate measurement of the gamma ray spectrum up to about 17 TeV, and indicated that the spectrum had an exponential cutoff at about 5 TeV . The flaring activity in Mrk 421 in early 2001 has now provided evidence for an exponential cutoff at about 4 TeV . The coincidence in the cutoffs (within observational uncertainties) for these two different extragalactic sources suggests that both are due to absorption via pair production on the EBL. In order to confirm this, it will of course be necessary to see similar cutoffs at lower energies for blazars at greater distances. There are already indications of this from 1426+42, a blazar at redshift $`z=0.13`$ (four times farther than Mrk 421), on which there is enough data from CAT , Whipple , and HEGRA  to begin to determine the spectrum from a few hundred GeV to several TeV. It will be very useful to measure the spectrum from this and other sources at comparable distances, such as PKS2155-304 at $`z=0.116`$, which will soon be possible with the next generation of atmospheric Cherenkov telescopes such as the H.E.S.S. array in Namibia, the CANGAROO-III array in Australia, and the VERITAS array in Arizona.

Assuming that the EBL is like the LateSF curve in Figs. 3,4, several authors (e.g. ) have argued that the TeV attenuation that this implies plus the observed spectrum of Mrk 501 leads to the requirement that the spectrum at the source have a strong upturn in its emitted flux above about 15 TeV, which would be very hard (although perhaps not impossible ) to understand, and they even suggest that one might have to abandon Lorentz invariance. However, if one instead assumes our Salpeter EBL in Fig. 4 and applies an SSC analysis to the HEGRA data from the Mrk 501 flares in 1997, the implied source spectrum is very reasonable, without an upturn at the high-energy end .

The compatibility of our new EBL calculations with the available data on TeV $`\gamma `$-ray attenuation is definitely worth further investigation. The results appear to be sensitive to the details of the models, raising the hope that they may be able to help answer important questions about star formation and dust reradiation, and also help to test the SSC modelling.

Figure 6 depicts the $`\gamma `$-ray “absorption edge,” the redshift of a source corresponding to an optical depth of unity, as a function of $`\gamma `$-ray energy. Travelling through the evolving extragalactic radiation field, $`\gamma `$-rays from sources at lower redshift suffer little attenution. The universe becomes increasingly transparent as $`E_\gamma `$ decreases, probing the background light at increasingly short wavelengths. (We are using the treatment of  to account for absorption of ionizing radiation by the Lyman alpha forest.) The models all have the same qualitative features, but differ significantly quantitatively. The location of the absorption edge is affected both by the assumed IMF and by the history of star formation. There is more absorption at most redshifts with the Kennicutt IMF because with a higher fraction of high mass stars, it is more efficient at producing radiation for a given stellar mass; there is more absorption nearby in the Late SF model because the starlight in this model is less diluted by the expansion of the universe. It is possible that measuring the transparency of the universe to $`\gamma `$-rays at $`100`$ GeV with a number of sources at various redshifts can provide a strong probe of star formation, although there are uncertainties due to extinction by dust.

## 5 Outlook

The semi-analytic modelling of the EBL described here follows the evolution of galaxy formation in time. Forward modelling is a more physical approach than backward modelling (luminosity evolution). Pure luminosity evolution (e.g., ) assumes that the entire evolution of the luminosity of the universe arises from galaxies in the local universe just becoming brighter at higher redshift by some power of $`(1+z)`$ out to some maximum redshift. It effectively assumes that galaxies form at some high redshift and subsequently just evolve in luminosity in a simple way. This is at variance with hierarchical structure formation of the sort predicted by CDM-type models, which appears to be in better agreement with many sorts of observations.

An alternative approach to modelling the EBL has been followed by Pei and collaborators , in which they find an overall fit to the global history of star formation subject to constraints from input data including the evolution of the amount of neutral hydrogen in damped Ly$`\alpha `$ systems (DLAS). Their first attempt , which was used as the basis for EBL estimates by , was somewhat misled by the sharp drop in the DLAS hydrogen abundance from redshift $`z3`$ to $`z2`$ reported in . With more complete data on DLAS (see, e.g., Fig. 14 of ) the $`z=3`$ point is lower and the neutral hydrogen abundance is almost constant from $`z=2`$ to 4. The latest paper by Pei et al.  takes a variety of recent data into account. Their approach is to follow the evolution of the total mass in stars, interstellar gas, and metals in a representative volume of the universe; they assume a Salpeter IMF. By contrast, the semi-analytic methods we use follow the evolution of many individual galaxies in the hierarchically merging halos of specific CDM models, here $`\mathrm{\Lambda }`$CDM. Despite the differences in approach, and the fact that  assumed $`\mathrm{\Omega }_m=1`$ and Hubble parameter $`h=0.5`$, their results are broadly similar to those from the semi-analytic approach (see their §4.4). In particular, their EBL is similar to our old results  for the Salpeter IMF. Our EBL results presented here are higher in the near-IR and more consistent with the direct determinations ; they are also higher in the mid-IR, probably mainly because of the warm dust and PAH features in our dust emission model. It will be interesting to see whether further development of the global approach of Pei et al. and of the semi-analytic approach lead to convergent results.

As our calculations show, the EBL, especially at $`\stackrel{<}{}1`$ $`\mu `$m and $`\stackrel{>}{}10`$ $`\mu `$m, is significantly affected by the IMF and the absorption of starlight and its reradiation by dust, as well as by the underlying cosmology. The cosmological parameters are becoming increasingly well determined by other observations. As data become available on $`\gamma `$-ray emission and absorption from sources at various redshifts, especially from the new generation of Atmospheric Cherenkov Telescopes and the new $`\gamma `$-ray satellites AGILE and GLAST, these data and their theoretical interpretation will help to answer fundamental questions concerning how and in what environments all the stars in the universe formed.

## Acknowledgments

I thank my collaborators Rachel S. Somerville, James S. Bullock, and Julien E. G. Devriendt. My work was supported by NASA and NSF grants at UCSC. I am grateful for a Humboldt Award, and I thank Leo Stodolsky for hospitality and Eckart Lorenz for enlightening discussions about $`\gamma `$-ray astronomy at the Max-Planck-Institut für Physik, München. I also thank Simon White for hospitality at MPI Astrophysics in Garching, and Heinz Völk for hospitality at MPI Nuclear Physics in Heidelberg. I thank Henric Krawczynski, Paolo Coppi, and Felix Aharonian for very helpful discussions of their SSC modelling. I thank Aldo Morselli for inviting me to give these lectures at the ISSS, and for his patience in waiting for me to send him the written versions.


text:
# Quark scattering at QGP/2SC interface

## I Introduction

The very interesting phenomenon arises at the junction between a normal-metal and a superconductor. The reflection of a conduction electron from the superconducting barrier, first elaborated in the paper , is known in the literature as the Andreev reflection. It was shown that the excitation in the normal metal outgoing from the reflection point has the sign of the momentum, the sign of the charge, and the sign of the effective mass opposite to that of the incoming electron (hole quantum numbers). This peculiar prediction was checked by tunneling experiments and ballistic (direct) observations .

In this letter we raise the interesting question of the interaction between the phase containing the free quarks (QGP) and the color-superconductor (CS). We shall show that the free quarks from the QGP phase falling at the CS barrier are reflected in the similar fashion as was described by Andreev in .

In the next section we develop the Bogolubov - de Gennes equations for superconducting phase. The third section contains the detailed discussion of the Andreev reflection. In the last chapter we make some comments about the influence of the Andreev reflection on the phenomena that can take place inside the Neutron Stars.

## II Bogolubov - de Gennes equations for 2SC phase

Our interest is the description of the phenomena that happens at the QGP/CS junction at moderate densities expected in the Neutron Stars. In this work we do not address the question of the gluon interaction with the CS phase thus we just concentrate on the fermionic degrees of freedom. The model of QGP is the gas of the free quarks. The quarks in QGP built the Fermi ladder, with the Fermi energy $`\mu `$ (of order $`400500`$ MeV in the Neutron Stars), where $`\mu `$ is a quark chemical potential. In the case of the two light quarks, which we consider here, the CS phase creates the 2SC superconductor (for the review of the subject see ). The essential physics of 2SC phase at moderate densities, is described by the effective four-fermion point interaction with the attractive pseudoscalar channel . This model leads to the Cooper instability of the Fermi sea which creates the new vacuum of condensed Cooper pairs. The physical excitations are fermionic quasiparticles separated from the vacuum by the energy gap $`|\mathrm{\Delta }|`$. The energy gap was calculated in the mean-field approximation in the variety of models giving the values $`40140`$ MeV which is about 10-30 per cent of the Fermi energy.

The effective hamiltonian at the mean-field level takes the form:

$$H=d^3x\left\{\psi _\alpha ^j\left(i\stackrel{}{\alpha }\stackrel{}{}+m\gamma _0\mu \right)\psi _\alpha ^j+\frac{\mathrm{\Delta }}{2}\psi _\alpha ^j(\tau _2)_{jk}t_{\alpha \beta }C\gamma _5\psi _\beta ^k\frac{\mathrm{\Delta }^{}}{2}\psi _\alpha ^j(\tau _2)_{jk}t_{\alpha \beta }C\gamma _5\psi _\beta ^k+(\mathrm{\Delta },m)\right\}$$
(1)
where $`m`$ is a quark mass and $`\psi _\beta ^j`$ is the field operator of the quark. The greek indices describe the color and the latin indices describe the isospin quantum numbers. The matrix $`\tau _2`$ is the antisymmetric Pauli matrix and $`t`$ is one of the three antisymmetric Gell-Mann matrices. For our calculations we chose it arbitrarily as the $`\lambda _2`$. The gap parameter $`\mathrm{\Delta }(\stackrel{}{r})`$ depends on the position. It vanishes inside the QGP phase and it takes a nonzero value in the 2SC phase. The equation of motion describing our system takes the form:

$`i\dot{\psi }_\beta ^j=\left(i\stackrel{}{\alpha }\stackrel{}{}+m\gamma _0\mu \right)\psi _\beta ^j+\mathrm{\Delta }(\stackrel{}{r})(\tau _2)_k^jt_\beta ^\sigma C\gamma _5\psi _\sigma ^{kT}`$ (2)
$`i\dot{\psi }_\alpha ^j=\stackrel{}{}\psi _\alpha ^ji\stackrel{}{\alpha }\psi _\alpha ^j\left(m\gamma _0\mu \right)+\mathrm{\Delta }(\stackrel{}{r})^{}\psi _\beta ^{kT}(\tau _2)_k^jt_\alpha ^\beta C\gamma _5`$ (3)
The equations (2) describe the interaction of the 6 “families” of quarks coupled through the non-zero value of the gap parameter. Using standard conventions (e.g. ), one can find that the up, ”red” quark couples to the down ”green” quark, the up ”green” to the down ”red”, whereas the ”blue” quarks remain free. For the first “family” the equations (2) read:

$`i\dot{\psi }_R^u=\left(i\stackrel{}{\alpha }\stackrel{}{}+m\gamma _0\mu \right)\psi _R^u\mathrm{\Delta }(\stackrel{}{r})C\gamma _5\psi _G^{dT}`$ (4)
$`i\dot{\psi }_G^d=\stackrel{}{}\psi _G^d\stackrel{}{\alpha }\psi _G^d\left(m\gamma _0\mu \right)\mathrm{\Delta }(\stackrel{}{r})^{}\psi _R^{uT}C\gamma _5`$ (5)
where $`u`$ and $`d`$ are isospin and $`R`$ and $`G`$ color degrees of freedom. The equations for the second “family” (up “green” and down “red” quarks) can be evaluate exactly the same way and the ”blue” quarks are not interesting in the first approximation. One can change the operator equations (4) to c-number equations by taking the expectation values $`f\mathrm{\Psi }_0|\psi _R^u|\mathrm{\Psi }_1`$ and $`g^{}\mathrm{\Psi }_0|\psi _G^d|\mathrm{\Psi }_1`$, where $`|\mathrm{\Psi }_0,|\mathrm{\Psi }_1`$ are ground and excited state respectively. In the QGP phase ($`|\mathrm{\Delta }|=0`$) the fields $`f`$ and $`g^{}`$ describe a “particle” and a “hole”. Let us now find the quasiparticle wavefunctions in the case of $`\mathrm{\Delta }=const`$. It is convenient to decompose the fermi fields in the plane wave basis. Remembering that the relativistic effects are suppressed in the presence of the Fermi sea one can neglect the contribution from the antiquarks:

$`f(t,\stackrel{}{r})={\displaystyle \underset{s}{}}u_s(\stackrel{}{q})\alpha _s(\stackrel{}{q})\mathrm{exp}(iEt+i\stackrel{}{q}\stackrel{}{r})`$ (6)
$`g^{}(t,\stackrel{}{r})={\displaystyle \underset{s}{}}d_s^{}(\stackrel{}{q})\beta _s^{}(\stackrel{}{q})\mathrm{exp}(iEt+i\stackrel{}{q}\stackrel{}{r}),`$ (7)
where $`s`$ describes the spin degrees of freedom, $`\alpha _s(\stackrel{}{q})`$ and $`\beta _s(\stackrel{}{q})`$ are some c-numbers. The energy $`E`$ describes the energy of the quark with respect to the Fermi surface in the absence of the gap parameter $`\mathrm{\Delta }`$. The Dirac bispinors $`u_s`$ and $`d_s`$ satisfy the equations:

$`(\stackrel{}{\alpha }\stackrel{}{q}+m\gamma _0\mu )u_s(\stackrel{}{q})=ϵ_qu_s(\stackrel{}{q})`$ (8)
$`d_s^{}(\stackrel{}{q})(\stackrel{}{\alpha }\stackrel{}{q}m\gamma _0+\mu )=d_s^{}(\stackrel{}{q})\overline{ϵ}_q,`$ (9)
and because momentum $`\stackrel{}{q}`$ can be in general the complex number one has to carefully distinguish between the left and the right eigenvectors. Inserting (6) into the equations (4), using bispinors algebraic relations and assuming constant value of the gap parameter one arrives at the Bogolubov - de Gennes equations for $`\alpha `$ and $`\beta `$ parameters:

$`E\alpha _{}(\stackrel{}{q})=ϵ_q\alpha _{}(\stackrel{}{q})+\mathrm{\Delta }\beta _{}^{}(\stackrel{}{q}),`$ (10)
$`E\beta _{}^{}(\stackrel{}{q})=ϵ_q\beta _{}^{}(\stackrel{}{q})+\mathrm{\Delta }^{}\alpha _{}(\stackrel{}{q})`$ (11)
and the similar set of equations with reverse spin projection on the momentum $`\stackrel{}{q}`$:

$`E\alpha _{}(\stackrel{}{q})=ϵ_q\alpha _{}(\stackrel{}{q})\mathrm{\Delta }\beta _{}^{}(\stackrel{}{q}),`$ (12)
$`E\beta _{}^{}(\stackrel{}{q})=ϵ_q\beta _{}^{}(\stackrel{}{q})\mathrm{\Delta }^{}\alpha _{}(\stackrel{}{q})`$ (13)
where $`ϵ_q=\sqrt{\stackrel{}{q}^2+m^2}\mu =\overline{ϵ}_q`$. We have two possible solutions of the uniform equations (10) (the case of (12) is similar):
for $`ϵ_q=\xi `$

$`\alpha _{}=\sqrt{{\displaystyle \frac{1}{2}}\left(1+{\displaystyle \frac{\xi }{E}}\right)}\mathrm{exp}\left(i{\displaystyle \frac{\delta }{2}}\right)`$ (14)
$`\beta _{}^{}=\sqrt{{\displaystyle \frac{1}{2}}\left(1{\displaystyle \frac{\xi }{E}}\right)}\mathrm{exp}\left(i{\displaystyle \frac{\delta }{2}}\right)`$ (15)
and for $`ϵ_q=\xi `$

$`\alpha _{}=\sqrt{{\displaystyle \frac{1}{2}}\left(1{\displaystyle \frac{\xi }{E}}\right)}\mathrm{exp}\left(i{\displaystyle \frac{\delta }{2}}\right)`$ (16)
$`\beta _{}^{}=\sqrt{{\displaystyle \frac{1}{2}}\left(1+{\displaystyle \frac{\xi }{E}}\right)}\mathrm{exp}\left(i{\displaystyle \frac{\delta }{2}}\right)`$ (17)
where $`\xi =\sqrt{E^2|\mathrm{\Delta }|^2}`$ and $`\delta `$ is a phase of the gap. Let us point out that only for $`E>|\mathrm{\Delta }|`$ the solution describes the propagating excitations. Finally the wavefunction describing the quasiparticle excitation of energy $`E`$ in the 2SC phase takes the form:

$`\psi (t,\stackrel{}{r})\left(\begin{array}{c}f\\ g^{}\end{array}\right)=\left(\begin{array}{c}\mathrm{exp}\left(i\frac{\delta }{2}\right)\sqrt{\frac{E+\xi }{2E}}\left(Au_{}(\stackrel{}{q_1})Bu_{}(\stackrel{}{q_1})\right)\\ \mathrm{exp}\left(i\frac{\delta }{2}\right)\sqrt{\frac{E\xi }{2E}}\left(Ad_{}^T(\stackrel{}{q_1})+Bd_{}^T(\stackrel{}{q_1})\right)\end{array}\right)\mathrm{exp}(iEt+i\stackrel{}{q_1}\stackrel{}{r})+`$ (22)
$`\left(\begin{array}{c}\mathrm{exp}\left(i\frac{\delta }{2}\right)\sqrt{\frac{E\xi }{2E}}\left(Cu_{}(\stackrel{}{q_2})+Du_{}(\stackrel{}{q_2})\right)\\ \mathrm{exp}\left(i\frac{\delta }{2}\right)\sqrt{\frac{E+\xi }{2E}}\left(Cd_{}^T(\stackrel{}{q_2})Dd_{}^T(\stackrel{}{q_2})\right)\end{array}\right)\mathrm{exp}(iEt+i\stackrel{}{q_2}\stackrel{}{r})`$ (25)
where $`q_{1,2}=\pm \sqrt{(\mu \pm \xi )^2m^2}`$ and $`A,B,C`$ and $`D`$ are arbitrary constants.

## III Andreev reflection

Let us consider the simple physical problem of the quark scattering on the superconducting plane surface placed at z=0 in space. If we approximate the boundary as the step function ($`\mathrm{\Delta }=0`$ for $`z<0`$ and $`\mathrm{\Delta }=const`$ for $`z>0`$) then we are looking for stationary solution of equations (4), independently for negative and positive $`z`$, supplemented with the condition that matches the wavefunction at $`z=0`$.

Let us consider in detail the situation when the “red” quark, with energy $`E`$ and spin projection $``$ falls at the plane boundary from the left. Then the wavefunctions (for $`\mathrm{\Delta }=0`$) take the form:

$$\psi _<(t,z)=\left(\begin{array}{c}u_{}(k)\mathrm{exp}ikz+Fu_{}(k)\mathrm{exp}ikz\\ Hd_{}^T(p)\mathrm{exp}ipz\end{array}\right)\mathrm{exp}\left(iEt\right)$$
(26)
The coefficients $`F`$ and $`H`$ describes the amplitude of the reflection of the $`u`$-particle and $`d`$-hole respectively. The first equation from (10) determine the value of the momentum of the $`u`$ quark as $`k^2=(\mu +E)^2m^2`$. In the similar fashion the second equation from (10) gives the value of the momentum of the reflected $`d`$ hole as $`p^2=(\mu E)^2m^2`$. For $`z>0`$ the quasiparticles excitation are described by the wavefunction $`\psi _>(t,z)`$ given by the expression (22) with $`B=D=0`$ (no spin flip in the scattering process). Let us note that for $`0<E<|\mathrm{\Delta }|`$, the momenta of the quasiparticles has to be complex number:

$$q_{1,2}=\pm \sqrt{\mu ^2m^2|\xi |^2\pm 2i|\xi |\mu }$$
(27)
Particularly simple form of momenta one can find in the massless limit. In that case $`q_{1,2}=\pm \mu +i|\xi |`$, where the signs are chosen as to describe the vanishing of the wavefunction in 2SC phase for $`z\mathrm{}`$. Now it is seen that for $`E<|\mathrm{\Delta }|`$ in the 2SC phase the quasiparticle wavefunction is exponentially decaying with the suppression factor $`|\xi |=\sqrt{|\mathrm{\Delta }|^2E^2}`$.

The continuity conditions matching the wavefunctions are of the form:

$`\psi _<(t,z=0)=\psi _>(t,z=0)`$ (28)
Using this condition one can find the amplitudes of the scattering process:

$`A=\sqrt{{\displaystyle \frac{2E}{E+\xi }}}\mathrm{exp}\left(i{\displaystyle \frac{\delta }{2}}\right)+O\left({\displaystyle \frac{1}{\mu }}\right)`$ (29)
$`C=\sqrt{{\displaystyle \frac{E\xi }{2E}}}{\displaystyle \frac{mE}{\mu ^2}}\mathrm{exp}\left(i{\displaystyle \frac{\delta }{2}}\right)+O\left({\displaystyle \frac{1}{\mu ^3}}\right)`$ (30)
$`F={\displaystyle \frac{m(E\xi )}{\mu ^2}}+O\left({\displaystyle \frac{1}{\mu ^3}}\right)`$ (31)
$`H=\sqrt{{\displaystyle \frac{E\xi }{E+\xi }}}\mathrm{exp}\left(i\delta \right)+O\left({\displaystyle \frac{1}{\mu }}\right)`$ (32)
In the limit where $`|\mathrm{\Delta }|,E,m<<\mu `$ only $`A`$ and $`H`$ contribute which means that only hole of the $`d`$-green quark can be reflected into the QGP phase and the only quasiparticle with momentum $`q_1`$ can propagate in the 2SC phase. It is interesting to point out that in the massless limit $`F=C=0`$ exactly. Additional insight we gain from the calculation of the probability current. From the equations (4) it follows that the probability current $`\stackrel{}{j}=\psi _R^u\stackrel{}{\alpha }\psi _R^u+\psi _G^d\stackrel{}{\alpha }\psi _G^d`$ is conserved. Using (22,29) one finds:

$$j_z=2\mu \{\begin{array}{ccc}0\hfill & \text{for}& \hfill E<|\mathrm{\Delta }|\\ \frac{2\xi }{E+\xi }+O\left(\frac{1}{\mu }\right)\hfill & \text{for}& \hfill E>|\mathrm{\Delta }|\end{array}$$
(33)
The result (33) has simple interpretation. If the quark $`u`$ have energy below the gap it can not excite quasiparticle inside the 2SC phase. In that case the matter waves does not propagate through the superconducting medium. However if the $`u`$ quark possesses energy above the gap it excites the quasiparticles with given transition coefficient. This result may be of importance for the physics of the Neutron Stars.

## IV Conclusions

In this paper we describe the phenomena of Andreev reflection at the junction between the cold quark-gluon plasma and superconducting 2SC phase. We found that the hole of different flavor and color than the incoming particle is reflected toward the QGP phase. Inside the 2SC phase the quasiparicles are excited, which for energy $`E`$ of the incoming particle below the gap $`|\mathrm{\Delta }|`$, penetrate the 2SC phase for depth of the order of $`1/\sqrt{|\mathrm{\Delta }|^2E^2}`$ whereas for the case of $`E>|\mathrm{\Delta }|`$ they propagate freely inside the superconducting phase. This result is similar to the situation one encounters in the condensed matter systems. This can be expected because in the high density QCD the purely relativistic phenomena are suppressed by the powers of quark chemical potential.

From the equation (33) one can see that the transport phenomena (like heat transport or density waves) is strongly affected by the presence of the interfaces inside the Neutron Stars. The suppression is exponential of the form $`\mathrm{exp}(|\mathrm{\Delta }|/T)`$ where $`T`$ is a temperature of the Neutron Star, usually much smaller than the expected superconducting gap. This phenomena certainly require more extensive analysisLet us note that probability current through the interface is not suppressed completely, because in the 2SC phase we have unpaired blue quarks..

The structure of the Neutron Stars can be complicated, thus more work is needed for better understanding of the Andreev reflection phenomena in superconducting QCD. In particular one has to consider other possible superconducting phases like CFL or other possible interfaces. Existence of such interfaces certainly influences the dynamics of matter flow inside the Neutron Stars.

Acknowledgement Let me specially thank Maciek A. Nowak for many important discussions. Also discussions with Krishna Rajagopal, Sanjay Reddy and Motoi Tachibana were of great value. This work was supported by Polish State Committee for Scientific Research, grant no. 2P 03B 094 19.


text:
# Flavour-symmetry breaking of the quark condensate and chiral corrections to the Gell-Mann-Oakes-Renner relation

## I Introduction

During the last two decades, QCD sum rules svz:79  and chiral perturbation theory ($`\chi `$PT) gl:84  have become indispensable tools in the phenomenology of particle physics. Both methods address the problem of dealing with QCD in the non-perturbative low-energy regime. Whereas in QCD sum rules the basic degrees of freedom are quarks and gluons and the non-perturbative region is approached from higher energies by including vacuum averages of composite operators, the condensates, $`\chi `$PT constitutes a low-energy expansion with the lowest lying hadrons as the fundamental fields. Nevertheless, a sound matching of QCD and low-energy effective theories still is an open problem of QCD phenomenology ppr:98 .

$`\chi `$PT is constructed such as to reflect the symmetry properties of QCD under chiral transformations. The QCD pattern of explicit chiral symmetry breaking is implemented by adding quark mass terms in the $`\chi `$PT Lagrangian, which introduces an additional free parameter $`B_0`$ at the leading order, besides the pion decay constant $`f_\pi `$. Since the $`\chi `$PT mass term contains the scalar quark density, it turns out that the parameter $`B_0`$ is related to the quark condensate.

The quark condensate in full QCD is not directly related to a physical quantity. Like the QCD coupling and the quark masses, its value depends on the renormalisation scale and on the renormalisation scheme employed in the calculation.111Contrary to the gluon condensate, the quark condensate, being an order parameter of chiral symmetry breaking, does not suffer from a renormalon ambiguity ben:99 . On the other hand, in $`\chi `$PT the information on the short-distance properties of the theory is hidden in the low-energy constants, and thus the question arises, how the quark condensates of QCD and $`\chi `$PT are related. The key to the answer lies in a low-energy theorem for the scalar and pseudoscalar correlation functions, which will be discussed in the next section.

In $`\chi `$PT, the dependence of the quark condensate on the renormalisation is reflected in an unphysical low-energy constant $`H_2^r`$, which appears at the next-to-leading order and is required for the renormalisation of the chiral expansion. $`H_2^r`$ thus comprises information on short-distance physics, which cannot be deduced from $`\chi `$PT alone. In section three, the constant $`H_2^r`$ will be estimated from flavour SU(3) breaking of the quark condensate. This flavour breaking can for example be extracted from QCD sum rules for the heavy-light $`B`$\- and $`B_s`$-meson systems.

Besides the flavour-symmetry breaking of the quark condensate, $`H_2^r`$ also governs the next-to-leading chiral correction to the well known Gell-Mann-Oakes-Renner (GMOR) relation gmor:68 . The GMOR relation constitutes the main source for numerical values of the quark condensate, assuming the light quark masses to be known. Making use of the estimate of $`H_2^r`$, in section four, also estimates for these chiral corrections and the quark condensate will be presented.

## II The scalar correlator

The central object which is investigated in the seminal version of QCD sum rules svz:79  is the two-point function $`\mathrm{\Psi }(p^2)`$ of two hadronic currents

$$\mathrm{\Psi }(p^2)i𝑑xe^{ipx}\mathrm{\Omega }|T\{j(x)j(0)^{}\}|\mathrm{\Omega },$$
(1)
where $`\mathrm{\Omega }`$ denotes the physical vacuum and in our case $`j(x)`$ will be the divergence of the flavour-changing vector (axialvector) current,

$$j(x)=^\mu (\overline{Q}\gamma _\mu (\gamma _5)q)(x)=i(Mm)(\overline{Q}(\gamma _5)q)(x).$$
(2)
$`M`$ and $`m`$ represent the masses of the quarks $`Q`$ and $`q`$ respectively, and the upper (lower) sign always corresponds to the scalar (pseudoscalar) case. Up to a factor $`(Mm)^2`$, $`\mathrm{\Psi }(p^2)`$ is just given by the scalar (pseudoscalar) two-point function. Furthermore, $`\mathrm{\Psi }(p^2)`$ satisfies a dispersion relation with two subtractions,

$$\mathrm{\Psi }(p^2)=\mathrm{\Psi }(0)+p^2\mathrm{\Psi }^{}(0)+p^4\underset{0}{\overset{\mathrm{}}{}}\frac{\rho (s)}{s^2(sp^2i0)}𝑑s,$$
(3)
where $`\rho (s)\text{Im}\mathrm{\Psi }(s+i0)/\pi `$ is the spectral function corresponding to $`\mathrm{\Psi }(s)`$.

A low-energy theorem relates the first subtraction constant $`\mathrm{\Psi }(0)`$ to the product of quark masses and quark condensates:

$$\mathrm{\Psi }(0)=(Mm)(\mathrm{\Omega }|\overline{Q}Q|\mathrm{\Omega }\mathrm{\Omega }|\overline{q}q|\mathrm{\Omega }).$$
(4)
The relation (4) can be derived by taking the divergence of the vector (axialvector) correlator, employing current commutation relations in addition bro:75 . Up to the next-to-leading order in the strong coupling constant, it was also verified in QCD perturbation theory bro:81 .

From the time-ordered product in eq. (1) and Wick’s theorem, one naturally obtains normal-ordered condensates. However, in the calculation of $`\mathrm{\Psi }(0)`$ also quartic mass corrections as well as higher-dimensional operators arise. On dimensional grounds, the operator corrections with dimension greater than four are accompanied by inverse powers of the quark masses thus being singular in the chiral limit. If the normal-ordered quark condensates are rewritten in terms of the non-normal-ordered minimally-subtracted condensates, the quartic mass as well as higher-dimensional operator corrections are cancelled, such as to yield the simple functional form of eq. (4) bro:81 ; sc:88 ; jm:93 ; cdps:95 .

The fact that the quark condensates in eq. (4) should be considered as non-normal-ordered implies that the subtraction constant $`\mathrm{\Psi }(0)`$ depends on the renormalisation scale and scheme. Whereas the product of quark mass times normal-ordered condensate is a renormalisation invariant quantity mqq , this is no longer true for the non-normal-ordered condensate. In the latter case, the renormalisation invariant involves additional quartic mass terms sc:88 . Numerically, the quartic mass terms only have some relevance for the strange quark. For the lighter up and down quarks they are completely negligible.

By considering QCD sum rules in which the subtraction constant $`\mathrm{\Psi }(0)`$ remains psi0 , it would in principle be possible to determine the quark condensate or flavour-breaking ratios of condensates. Here, care has to be taken that the quartic mass corrections on the left-hand side of eq. (3), which cancel the renormalisation dependence of $`\mathrm{\Psi }(0)`$ are included properly. Nevertheless, even though the various quantities related to $`\mathrm{\Psi }(0)`$ generally were found to be reasonable, since the sum rules in question are subject to very large higher-order perturbative corrections jop:01 , these results should be interpreted with caution.

## III The quark condensate in $`\chi `$PT

The status of the quark condensate in $`\chi `$PT can be clarified by calculating the scalar (pseudoscalar) correlator and the corresponding low-energy theorem also in this framework gl:84 . It is then found that the functional form of $`\mathrm{\Psi }(0)`$ exactly resembles eq. (4), and thus the quark condensate in $`\chi `$PT should be identified with the non-normal-ordered vacuum average of quark-antiquark fields.

As has been already discussed in the introduction, in $`\chi `$PT the quark condensate, as well as flavour-breaking ratios, depend on the unphysical low-energy constant $`H_2^r`$. An example of this dependence is displayed by the ratio

$$\frac{\overline{s}s}{\overline{q}q}=1+3\mu _\pi 2\mu _K\mu _\eta +\frac{8\mathrm{\Delta }_{K\pi }}{f_\pi ^2}(2L_8^r+H_2^r),$$
(5)
where $`\overline{q}q[\overline{u}u+\overline{d}d]/2`$ represents the isospin average of the up- and down-quark condensates, and in the following, the vacuum state $`\mathrm{\Omega }`$ is omitted. Furthermore, the $`\mu _P`$ are chiral logarithms which take the form

$$\mu _P=\frac{M_P^2}{32\pi ^2f_\pi ^2}\mathrm{ln}\frac{M_P^2}{\nu _\chi ^2},$$
(6)
with $`\nu _\chi `$ being the chiral renormalisation scale, $`\mathrm{\Delta }_{K\pi }=M_K^2M_\pi ^2`$, and $`L_8^r`$ is a physical low-energy constant in the next-to-leading order chiral Lagrangian gl:84 .

The relation (5) offers the possibility to determine the low-energy constant $`H_2^r`$, since the ratio $`\overline{s}s/\overline{q}q`$ can be obtained independently, e.g. in the framework of QCD sum rules. A rather reliable place to determine $`\overline{s}s/\overline{q}q`$ are sum rules for the leptonic decay constants of the heavy $`B`$\- and $`B_s`$-mesons jl:01  (and references therein), since the deviation of the flavour SU(3) breaking ratio $`f_{B_s}/f_B`$ from one is very sensitive to the quark condensate ratio $`\overline{s}s/\overline{q}q`$.

The ratio of the heavy meson decay constants $`f_{B_s}/f_B`$ is, unfortunately, not yet known experimentally. Thus we have to resort to a different framework to find results on this ratio, independent of QCD sum rules. Such results have been obtained in the framework of lattice QCD and a recent average yielded $`f_{B_s}/f_B=1.16\pm 0.04`$ ber:00 .222More recent results can be found in refs. fblat . Employing this value in the QCD sum rules, the corresponding flavour breaking ratio of the quark condensates turns out to be

$$\overline{s}s/\overline{q}q=0.8\pm 0.3.$$
(7)
Within the uncertainties, this result is in agreement to previous determinations ssoqq . Nevertheless, older results have not been included in eq. (7), since in these investigations normal-ordered condensates have been used and the corresponding quark mass corrections have not been included. A proper treatment of the determination of $`\overline{s}s/\overline{q}q`$ from QCD sum rules for other channels will be left for future work.

In addition to the ratio $`\overline{s}s/\overline{q}q`$, for the determination of $`H_2^r`$ also the low-energy constant $`L_8^r`$ is required. This constant can be determined from flavour-symmetry breaking of meson masses and decay constants, which at next-to-leading order in $`\chi `$PT satisfies the relation gl:84

$$\frac{1}{(R+1)}\frac{M_K^2}{M_\pi ^2}+\frac{f_K}{f_\pi }=\frac{3}{2}+\frac{3}{4}\mu _\pi \frac{1}{2}\mu _K\frac{1}{4}\mu _\eta +\frac{8\mathrm{\Delta }_{K\pi }}{f_\pi ^2}L_8^r.$$
(8)
Here $`R`$ represents the quark-mass ratio $`Rm_s/\widehat{m}=24.4\pm 1.5`$ leu:96 , and the value $`f_K/f_\pi =1.22\pm 0.01`$ pdg:00  will be used. Employing the relation (8), the constant $`L_8^r`$ is found to be

$$L_8^r(\nu _\chi =M_\rho )=(0.88\pm 0.24)10^3,$$
(9)
where the pion decay constant has been varied in the range $`f_\pi =92.4\pm 10.0\text{MeV}`$. Although the experimental uncertainty for $`f_\pi `$ is much smaller pdg:00 , the chosen error reflects the difference between the physical value and the leading order decay constant $`f_0=82\text{MeV}`$, such that uncertainties from higher orders in the chiral expansion should be largely taken into account. The result (9) is in agreement to the original value obtained by Gasser and Leutwyler gl:84  at the scale $`\nu _\chi =M_\eta `$.

The finding of eq. (9) for $`L_8^r`$ can also be compared with the very recent value $`L_8^r(\nu _\chi =M_\rho )=(0.62\pm 0.20)10^3`$ abt:01 , which has been obtained in $`\chi `$PT at $`𝒪(p^6)`$. Within the uncertainties both results are compatible, but it is seen that including the next-next-to-leading corrections in the chiral expansion effectively lowers the constant $`L_8^r`$. Since the presented analysis stays at the next-to-leading order, for consistency also the corresponding value (9) for $`L_8^r`$ has been utilised in the following.

Inserting the results of eqs. (7) and (9) into (5), it is possible to deduce an estimate on the unphysical low-energy constant $`H_2^r`$:

$$H_2^r(\nu _\chi =M_\rho )=(\mathrm{\hspace{0.17em}3.4}\pm 1.5)10^3.$$
(10)
The negative value for $`H_2^r(\nu _\chi =M_\rho )`$, found in eq. (10), is in contrast to estimates of the same quantity in the framework of $`\chi `$PT with explicit inclusion of resonance fields egpr:89 . Assuming saturation of the low-energy constants in the chiral Lagrangian by meson resonances, reasonable values for the physical constants $`L_i^r`$ are obtained, whereas in the case of $`H_2^r`$, saturation with scalar resonances yields the relation $`H_2^r=2L_8^r`$, which is clearly violated by the result (10).333In fact, employing the relation $`H_2^r=2L_8^r`$, lead the authors of abt:01  to deduce $`\overline{s}s/\overline{q}q=1.69`$. Since $`H_2^r`$ and $`2L_8^r`$ have the same dependence on the chiral scale $`\nu _\chi `$, it is even impossible to find a particular scale at which the relation $`H_2^r=2L_8^r`$ would hold.

## IV Chiral corrections to the Gell-Mann-Oakes-Renner relation

The unphysical low-energy constant $`H_2^r`$ also appears in the next-to-leading order chiral corrections to the GMOR relations gmor:68 ; gl:84 ,

$`4\widehat{m}\overline{q}q`$ $`=`$ $`\mathrm{\hspace{0.17em}2}f_\pi ^2M_\pi ^2(\mathrm{\hspace{0.17em}1}\delta _\pi ),`$ (11)
$`(m_s+\widehat{m})[\overline{s}s+\overline{q}q]`$ $`=`$ $`\mathrm{\hspace{0.17em}2}f_K^2M_K^2(\mathrm{\hspace{0.17em}1}\delta _K),`$ (12)
where $`\widehat{m}(m_u+m_d)/2`$ is the isospin average of the light up- and down-quark masses, and $`\delta _\pi `$ as well as $`\delta _K`$ comprise the higher-order chiral corrections. The left-hand sides of eqs. (11) and (12) just correspond to minus $`\mathrm{\Psi }(0)`$ in the pseudoscalar case with $`(ud)`$ and $`(sq)`$ quark flavours respectively. At the next-to-leading order, these corrections have been found to be gl:84 :

$$\delta _\pi =4\frac{M_\pi ^2}{f_\pi ^2}(2L_8^rH_2^r)\text{and}\delta _K=\frac{M_K^2}{M_\pi ^2}\delta _\pi .$$
(13)
The results for $`H_2^r`$ and $`L_8^r`$, obtained in the last section, now allow for an estimation of the chiral corrections $`\delta _\pi `$ and $`\delta _K`$, which take the values

$$\delta _\pi =0.047\pm 0.017\text{and}\delta _K=0.61\pm 0.22.$$
(14)
With roughly 5%, the size of the correction $`\delta _\pi `$ appears very reasonable, whereas $`\delta _K`$ in the strange quark sector turns out rather large, due to the enhancement factor $`M_K^2/M_\pi ^2`$, although the uncertainties are also quite big. If, on the other hand, the relation $`H_2^r=2L_8^r`$ which follows from scalar resonance saturation egpr:89  would be correct, both chiral corrections $`\delta _\pi `$ and $`\delta _K`$ would have to vanish.

The knowledge of the light quark masses also allows for a determination of the light quark condensate from the first GMOR relation (11). Using the recent result $`(m_u+m_d)(2\text{GeV})=8.1\pm 1.4\text{MeV}`$ jop:01 , together with the estimate (14) for $`\delta _\pi `$, one arrives at

$$\mathrm{\Omega }|\overline{q}q|\mathrm{\Omega }(2\text{GeV})=(267\pm 16\text{MeV})^3,$$
(15)
which can be considered as an update of previous determinations of the light-quark condensate $`\overline{q}q`$. Since it is still more common to quote the quark condensate at a scale of $`1\text{GeV}`$ the corresponding value $`\overline{q}q(1\text{GeV})=(242\pm 15\text{MeV})^3`$ is also provided. The latter results are in good agreement with direct determinations of $`\overline{q}q`$ from QCD sum rules dn:98  and lattice QCD qqlat .

## V Conclusions

From the structure of the low-energy theorem of eq. (4) for the scalar (pseudoscalar) two-point function, it is possible to infer that the quark condensate in $`\chi `$PT should be identified with the non-normal-ordered vacuum average of quark-antiquark fields. Contrary to the corresponding case of normal-ordered condensates, the product of quark mass times non-normal-ordered condensate, and thus $`\mathrm{\Psi }(0)`$, are not renormalisation invariant. In $`\chi `$PT the dependence of $`\mathrm{\Psi }(0)`$ on the short-distance renormalisation scale and scheme is reflected in the dependence on the unphysical low-energy constant $`H_2^r`$.

Because the low-energy constant $`H_2^r`$ depends on the renormalisation at short-distances, its value cannot be deduced from relations to physical quantities alone. One possible quantity to determine $`H_2^r`$ is the SU(3) flavour-breaking ratio of the quark condensates $`\overline{s}s/\overline{q}q`$. The estimate presented in eq. (10) has been obtained on the basis of $`\overline{s}s/\overline{q}q`$ from QCD sum rules for the heavy $`B`$\- and $`B_s`$-mesons, making additional use of results on the ratio $`f_{B_s}/f_B`$ from lattice QCD.

The obtained value for $`H_2^r`$ also allows to estimate the chiral corrections to the GMOR relations of eqs. (11) and (12). It is found that the first GMOR relation in the light up- and down-quark sector receives a moderate negative correction of the order of 5%, whereas the corresponding correction in the strange-quark sector is enhanced by a factor $`M_K^2/M_\pi ^2`$, and thus turns out to be rather large. Together with values for the light quark masses, the first GMOR relation makes it possible to determine the light quark condensate $`\overline{q}q`$, and an updated value has been presented in the last section.

Since the low-energy constants of $`\chi `$PT contain the physics of energies higher than the pseudoscalar meson masses, it is possible to derive effective Lagrangians including higher resonance fields, which then yield contributions to the chiral constants egpr:89 . For the physical constants $`L_i^r`$, very reasonable estimates were obtained in this manner. For the unphysical constant $`H_2^r`$, however, it was shown above that the result $`H_2^r=2L_8^r`$ from the resonance Lagrangian is not satisfied. It should already be clear that this relation has to receive additional contributions, since $`H_2^r`$ depends on the short-distance renormalisation, whereas this is not the case for $`L_8^r`$.

Nevertheless, the uncertainties in the ratio $`\overline{s}s/\overline{q}q`$ are still sizeable, so that no definite conclusions can be reached at this stage. Further work should be devoted to the extraction of $`\overline{s}s/\overline{q}q`$ from QCD sum rules, particularly paying attention to the proper definition of the quark condensate as well as to the corresponding higher-order strange-mass corrections. In addition, an analysis analogous to the one presented here at the next-next-to-leading $`𝒪(p^6)`$ in $`\chi `$PT could provide further useful information on the questions raised above, since also the uncertainty resulting from these higher orders plays an important role.

###### Acknowledgements.

It is a great pleasure to thank J. Prades, H. G. Dosch, A. Pich and E. de Rafael for helpful discussions. The author would also like to thank the Deutsche Forschungsgemeinschaft for their support.