text:
# Molecules in Southern Dense Cores I: Mopra and SEST Millimetre-Wave Observations

## 1 Introduction

Molecular-line emission provides a powerful method for investigating the events associated with star formation. Molecular emission can be used to find the temperature, density and kinematics of gas from which stars form, while emission from different molecules traces gas of differing physical conditions. The theoretical and observational underpinning for this type of study has been the work of many investigators over a thirty year period; a good review can be found in van Dishoeck & Blake (1998).

The southern part of the Galactic plane includes the innermost parts of the galaxy as well as major dense molecular clouds and star-forming regions, many of which are difficult to observe with northern hemisphere telescopes due to their southern declinations. Consequently molecular line emission, molecular abundances and physical conditions have not been extensively investigated in clouds at declinations south of -30 degrees.

In the radio spectrum, molecular lines occur primarily at millimetre and submillimetre wavelengths. The commissioning of a millimetre wave receiver on the Mopra 111The Mopra telescope is part of the Australia Telescope which is funded by the Commonwealth of Australia for operation as a National Facility managed by CSIRO. antenna of the Australia Telescope in 1994 provided the opportunity to observe molecular transitions at 3 mm with high velocity resolution, making it possible to produce a statistical database of molecular spectral-lines and molecular abundances in southern molecular clouds.

In the study reported here, a sample of 24 molecular clouds with declinations south of -30 degrees has been observed with the Mopra Telescope in the 3-mm transitions of the molecules CO, CS, C2H, CH3OH, HCN, HNC, HCO+, HC3N, OCS, HNCO and SO and several of their isotopomers. In addition, complementary observations of 2-mm and 3-mm transitions have been obtained with the Swedish-ESO Submillimetre Telescope (SEST) 222The SEST is operated jointly by the Swedish Natural Science Research Council and the European Southern Observatory.

## 2  Observations and Equipment

### 2.1 Equipment

#### 2.1.1  The Mopra Telescope

The 22-m ‘Mopra’ antenna of the Australia Telescope (operated by the Australia Telescope National Facility, CSIRO and the University of New south Wales) is situated near Coonabarabran in New South Wales, Australia. The telescope is at latitude 31 degrees south and has an elevation of 850 metres above sea level. In 1994 the telescope was fitted with a millimetre-wave receiver allowing it to operate at frequencies between 85 and 116 GHz. The antenna has a ‘shaped’ parabolic surface area to increase its forward gain and until 1999 it had a solid reflecting surface extending over a diameter of only 15 metres. Thus at millimetre wavelengths the reflecting surface had an effective diameter of 15 m, resulting in a half-power-beam-width that varied from 47 arcsec to 35 arcsec across the range 85 - 116 GHz. The antenna was resurfaced in 1999, 333 Resulting from funding contributions from the University of New South Wales, the solid surface was extended over the full diameter in 1999, extending the millimetre-wave collecting area to 22 metres. extending the effective millimetre-wave diameter to 22 metres, giving a half-power beamwidth varying from 35 to 24 arcsec over the same range of frequencies.

The Mopra mm-wave receiving system includes a cryogenically-cooled SIS (Superconductor Insulator Superconductor) amplifier that can be tuned to either single or double side-band mode. For single sideband operation, the response of the unwanted sideband can be decreased to at least 10 dB below the other sideband response using a calibration signal injected into both sidebands. A polarization splitter is used to separate the incoming signal into two channels that can be tuned separately, allowing the observation of two frequencies simultaneously. The ‘backend’ is a digital autocorrelator, providing two simultaneous outputs with available band widths ranging from 4 to 256 MHz, which can be split into a maximum of 4096 channels. For the duration of this survey, the only observing method available was position switching: i.e. the observations were made by switching between the on-source position and a reference position with insignificant emission.

Both the antenna and the digital correlator have features in common to those used in the Australia Telescope Compact Array, and further details can be found in the Special Edition of the Journal for Electrical and Electronic Engineering, Australia, Vol 12, No 2 (1992). Technical details specific to the Mopra antenna can be found at www.narrabri.atnf.csiro.au/mopra/mopragu.pdf.

An important aspect of the study reported here (involving observations collected between 1994 and 2000) was the formulation of an intensity calibration scheme for the Mopra telescope. This calibration scheme is discussed in the appendix.

#### 2.1.2  The SEST Telescope

The SEST is situated at the European Southern Observatory at La Silla in Chile, about 600 km north of Santiago. The SEST has a Cassegrain design and a parabolic collecting surface with a diameter of 15 metres, giving a half power beam width ranging from 57 arcsec at 86 GHz to 15 arcsec at 346 GHz.

The SEST is fitted with three SIS receivers that collectively cover the frequencies 78 - 117 GHz, 128 - 170 GHz, 215 - 270 GHz and 320 - 363 GHz. A polarization splitter allows the observation of two different frequencies simultaneously, although the combinations of frequencies are limited to 3-mm and 2-mm, or 3-mm and 1.6-mm.

Three Acousto-Optical Spectrometers (AOS) provide backends to the receivers. The high resolution AOS has a maximum bandwidth of 86 MHz, which can be divided into up to 2000 channels.

For more instrumentation details see Booth et al. (1989) or The SEST Handbook (1997).

### 2.2 Selection of the Survey Sample

The aim of the selection process was to select around 25 southern molecular clouds with strong molecular line emission and evidence of complex chemistry for detailed study, and involved a two-step process. Firstly, previous molecular-line surveys of the southern Galactic plane (Whiteoak & Gardner 1974, Whiteoak & Gardner 1978, Gardner & Whiteoak 1978, Batchelor et al. 1981, Dickinson et al. 1982, Gardner & Whiteoak 1984, Peters et al. 1986) were used to choose 50 candidate molecular clouds south of declination -30 degrees. Secondly, these candidates were observed with the Mopra radio telescope between September and November 1994, in the 3-mm transitions of CS, HCN and HCO+ at frequencies of 98.0, 88.6 and 89.2 GHz respectively. These transitions were selected firstly because these molecules are generally abundant in molecular clouds and secondly because they were within the limited range of frequencies available on the Mopra telescope at that time. These observations were used to locate the line emission peaks, and select 24 of the brighter clouds for detailed study. The sample included several (5) dark clouds and also some objects showing unusual intensity ratios of the three transitions. Note that these observations were made during commissioning of the receiver, and before the development of better positional calibration and intensity calibration. Consequently, later observations of these transitions were obtained for use in the quantitative study (see table 3).

The list of 24 Galactic molecular clouds selected is given in table 1, with the observed positions in both B1950 and J2000 coordinates. The two molecular clouds Orion MC and M 17, which are near the celestial equator and well studied by radio telescopes in the northern hemisphere, have been included for comparison and calibration. Four separate pointing positions were chosen within the complex NGC 6334, and two positions within G291.3-0.7 (NGC 3576, RCW 57). Therefore, including the two reference clouds (Orion and M 17), there are 30 target positions in table 1.

The exact positions were chosen based on 13CO 1-0 observations at 110 GHz with the Mopra telescope in September 1995. The clouds were mapped over a region of a few square arcminutes to determine the position of the peak, and the angular size of the 13CO emission. The beam size was 37 arcsec, and the pointing errors 10-15 arcsec, with occasional discrepancies of up to 20 arcsec; the positions listed in table 1 are believed to be within 20 arcsec of the 13CO condensation peak. Later observations centred on these positions (table 3) have pointing errors better than 10 arcsec.

### 2.3 Survey Observations

#### 2.3.1  Mopra Observations

Molecular transitions between 85 and 116 GHz were observed with the Mopra telescope between September 1994 and September 1996 with the 15-metre solid millimetre surface, and in March and September 2000, after the resurfacing to 22 metres. Emission from 11 different molecules was observed (CO, CH3OH, C2H, CS, HCN, HC3N, HCO+, HNC, HNCO, OCS and SO), generally in more than one transition. Observations of more rare, isotopically substituted species of the molecules CO, CS, HCN, HCO+ and HNC were also obtained (13C for 12C and 18O for 16O). The actual transitions observed are listed in table 2, together with the date of observation, the frequency and the telescope used.

The single-sideband receiver temperature varied between 85 and 130 K. All observations were obtained with a bandwidth of 64 MHz divided into 1024 channels, giving a velocity resolution between 0.44 and 0.32 km s-1 after Hanning smoothing.

Observations were made in position switching mode, generally with a separation of 30 arcmin; however, for some extended sources such as the dark clouds, 60-arcmin separations were used. Periodic observations of SiO masers of known position were used to correct the pointing model, giving a pointing accuracy of better than 10 arcsec for observation obtained in 1995 or later. This includes all observations reported in table 3.

#### 2.3.2  SEST Observations

The SEST was used in June 1999 to observe transitions between 84 and 147 GHz. The system temperatures were 170 - 250 K at 3 mm and 230-400 K at 2 mm wavelength bands. The beamwidth varied from 59 arcsec at 85 GHz to 35 arcsec at 147 GHz. Note that although the effective diameter of SEST was the same as Mopra (15 m) the Mopra dish was ’shaped’ with different illumination resulting in a smaller beamwidth for a given frequency. A bandwidth of 42 MHz was used, divided into 1000 channels, giving velocity resolution between 0.30 and 0.17 km s-1 after Hanning smoothing.

The observations were made using two different observing modes. Generally, dual-beam switching mode was utilised with a throw of 11 arcmin between the on-source and reference positions. However, for extended sources, where line emission was likely to be present in the reference position, frequency switching was used, with a frequency throw of 4 MHz. Periodic observations of SiO masers were used to check focus and pointing, giving a pointing accuracy better than 5 arcsec. Intensity calibration of the observations was performed on-line by the chopper wheel method (see Appendix), with the spectra obtained being on the TA* scale.

### 2.4 Data Reduction

All spectra were reduced using the CLASS (Continuum and Line Analysis Single-dish Software) package of the GILDAS (Grenoble Image and Line Data Analysis Software) working group.

The spectra were Hanning-smoothed to increase the signal-to-noise ratio before baseline subtraction and Gaussian fitting to features. The spectral profiles generally had complex shapes, so one or more Gaussian components were fitted.

For Mopra spectra, the peak and integrated intensities of the Gaussian components were then corrected for atmospheric attenuation and scaled to the SEST corrected antenna temperature ($`T_A^{}`$) scale (see Appendix for more details on intensity calibration of Mopra). The SEST corrected antenna temperatures ($`T_A^{}`$) were then converted to main-beam brightness temperatures ($`T_{MB}`$) using the SEST main beam efficiency ($`B_{eff}/F_{eff}`$) which varied from 0.75 to 0.66 between 87 and 147 GHz (Booth et al., 1989, SEST Observers Handbook).

## 3 Results

Table 3 lists the observed molecular transitions, and the results of the Gaussian-component fits for each molecular cloud. An excerpt from table 3, showing results for the molecular cloud G265.5+1.5, is given in the body of this paper. The complete table can be found in the electronic version of the paper and on the VizieR data base of astronomical catalogues at the Centre de Donn$`\stackrel{´}{e}`$es astronomiques de Strasbourg (CDS) website. The first two columns list the molecule and the transition observed. The next three columns are the Gaussian fit parameters (central velocity, full-width at half maximum and peak expressed as main beam brightness temperature TMB). Multiple Gaussians were used to fit complex profiles. However, interpreting components in terms of separate clouds is difficult because in some cases, double-peaked profiles fitted as two Gaussians may be a broad emission line with a self-absorption dip. The velocities of the HCN hyperfine components are for the J = 2-1 component with the F = 1-0 and F = 1-1 line centres fitted using the same velocity scale. These latter hyperfine transitions are shown in brackets in the table. (This is because the components are often blended, so it was not feasible to velocity correct the hyperfine components separately.) The 6th column shows the integrated line emission ($``$ TMB dv) for individual Gaussian components, while the 7th column shows the total line emission in the transition, obtained by adding the individual Gaussian components; in the case of HCN this total includes the emission in the three hyperfine components. The main beam brightness temperature (TMB), is the true brightness temperature of the source, if the source just fills the main beam. If the source angular size is known, then this temperature can be deconvolved with the source size to give the actual source brightness temperature.

The 8th column gives the uncertainty in the integrated brightness temperature (column 7). For Mopra observations, the scatter of repeated observations of a source was used to make an estimate of the accuracy of the integrated brightness temperature of the lines. The error in integrated brightness was considered to be made up of two parts adding in quadrature - a constant percentage error B, which dominates at high brightness, and a constant absolute error A, due to noise in the spectra, which dominates at low brightness. From the scatter of repeated observations of strong lines ($``$TMB dv $`>`$ 5 K km s-1) we determine B = 12 %, and from repeated observations of weak lines ($``$TMB dv $`<`$ 5 K km s-1) we determine A = 0.24 K km s-1. For a typical linewidth of 4 km s-1, this value of A corresponds to peak brightness 0.06 K and is consistent with the contribution to line fitting from the noise in the spectra.

The calibration uncertainties of the SEST observations are assumed to be 10 % at 3 mm and 15 % at 2 mm (SEST Observers Manual 2003; Heikkila 1998), and the total uncertainty the sum in quadrature of the calibration uncertainty and the formal errors in fitting the lines using the CLASS fitting routine (GAUSS). The Mopra and SEST uncertainties turn out to be similar (listed in column 8 of table 3).

## 4 Discussion

Some general properties of the line emission and profiles are discussed in this section, including the ratios of the emission from different isotopomers of several molecules. In addition, the line intensities have been used to calculate molecular abundances using both LTE and LVG models, which we present in a separate paper (Hunt Cunningham et al. 2004).

### 4.1 Spectral Line Intensities and Profiles

Figure 1 shows CS 2-1 spectra towards each position listed in table 1. Each spectrum is plotted over a velocity interval of 30 km s-1. This transition was chosen to illustrate the characteristic line profile in each molecular cloud as it is ubiquitous, generally has good signal to noise, and tends to have a ‘typical line profile. On the other hand, the CO spectra for each molecular cloud often show emission features covering a significantly wider range of velocities than found for other molecules, as can be seen in figure 2, which shows spectra for a number of different transitions towards the molecular cloud G265.1+1.5 (RCW 36). The spectra towards this molecular cloud illustrate features common to many of the molecular clouds.

Self-absorption is evident in the 12CO spectra towards about a third of the molecular clouds, and can also be seen in some spectra of 13CO and C18O. Transitions that are believed to be self-absorbed are marked in table 3 by an asterisk. Self-absorption is also evident in some clouds in transitions of CS, HCN, HCO+ and SO, but not generally in transitions of OCS, HC3N and CH3OH, implying that these latter transitions are more likely to be optically thin. In addition to the presence of self-absorption, the spectra, particularly those of 12CO, often show a complex emission/velocity structure that cannot be adequately fitted by a single Gaussian. As expected, emission from the dark clouds in the sample (HH 46, the Coal Sack, the Chameleon, Lupus and Corona Australis) is distinguished by narrow line widths compared to other clouds in the sample (see figure 1).

High-velocity wings and asymmetric line profiles are seen in the spectra of a number of molecular clouds (e.g spectra of CO, CS, HCN and SO observed towards the positions NGC 6334(CO) and NGC 6334 (N)) suggesting that outflows, presumably from young stellar objects (YSOs), are present. Figure 3 shows some examples of spectra with line wings and asymmetries.

### 4.2 Isotopomer Intensity Ratios

The 12CO emission is found to be optically thick for all molecular clouds in this sample. Table 4 lists the emission intensity ratios for 12C:13C isotopomers of CO, CS, HCN, HNC and HCO+, as well as the sample standard deviation ($`\sigma _{n1}`$), the minimum and maximum values of the ratio, and the number of molecular clouds for which the ratio has been calculated. The terrestrial abundance ratio of the atoms 12C:13C is 89, the ratio of the abundances of the molecules 12CO:13CO in nearby molecular clouds is found to be around 50; however, significant variations in this ratio have been found in individual molecular clouds, and it is likely that the ratio varies with galactocentric distance. Langer and Penzias (1990) have found a value of 57 for molecular clouds at the Sun’s distance from the Galactic Centre, with values ranging from 25 near the centre of the Galaxy to 70 at a distance of 12 kpc. In contrast, the mean ratio for the intensity of 12CO emission (hereafter I12) to emission from 13CO (hereafter I13) in table 4 is 3.1 $`\pm `$ 0.2 for this sample of southern molecular clouds, implying that the 12CO emission is optically thick for all molecular clouds of the sample. This is very close to the value for I12: I13 of 3.0 $`\pm `$ 0.9 found for a sample of northern Galactic-plane GMCs between l = 34o and l = 51o (Polk et al. 1988) suggesting that there is no systematic difference in this ratio between the northern and southern Galactic plane. Doty & Neufeld (1997) have modelled I12: I13 for the dense quiescent star-forming cores of molecular clouds, assuming an approximately terrestrial abundance ratio of 100:1 for 12CO:13CO. They have predicted an I12: I13 ratio of between 2 and 6, which is consistent with the findings of this study, and again suggests that the 12CO emission has high optical depth for all clouds in this sample.

Table 5 lists the same values for the intensity ratio of 13C: 18O isotopomers. The mean ratio of I13 to emission from C18O (hereafter I18) is 6.9 $`\pm `$ 0.4. This is somewhat higher than the value for this ratio of 5.5 that would be expected for optically thin emission if terrestrial ratios for atomic 12C:13C and 16O: 18O were assumed. The values for I13: I18 for individual clouds range from 2.5 to 15.2, so in some molecular clouds the ratio is far higher than expected from the terrestrial ratio of the isotopes.

The ratio I12: I18 found in this study is 25 $`\pm `$ 3, with a range of 4.2 to 98. Doty & Neufield (1997), assuming terrestrial ratios, predicted that I12: I18 should be between 5 and 20, assuming that 12CO and C18O were coextensive. A possible reason for the high ratios found in this study is that the emission from 12CO arises from a larger volume of gas than does emission from C18O. This suggestion is supported by the average line-widths for 12CO and C18O emission for the sample. The average line-width for 12CO is 8.7 km s-1, while the average line-width for C18O is 4.1 km s-1, suggesting that 12CO emission arises from a larger volume of the gas (assuming that line-width is proportional to radius, as if will be if there is a linear velocity gradient across the molecular cloud). However, note that this difference in average line widths may be due to optical depth effects. If the 12CO emission is optically thin in the wings of the line profiles then the C18O emission in the line wings will not be detected, assuming a terrestrial abundance of 12CO:C18O of 490:1.

The average line-width of 13CO is 5.3 km s-1, similar to that for C18O, suggesting that the C18O and 13CO emission arise from a similar volume of gas. If this is the case, then the optical depth of the 12CO emission is likely to be even higher than that implied by the emission ratios in Table 4. The critical densities for 13CO and C18O are lower than that for 12CO, implying that the difference in volume of the emitting gas is due to a lower abundance of the more rare isotopomers in outer regions of the molecular cloud rather than to varying excitation conditions. The unexpectedly high values of I13:I18 can be understood if the 13CO and C18O emission are not completely co-extensive, with the C18O emission arising from a smaller volume of gas. The relative average line widths (4.1 km s-1 for C18O, 5.3 km s-1 for 13CO) provide some support for this assertion. It is possible that self-shielding (from dissociating UV radiation both within and incident on the molecular cloud) in the more abundant isotopomers accounts for their apparently wider distribution within the molecular gas (Langer & Penzias 1990).

The average ratio of emission from H13CO+ to that from HC18O+ is 15.3 $`\pm `$ 3.5 (Table 5). Once again, this is far higher than would be expected from optically thin emission from both isotopomers if terrestrial ratios of atomic 12C: 13C and 16O: 18O apply. The high ratio suggests that the HC18O+ emission may arise from a smaller volume of gas than the H13CO+ emission, once again possibly due to self-shielding in the more abundant isotopomer (H13CO+). The higher ratio for emission from H13CO+:HC18O+ compared with that from 13CO:C18O may also support this suggestion. HC18O+ can be expected to be far less abundant than C18O, and consequently full self-shielding from dissociating UV radiation will occur at a significantly smaller cloud radius. On the other hand, the high ratio may be due to chemical fractionation effects. Chemical fractionation (the preference for a molecule consisting of one isotope over another in a chemical reaction) can lead to different isotope ratios in the molecular products of a reaction with respect to the reactants. More complex products have a higher ratio of the lighter to the heavier isotope (e.g a higher 16O: 18O ratio) (Engel & Macko 1997), so HC18O+ may be less abundant with respect to its 13C isotopomer than C18O.

The average emission ratios of the 12C to 13C isotopomers of CS, HCN, HNC and HCO+ (Table 4) suggest that the 12C isotopomer of these molecules is generally optically thick in dense molecular clouds, although less so than 12CO. A surprising result is that the minimum value for the emission ratio HCN:H13CN is 0.9, and the minimum value for HCO+:HC13O+ is 0.3. Both these low values occur in the same molecular cloud, the Chameleon dark cloud. No other molecular clouds have values less than unity for these ratios. In fact, all ratios for emission from 12C isotopomers to 13C isotopomers in the Chameleon dark cloud are low compared to other molecular clouds. The reason for this is not clear. It is possible that the Chameleon dark cloud has a low elemental ratio of 12C to 13C, but further observations of the Chameleon dark cloud will be required to clarify the situation.

## 5 Summary

Transitions of 11 different molecules and several of their isotopomers have been observed in a sample of 26 southern molecular clouds. The observational results are as follows:

1. The CO spectra often show a significantly wider range of velocities than found for other molecules.

2. Self-absorption is present in about a third of the CO profiles, and in some profiles of 13CO, C18O, CS, HCN, HCO+ and SO.

3. Self-absorption is not generally observed in transitions of OCS, HC3N and CH3OH, suggesting that these transitions are optically thin in most molecular clouds of the sample, or that these transitions only exist in the higher-excitation, inner parts of the clouds.

4. CO is optically thick in all molecular clouds of the sample. The average emission ratio for 12CO:13CO is $`3.1\pm 0.2`$, very close to values found in surveys of northern hemisphere molecular clouds, and predicted by models.

5. The average emission ratio for 13CO:C18O is $`6.9\pm 0.4`$, higher than would be expected for optically thin emission if Solar System ratios for 13C:18O were assumed. This may be because C18O occupies a smaller volume of gas than does 13CO, due to self-shielding of 13CO.

6. The average emission ratio of H13CO+:HC18O+ is $`15.3\pm 3.5`$, also higher than would be expected for optically thin emission if Solar System ratios for 13C:18O were assumed. Self-shielding in H13CO+ is a possible explanation, however chemical fractionation effects may also be responsible for the high ratio.

7. CS, HCN, HNC and HCO+ are generally found to be optically thick for the molecular clouds of the sample.

8. The Chameleon molecular cloud has 12C:13C isotopomer ratios of less than unity for HCN and HCO+, and has low 12C:13C isotopomer ratios for all molecules. Further observations will be required to clarify the reasons for this situation.

## Acknowledgments

This research forms part of the PhD thesis of MRHC undertaken at the University of Western Sydney. MRHC acknowledges support from an Australian Postgraduate Award and thanks DIST Access to Major Research Facilities Scheme for travel support to SEST. MRHC also gratefully acknowledges support from the ATNF, including free accommodation and use of computing facilities for the duration of this work. The Mopra telescope is part of the Australia Telescope which is funded by the Commonwealth of Australia for operation as a National Facility managed by CSIRO.

## Appendix: Calibration of the Mopra Telescope

The accurate intensity calibration of the observations reported here is important, as they have been used to calculate abundances (Hunt Cunningham, in preparation). At millimetre wavelengths there is significant emission and absorption of radiation by the atmosphere, and the component due to water vapour in the atmosphere can vary over timescales as short as a few minutes. One way to correct for these optical depth fluctuations is the ‘chopper wheel’ method, where an ambient temperature blackbody source is periodically inserted into the antenna beam as a reference (see e.g. Rohlfs & Wilson 1996; Kutner & Ulich 1981). This method as applied to both the Mopra telescope and the SEST is outlined here. More information on the terms and definitions used below can be found in Rohlfs & Wilson (1996).

The antenna temperature directly measured at the telescope ($`T_A`$) needs to be corrected by various factors to account for atmospheric attenuation, aperture efficiency, far-sidelobe contributions (back scattered radiation) and electronic losses and noise (receiver temperature). A quantity known as the calibrated antenna temperature ($`T_A^{}`$) is corrected for atmospheric attenuation, backscattered radiation and receiver temperature, while what is known as the corrected antenna temperature ($`T_A^{}`$) is also corrected for the forward efficiency of the telescope ($`F_{eff}`$), so $`T_A^{}=T_A^{}/F_{eff}.`$

The forward beam includes the sidelobes in the forward direction, so a better parameter to use is the main beam brightness temperature ($`T_{MB}`$), where

$$T_{MB}=T_A^{}/B_{eff}=T_A^{}F_{eff}/B_{eff}.$$

The SEST observations were made with on-line chopper wheel calibration system, and the observations were corrected by the on-line software to the corrected antenna temperature $`T_A^{}`$, so only the (frequency dependent) correction $`F_{eff}/B_{eff}`$ needed to be applied to get $`T_{MB}`$.

The Mopra system did not have such on-line chopper wheel or software, so a ‘blackbody paddle’ was inserted every few minutes to measure the gated total power (GTP) which could then be used off-line to calibrate the antenna temperature $`T_A`$ for atmospheric attenuation and instrumental parameters.

Let $`T_R`$ be the receiver temperature, $`T_S`$ the source temperature, $`T_a`$ the ambient temperature of the atmosphere and $`\tau `$ the optical depth of the line of sight through the atmosphere ($`\tau =\tau _0\mathrm{sec}(Z)`$, where $`Z`$ is the zenith angle and $`\tau _0`$ the optical depth vertically through the atmosphere). The on-source signal $`S`$ is then

$$S=C[T_R+T_a(1\mathrm{exp}(\tau ))+T_S\mathrm{exp}(\tau )],$$

where $`C`$ is a constant, and the off-source reference $`R`$ is

$$R=C[T_R+T_a(1\mathrm{exp}(\tau ))].$$

The usual position switching, subtracting the reference from the signal, gives

$$SR=CT_S\mathrm{exp}(\tau )$$

correcting for the receiver temperature and the atmosphere emission, but not for the atmospheric absorption of the source emission.

The paddle reference can be considered as a black body, at the ambient temperature $`T_a`$. The power $`P_r`$ at the off-source position is

$$P_r=C^{}[T_R+T_a(1\mathrm{exp}(\tau ))]$$

where $`C^{}`$ is another constant, and the power $`P_p`$ with the black-body paddle in the telescope beam is

$$P_p=C^{}[T_R+T_a]\text{ so }P_pP_r=C^{}T_a\mathrm{exp}(\tau ).$$

This can be used to correct the observations for atmospheric attenuation by

$$(SR)/(P_pP_r)=(C/C^{})(T_S/T_a)\text{ or }$$

$$T_S=(C^{}/C)T_a(SR)/(P_pP_r)$$

where $`C`$ and $`C^{}`$ are constants of the telescope hardware and software to be determined and calibrated out.

A complication of the Mopra spectral calibration here is that the observing program SPECTRA did not output the difference $`SR`$ of the signal and reference spectra, but the quotient

$$Q=K[(S/R)1]T_R=KT_R(SR)/R$$

where $`K`$ is a constant required from the instrumental calibration.

The standard calibration method equation is modified to

$$QP_r/(P_pP_r)=KT_R(SR)P_r/[R(P_pP_r)]$$

$$=KT_RT_S/T_a\text{ or }$$

$$T_S=Q(T_a/KT_R)P_r/(P_pP_r).$$

The Mopra telescope was not calibrated absolutely, for the period that these observations were made, so the factor $`K`$ was significant, and we had to do a relative calibration by comparing Mopra and SEST observations of the calibrator sources Orion MC and M17. The integrated line intensities (in K km s-1) were used for 43 lines between 86 and 115 GHz observed at Mopra several times over a range of weather conditions. Before calibration, the scatter of repeated observations in $`Q`$ about the mean for that line was rms 17 % with some data points more than 30 % from the mean, indicating the effect of variable atmospheric absorption. After calibration, the scatter about the mean was 7 %, showing the effectiveness of the correction for absorption. Figure 4 shows the uncalibrated Mopra integrated antenna temperatures ($`Q𝑑v`$) plotted with the SEST integrated, corrected antenna temperatures ($`T_A^{}𝑑v`$), while figure 5 shows the same data after the calibration procedure outlined above.

The Mopra spectra quotient $`Q`$ values were converted to the SEST corrected antenna temperature scale $`T_A^{}`$ with a factor $`1.62\pm 0.06`$ accounting for the instrumental factors in the Mopra SPECTRA $`Q`$ scale.


format: text

text:
# Low-mass supernovae in the early Galactic halo: source of the double 𝑟/𝑠-process enriched halo stars?

## 1 Chemical peculiarities in halo stars

The metallicity distribution of the stars in the Galactic halo ranges from \[Fe/H\]$`=1`$ for the most metal-rich stars, to \[Fe/H\]$`=5.3`$ at the minimum (HE 0107$``$5240: Christlieb et al. (2002)). The first halo stars formed from almost unprocessed primordial gas. A short formation time is favoured, because chemical enrichment due to supernovae is rapid. The Milky Way disk formed out of gas pre-enriched to about \[Fe/H\]$`1`$. The pre-enrichment may have been due to the Bulge formation (Renzini, 2003) but the most metal-rich stars in the halo show that this level of enrichment was already reached by the end of the halo formation. The global metallicity of the Universe reached \[Fe/H\]$`1`$ at $`z=3`$ (Renzini, 2003).

Some halo main-sequence stars show evidence for abundance alterations. Because main-sequence stars should not have significantly altered their surface abundances, and pollution from planetary companions is unlikely given the dearth of planets at subsolar metallicity (Bodaghee et al., 2003), pollution from a companion star by mass transfer is considered most likely. In some cases high overabundances of s-proces elements (e.g. Pb) are found (Aoki et al., 2000; Van Eck et al., 2001; Lucatello et al., 2003). These elements are formed in Asymptotic Giant branch (AGB) stars (Van Eck et al., 2003), which first appeared $`10^8`$ yr after the halo formation. Other metal-poor stars show evidence for large enhancements of r-process elements, a.o. europium (Sneden et al., 2000) and gold (Cowan et al., 2002). These form in type-II supernovae, and are also explained via pollution from a companion star.

The discovery that several stars show enhancements of both r-process and s-process elements (Hill et al., 2000; Cohen et al., 2003) is puzzling, as they require pollution from both an AGB star and a supernova.

In this paper we will present a calculation which indicates that at low metallicity, inefficient mass loss may allow the degenerate cores of AGB stars to reach the Chandresekhar mass. This would lead to a AGB supernovae, polluting low-mass companions by both s-process elements (enriched in the AGB envelope) and r-process elements formed in the subsequent explosion. The next section discusses the double-polluted stars. Section 3 discusses AGB mass loss at low $`Z`$, and Section 4 derives initial-final masses for these stars and discusses the possibility of AGB supernovae.

## 2 Double polluted halo stars

Three stars are known to show the double enhancements: HE2148$``$1247 (Cohen et al., 2003), CS22948$``$027 (Hill et al., 2000; Preston & Sneden, 2001) and CS29497$``$034 (Hill et al., 2000), with \[Fe/H\]=$`2.3`$,$`2.45`$, and $`2.90`$ respectively. CS22898$``$027 (Preston & Sneden, 2001; Aoki et al., 2002), at \[Fe/H\]$`=2.25`$, may also show double enrichment. The metallicities are at the lower range of the gaussian halo distribution (Ryan & Norris, 1991). The process leading to the double enhancement appears to require \[Fe/H\]$`<2`$. A lower limit to the metallicities of the double enriched stars is not evident. Although there is a tail in the halo distribution extending to \[Fe/H\]$`=4`$ (Ryan & Norris, 1991) or lower (Christlieb et al., 2002), there are few stars in this tail.

Cohen et al. (2003) discuss possible scenarios for double enhancements. A triple system including one massive pre-supernova and one pre-AGB star, plus the remaining low-mass star, appears very unlikely. They therefore suggest reverse mass-transfer from the low-mass star to the post-AGB white dwarf, leading to accretion-induced collapse, and a type-Ia supernova. A third star in a close orbit around the AGB primary may instead be considered for the mass donor, leading to a more traditional type-Ia scenario, but this again leads to an unlikely close triple system. There is also no evidence that type-Ia SNe manufacture significant amounts of r-process elements, arguing against both these scenarios.

Pollution from an AGB companion (s-process elements) occurs via wind accretion. The peak transfer efficiency ocurs for binaries with periods around 3000 days (Pols et al., 2003). Han et al. (1995) stress the additional importance of a common envelope phase in mass transfer, but orbital eccentricities show that all but the closest binaries ($`d<1\mathrm{AU}`$) avoid a common envelope (Karakas et al., 2000; Pols et al., 2003). Roche lobe overflow is also unlikely (Han et al., 2002), firstly because it tends to happen on the first giant branch (before the onset of the s-process enhancements) and secondly because it is unstable for donor stars more massive than the secondary, as would have been the case here. For the more likely wider orbits, the wind interaction has little or no effect on the primary AGB star, which continues to evolve as in the case of a single star.

The wide orbit argues against the accretion-induced collapse, because of the need to transfer a large amount of mass from a distant low-mass main-sequence star in order to reach the Chandresekhar mass. A scenario which does not require either a triple system or reverse mass transfer, and explains the metallicity dependence of the double-enrichment process, could be of interest.

## 3 AGB mass loss at low metallicity

### 3.1 Mass loss expectations

Stars on the AGB burn helium and hydrogen in shells around their inert C/O core. Mass loss increases during the AGB: once the mass-loss rate in the wind significantly exceeds the nuclear burning rate ($`10^7\mathrm{M}_{}\mathrm{yr}^1`$), the evolution comes to a sudden end (Willson, 2000). The remaining envelope is quickly removed and nuclear burning will cease. The C/O core at this point forms the subsequent white dwarf.

The dependency of the mass loss on the stellar parameters is not well known. Several formalisms have been proposed in the literature, depending typically on the radius, mass and luminosity of the star. The one used most extensively is the Bloecker (1995) relation

$$\dot{M}=\mathrm{4.8\hspace{0.17em}10}^9\frac{LR}{M}\frac{L^{2.7}}{M}\mathrm{M}_{}\mathrm{yr}^1,$$
(1)
with quantities in solar units. Such relations are generally derived for Galactic stars. An explicit metallicity dependence is not included.

The mass loss of AGB stars is a two-step process. First, pulsations extend the atmosphere and drive a small mass loss. Secondly, dust forms in the extended atmosphere: radiation pressure on the dust now drives the large mass-loss rates observed. At low metallicity, the efficiency of the dust formation is reduced and this will limit the mass-loss rates. Bowen & Willson (1991) derive mass-loss rates at low $`Z`$ from theoretical pulsation models. They find that at low metallicity (\[Fe/H\]$`<1`$) the dust does not play a significant role and the wind becomes purely pulsation driven. Stars with lower metallicity are found to have lower mass-loss rates. This will allow the core to grow to a larger mass and reach higher final masses.

The AGB wind depends on the stellar radius. As low-metallicity stars have smaller radii (Iben (1984) gives for AGB stars the relation $`RZ^{0.088}`$), this implies a further metallicity dependence of the mass loss, even for dust-free, pulsation driven winds.

### 3.2 Observational Evidence

There are few direct observations of mass loss at low metallicity. The best data set comes from ISO observations of AGB stars in the Magellanic Clouds (Trams et al., 1999), from which both mass-loss rates and luminosities were derived (van Loon et al., 1999). Comparison with Galactic stars is not straightforward, because within the Galaxy distances tend to be poorly known and the luminosities are uncertain. The ISO data shows that both LMC and SMC stars reach similar mass-loss rates to those shown by Galactic stars near the tip of the AGB (van Loon, 2000), although the unknown dust-to-gas ratios and expansion velocities could hide some difference. The luminosity of the LMC and SMC stars tends to be high, reaching up to $`L=5\times 10^4\mathrm{L}_{}`$.

Mass loss at even lower metallicities (\[Fe/H\]$`<0.7`$) is very poorly studied. However, the superwind phase will lead to the formation of a planetary nebula (PN) and these can be observed to very large distances. There is a good relation between the number of PNe and the luminosity of the host galaxy: for the most populated systems (Bulge, LMC), $`\mathrm{log}[N(\mathrm{PN})]\mathrm{log}L_V/L_{}6.9`$ (Magrini et al., 2002). However, local galaxies with low $`Z`$ show some evidence for a deficiency in the number of PN (Magrini et al., 2003). This becomes apparent only for \[Fe/H\]$`<1`$ (Fig. 1), well below the LMC or SMC metallicity. Plotting the deficiency as function of ratio of C-type over M-type AGB stars (Fig. 1, right panel) shows a clear relation: this ratio is an effective tracer of the metallicity of the AGB population (Groenewegen, 1999).

The deficiency of PNe suggests that the mass-loss rates do not reach as high values for metallicities \[Fe/H\]$`\begin{array}{c}<\hfill \\ \hfill \end{array}1`$. Lower peak mass-loss rates will give rise to less dense PNe which will fade faster. Conversely, the lack of a relation at \[Fe/H\]$`>1`$ suggests that all these AGB stars reach similar mass-loss rates (although not necessarily at the same luminosities) regardless of metallicity. This change of behaviour at \[Fe/H\]$`=1`$ is in excellent agreement with the predictions of Bowen & Willson (1991).

## 4 AGB supernovae

We can estimate the effect of metallicity on the initial–final-mass. The goal is to obtain qualitative estimates: approximate, simple relations will suffice. In a more detailed calculation, more complicated relations could be substituted but these introduce their own uncertainties. The largest uncertainties comes from the mass-loss estimates, and not including hot bottom burning.

The relation between the luminosity corresponding to a certain mass-loss rate, and the metallicity, is taken from Bowen & Willson (1991). Between $`Z=Z_{}`$ and $`Z=0.1\times Z_{}`$, this luminosity increases by a factor of 1.3, in part because the wind is no longer dust driven. At lower Z, the luminosity required to support the same mass-loss rate scales as (Bowen & Willson, 1991) :

$$\mathrm{log}L(Z)/L(Z_{})=0.120.13\mathrm{log}Z/Z_{};Z/Z_{}0.1.$$
(2)
The effect on the final mass of the star is derived from the core mass–luminosity relation, which we take from Boothroyd & Sackmann (1988):

$$\frac{L}{L_{}}=59250\left(\frac{M_\mathrm{c}}{M_{}}0.495\right)$$
(3)
The increase of the final mass is very small for low core mass, but becomes significant for core masses corresponding to the heaviest white dwarfs. We use an approximate present initial–final mass relation corresponding to solar metallicity:

$$M_\mathrm{f}=0.5+\frac{1}{12}M_\mathrm{i}.$$
(4)
The true relation is shallow at low core mass and steepens at high core mass (Vassiliadis & Wood, 1993), so this is a significant simplification which may be more realistic at high than at low $`M_\mathrm{c}`$.

The result is shown in Fig. 2. At solar metallicity, the most massive AGB stars reach final masses around 1 M. Already at $`Z/Z_{}=0.1`$, the most massive remnants approach the Chandresekhar limit. \[This should be taken with caution as the most massive AGB stars may show higher luminosity than predicted from their core masses, due to the effects of hot bottom burning (Bloecker & Schoenberner, 1991).\] At $`Z/Z_{}<0.01`$, the Chandresekhar mass is reached for 5 M stars and at $`Z/Z_{}<0.001`$, for 4 M.

Carbon ignition in the degenerate core will cause these stars to explode. Iben & Renzini (1983) named such objects ’type 1.5 supernovae’; physically they could be considered an extension of type II SNe. The early work did not include the superwind which was later found to terminate the AGB before these events could occur. However, an absence of a superwind at the lowest metallicities means that AGB supernovae could still play a role in the early evolution of the Universe.

## 5 Discussion

### 5.1 Double enriched stars

We have presented calculations which indicate the possibility that massive AGB stars in the early halo would have given highly massive remnants, leading to a population of AGB supernovae. In the context of the mixed r/s-enhanced halo stars, this suggests that the same companion star was responsible for the double pollution. For an initial mass in excess of 3–4 M, the companion would first have developed into an AGB star. The wind of this star would have carried the s-process elements. After the core reached the Chandresekhar limit, a supernova explosion occured. r-Process elements formed during the explosion also polluted the companion star. The explosion may disrupt the binary system and so this model does not automatically predict that all such stars are binaries. However, if a binary, the companion should be a supernova remnant with appropriate mass.

Compared to the previously poposed models for the double enriched stars, the present model has the advantage of not requiring a close triple but only a binary, and it includes an explanation for the low metallicities of the doubly enriched stars. r-Process elements appear to be produced in low-mass (10 M) type-II supernovae (Pagel & Tautvaisiene, 1995): the similar masses suggest their production could also be expected in AGB supernovae.

### 5.2 Lithium

Preston & Sneden (2001) find that CS 22898$``$27 has a more or less normal lithium abundance for its metallicity, and wonder whether this is easily explained assuming pollution from an AGB star. AGB stars experience lithiu-rich phases during hot bottom burning. However, the lithium is easily destroyed again and at different phases, the surface abundances can vary from very high to strongly depleted (Sackmann & Boothroyd, 1992): the lithium-rich phase is short-lived. The total lithium yield of low metallicity, hot bottom burning stars can be very close to the original Spite plateau abundance (Ventura et al., 2002). Thus, the normal lithium abundance can be a coincidence.

For hotter s-process enriched stars, the lithium 6708Å line coincides with a CeII line and derived abundances may be upper liimits in such cases (Reyniers et al., 2002).

### 5.3 Supernova rates

The effect on the supernova rate could be significant. Assuming a Salpeter IMF with $`\alpha =2.35`$, the ratio of stars exploding as SN on the AGB, compared to the normal type-II’s is shown in Fig. 2. At the extreme $`Z`$, the contribution of these AGB-SN exceeds the normal type-II SN by a factor of a few. The AGB stars will however explode much later and would therefore play less of a role in the fast enrichment of the gas.

### 5.4 Globular clusters

Old low-$`Z`$ population should show the effect of the higher final masses, and contain a population of heavy white dwarfs and an overpopulation of neutron stars. This would affect globular clusters: dynamical evidence for the presence of a high-$`M/L`$ population of heavy remnants in the core of M15 (\[Fe/H$`=2.25`$) is presented by Phinney (1992), Dull et al. (1997) and Baumgardt et al. (2003). A high $`M/L`$ was also found in the core of NGC 6752 (\[Fe/H$`=1.56`$), which could be caused by an excess population of heavy white dwarfs or neutron stars (Ferraro et al., 2003).

## 6 Conclusions

From our present knowledge of AGB mass loss, the early Universe may have contained a population of supernovae which is not found in the local, high metallicity Universe. These AGB-SN are possible sites for the r-process, and provide a simple explanation for the small sample of low-metallicity, Galactic halo stars which show evidence for both s-process and r-process enrichment.

The absence of dust driven winds at \[Fe/H\]$`=1`$ is supported both by models (Bowen & Willson, 1991) and observations (this paper). The prediction that stellar remnants are more massive at lower metallicities appears robust, Regardless of whether these remnants reach the Chandresekhar limit. The presence of high-mass white dwarfs in globular clusters is likely, and could explain the evidence for an excess $`M/L`$ in the cores of both M 15 and NGC 6752.


text:
# Sub-millimeter emission from type Ia supernova host galaxies at 𝑧=0.5

## 1 Introduction

In recent years, deep surveys have found galaxies and QSOs at redshifts greater than six, making it possible to trace the cosmic history of star formation and AGN activity from the epoch of reionization to the present day. In practice however, achieving these goals is difficult, as translating observations of light from galaxies into measures of stellar content and formation histories is affected by the obscuring effects of dust. Whilst dust obscuration is known to have a significant effect locally (e.g. Tresse & Maddox (1998)), the role of dust at higher redshifts is even more important. Chemical evolution models (Calzetti & Heckman, 1999; Pei et al, 1999) predict that dust obscuration should peak at 2-3 times the current value at $`1<z<2`$, before slowly declining. Observationally, the importance of dust obscuration at high redshifts is exemplified by the very different cosmic star formation histories inferred from optical and infrared surveys. Infrared surveys produce star formation histories that rise more sharply from $`z=0`$ to $`z=1`$, and show a slower decline at higher redshifts, than the star formation histories inferred from optical surveys (Lilly et al, 1996; Madau et al, 1996; Steidel et al, 1999; Rowan-Robinson et al, 1997; Hughes et al, 1998).

These surveys however do not examine extinction evolution in inactive galaxies. Blank-field infrared surveys are biased towards systems containing a dusty starburst or AGN, whereas Lyman Break galaxies are selected on rest-frame UV flux, which may bias them towards having high rates of unobscured star formation. In this paper, we present preliminary results from a program to use deep sub-millimetre observations to investigate the evolution of sub-mm emission and dust opacity with redshift in what should be normal galaxies; the host galaxies of distant type Ia supernovae. Observations and data analysis are described in §2, and results are presented in §3. We discuss mean extinction levels and dust properties in §4, and also mention some implications for deep sub-millimetre surveys and supernova cosmology. We assume $`\mathrm{\Omega }_0=1`$, $`\mathrm{\Omega }_m=0.3`$, $`\mathrm{\Lambda }=0.7`$, and $`H_0=70`$.

## 2 Observations

The hosts of the type Ia supernovae discovered by the Supernova Cosmology Project (SCP, Perlmutter et al (1999)) and the High Redshift Supernova Search Team (HZT, Riess et al (1998)) are ideal targets for studying the evolution of dust opacity in normal galaxies; the galaxies have spectroscopic redshifts, the galaxies are selected independent of their apparent magnitudes as the discovery is based on the supernova rather than the host, and type Ia supernovae come from evolved stellar populations. Furthermore, since the supernova spectrum must be disentangled from the host galaxy spectrum, none of the host galaxies contain the signature bright emission lines of a starburst or AGN (Riess et al (1998); Perlmutter et al (1999), A Riess priv. comm., G Aldering priv. comm.).

We observed seventeen host galaxies of supernovae discovered by the SCP and the HZT. The only selection criteria were that the host galaxy redshifts should lie in a narrow interval around z=0.5, so that the fluxes could be coadded without $`k`$-correcting to a common wavelength, and that the supernova in the host should be spectroscopically confirmed as a type Ia. The host galaxies were observed at $`450\mu `$m and $`850\mu `$m using the Sub-mm Common User Bolometer Array (SCUBA) on the James Clerk Maxwell Telescope (JCMT). All observations were performed in photometry mode. The atmospheric conditions were of exceptional quality, with a mean $`850\mu `$m sky opacity of $`\tau _{850}0.15`$. We note the possibility that sub-mm flux from any of the objects in our sample could in principle arise from a chance-aligned high redshift source, however assuming a sub-mm source density from recent blank-field sub-mm surveys (Eales et al, 1999; Scott et al, 2002; Fox et al, 2002; Borys et al, 2003) and a conservative search radius of $`2\mathrm{}`$ (the jiggle offset on photometry mode observations with SCUBA) gives a probability that the flux from any one object in the sample comes from a chance-aligned sub-mm source of approximately 1 in 400. We therefore conclude that any sub-mm flux seen in our sample is most likely due to the target, rather than a background source. The data were reduced using the SCUBA User Reduction Facility (SURF) pipeline using standard methods.

For one host galaxy, that of SN1997ey, we obtained Hubble Space Telescope (HST) imaging from the HST data archive. Observations were taken using the Space Telescope Imaging Spectrograph (STIS) using a clear filter. The data were reduced using the IRAF reduction package calstis. The final image is presented in Figure 1.

## 3 Results

The sample, their redshifts, and the sub-mm fluxes are presented in Table 1. One object (the host galaxy of SN 1997ey) is detected strongly at both $`450\mu `$m and $`850\mu `$m, one object (the host of SN2000eh) is detected at $`850\mu `$m only, the rest of the sample are undetected individually.

As our sample is randomly selected over a narrow redshift range, we can coadd the observed sub-mm fluxes to measure the mean sub-mm flux of normal galaxies at $`z0.5`$. We only consider the $`850\mu `$m fluxes, as the $`450\mu `$m calibration errors are large. The mean observed $`850\mu `$m flux for all 17 galaxies is 1.55mJy$`\pm `$0.31mJy. The host galaxy of SN1997ey is however significantly brighter than the others, and it is probable that this source is an outlier, and that the flux of this galaxy will bias the mean value significantly. Assuming that the ‘true’ source flux distribution is a Gaussian, then the distribution of fluxes for the whole sample is consistent with a gaussian at only the $`5\sigma `$ level, whereas excluding the host of SN1997ey makes the distribution of fluxes consistent with a gaussian at $`<2\sigma `$, and excluding any of the remaining sources does not make the fit statistically better. We therefore also quote the error weighted mean flux of the 16 galaxies with the host of SN1997ey excluded, which is 1.01mJy$`\pm 0.33`$mJy, and use this value in the following analysis.

## 4 Discussion

We first review the origin of sub-mm flux from extragalactic sources, and the relation between sub-mm flux and optical extinction. In many cases sub-mm emission from extragalactic sources arises in starbursts; compact regions ($`<1`$kpc across) containing $`10^8M_{}`$ of dust, heated by young massive stars forming at a rapid rate ($`>50M_{}`$yr<sup>-1</sup>). Since sub-mm emission measures physical dust volume rather than the volume of space the dust is dispersed in, the same level of sub-mm emission can arise from a comparable mass of dust distributed in the disk of a galaxy, heated by the general interstellar radiation field. This type of sub-mm emission, known as ’cirrus’ emission, is seen locally, both in Galactic sources (Casey, 1991), and in nearby galaxies (Bianchi, Davies & Alton, 1999; Dunne et al, 2000). Furthermore, there is evidence that cirrus emission may be important in blank-field sub-mm sources (Efstathiou & Rowan-Robinson, 2003). As our targets contain a type Ia supernova, known to come from evolved stellar populations, and because none of the targets contain optical signatures of starbursts or AGN, we are confident that any sub-mm flux from our targets will be cirrus emission.

The relation between cirrus sub-mm emission and optical extinction has been discussed by Hildebrand (1983), Casey (1991) and Bianchi, Davies & Alton (1999). It turns out that the integrated sub-mm flux from a galaxy, with certain assumptions, scales linearly with $`A_V`$. We summarize these arguments here. If the dust responsible for the sub-mm emission is also responsible for the optical extinction, then the ratio of the extinction efficiency in the V band to the emission efficiency in the sub-mm at some wavelength $`\lambda `$ is equal to the ratio of the optical depths:

$$\frac{Q_{ext,V}}{Q_{em,\lambda }}=\frac{\tau _V}{\tau _\lambda }$$
(1)
In the optically thin case, $`\tau _\lambda `$ can be written:

$$\tau _\lambda =\frac{I_\lambda }{B_\lambda (T)}$$
(2)
where $`I_\lambda `$ is the sub-mm intensity and $`T`$ is the dust temperature. Solving for $`A_V`$ gives (Bianchi et al, 1998):

$$A_V\frac{Q_{UV}}{Q_{\lambda _0}}\left(\frac{\lambda }{\lambda _0}\right)^\beta \frac{I_\lambda }{B_\lambda (T)}$$
(3)
where $`\lambda `$ is the observed wavelength, $`\lambda _0`$ is a reference wavelength, and $`Q_{\lambda _0}`$ is the extinction efficiency at the reference wavelength.

There are two important assumptions in this derivation, which may affect our results. Firstly, even though the fluxes lie within the Rayleigh-Jeans tail, the sub-mm flux still scales linearly with dust temperature, meaning that a mean increase in temperature could be responsible for increasing the sub-mm flux, rather than an increase in optical extinction. However, even assuming pure luminosity evolution between the local Universe and $`z=0.5`$, the mean dust temperature will only increase as the 1/5th power of $`A_V`$ (i.e. increasing with the mean surface brightness). Furthermore, the predicted rise in dust temperature between z=0 and z=0.5 from chemical evolution models (Pei et al, 1999) is probably of the order 1K or less, insufficient to produce a significant rise in sub-mm flux. Secondly is the gas-to-dust ratio. When using this formula to compare optical extinctions between different objects, it is assumed that the gas-to-dust ratio for both objects is the same. When comparing local galaxies to galaxies at $`z=0.5`$ this criterion will not be satisfied; globally we expect significantly less star formation to have occurred by $`z=0.5`$, which may lead to a higher mean gas-to dust ratio than locally. Based on indirect evidence (Pettini et al, 1999), we expect this rise to be no more than $`10\%20\%`$ at $`z=0.5`$, although at higher redshifts the change may be larger (Pei et al, 1991). As there is no published measure of the gas-to dust ratio in normal galaxies at $`z=0.5`$, we have assumed that the change from the local mean value is negligible.

### 4.1 Dust opacity evolution in normal galaxies

To establish if the mean sub-mm flux of our sample differs from the mean sub-mm flux of local normal (i.e. inactive) galaxies, we would like to compare to a local sub-mm survey of galaxies selected on optical flux. Currently however, no such survey has been published. Therefore, to estimate the expected observed mean $`850\mu `$m flux at $`z=0.5`$ from normal galaxies if there had been no evolution in dust opacity between the local Universe and $`z=0.5`$, we consider the local sub-millimetre survey of Dunne et al (2000). We have redshifted the mean $`850\mu `$m flux from this survey to z=0.5, with $`k`$-corrections derived using the cirrus models from Efstathiou & Rowan-Robinson (2003). The resulting predicted observed $`850\mu `$m flux at $`z=0.5`$ is 0.8mJy $`\pm 0.1`$mJy. This survey is however of galaxies from the IRAS Bright Galaxy Sample, which selects against galaxies with small dust masses, and is composed of spiral, irregular and interacting systems. The morphologies of the galaxies discovered by the supernova cosmology teams are however known to be $`70\%`$ spirals and irregulars, and $`30\%`$ ellipticals (Farrah et al, 2002b; Sullivan et al, 2003). As our sample is randomly selected from the SNeIa discovered by both the SCP and the HZT, it is reasonable to assume that the same morphological mix is present in our sample. In order to compare our value and the ‘no evolution’ value, we must correct the ‘no evolution’ value to reflect the morphological mix in our sample. If we assume that sub-mm emission from ellipticals is negligible, then this gives a predicted ‘no evolution’ flux estimate at $`z=0.5`$ of 0.56mJy $`\pm 0.1`$mJy. Comparing our observed mean $`850\mu `$m flux of 1.01mJy $`\pm 0.33`$mJy (i.e. excluding the host of SN1997ey) with this value, implies a rise in sub-mm flux from normal galaxies at $`z=0.5`$ compared to locally in the range $`25\%135\%`$ ($`1\sigma `$ errors). Therefore, our results imply that sub-mm emission from normal galaxies at $`z=0.5`$ is comparable to or higher than sub-mm emission from local normal galaxies.

We can use this mean sub-mm flux to estimate the mean value of $`A_V`$ in spiral and irregular galaxies at $`z=0.5`$. The most recent estimate of mean extinction levels in local normal galaxies (i.e including only an optically thin cirrus component) is that of Rowan-Robinson (2003), who derives $`E(BV)=0.1`$, and therefore $`A_{V,local}=0.31`$. Hence, we obtain $`A_{V,z=0.5}=0.56\pm 0.17`$. This result is in good agreement with indirect observational evidence based on the star formation history of the Universe (Lilly et al, 1996; Rowan-Robinson et al, 1997). It is also in agreement with model predictions for the evolution in dust opacity with redshift (Calzetti & Heckman, 1999), which predict that the change in mean reddening between the local Universe and z=0.5 will be $`\mathrm{\Delta }E(BV)=0.05`$, with a corresponding increase in extinction of $`A_V0.15`$, and with chemical evolution models (Pei et al, 1999) which model the evolution of gas and dust mass with redshift.

### 4.2 Dust properties in spiral galaxies at z=0.5

We now examine the dust properties in the host galaxy of SN1997ey, which was detected at $`450\mu `$m and $`850\mu `$m. One possibility is that this galaxy harbours a dust enshrouded starburst that produces the sub-mm flux; such a system would be unsuitable for examining the properties of dust in normal galaxies at $`z=0.5`$. This is however very unlikely, for two reasons. Firstly, as described earlier, none of the optical spectra of the supernova host galaxies from either the SCP or the HZT, from which this galaxy is selected, show signs of an AGN or starburst. Secondly, the HST STIS image for this galaxy, taken after the supernova had faded and presented in Figure 1, shows a disk galaxy with no signs of morphological disturbance that are seen in local (Surace et al, 1998; Farrah et al, 2001; Bushouse et al, 2002) or high redshift (Farrah et al, 2002a) starburst galaxies. We therefore conclude that the observed sub-mm flux comes from the host galaxy of SN 1997ey, and that this galaxy is a normal disk galaxy.

As there exists sub-mm data at only two wavelengths, the properties of the dust in this galaxy are best estimated using a simple greybody formula. Figure 2 shows a reduced $`\chi ^2`$ distribution plot of $`\beta `$ and $`T`$ for the dust in the galaxy. Two particularly interesting cases are for $`\beta 2`$, the cirrus or ‘Milky-Way’ type dust solution, and for $`\beta 0`$, the ’grey’ dust solution. Whilst ‘grey’ dust has previously been proposed to exist at high redshifts (Aguirre, 1999), it is notable that, for $`\beta 2`$, the dust temperature range is comparable to that seen locally, and consistent with the models of Pei et al (1999). We therefore find no observational evidence that the dust at high redshift differs fundamentally in nature from the dust seen in the local Universe, and conclude that the dust in the host galaxy of SN1997ey is normal cirrus dust.

### 4.3 Normal galaxies and sub-millimeter surveys

The detection of one object in our sample, with an $`850\mu `$m flux of $`8`$mJy, highlights an alternative interpretation for the nature of some of the sources found in blank-field sub-mm surveys. We have observed 17 galaxies which have $`B23`$, and detected one of them at an $`850\mu `$m flux of $`8`$mJy. There are approximately 3000 galaxies per sq deg at this magnitude, so our detection rate is $`200`$ per sq deg. Whilst statistics based upon one object are obviously not trustworthy, it is interesting that this detection rate agrees surprisingly well with the detection rate reported in previous blank-field sub-mm surveys to comparable depths (Scott et al, 2002; Fox et al, 2002), who find 19 objects to 8 mJy in 260 sq arc min, or 270 per sq deg.

The nature of our source, and the assumed nature of blank-field survey sub-mm sources, are different. Blank-field survey sub-mm sources are generally thought to lie at $`z>1`$, and to be galaxies undergoing extreme bursts of dust-shrouded star formation (Fox et al, 2002). The redshift, luminosity and star formation rates are however in many cases estimated using only a limited set of photometric data. Our detection, which has a comparable flux to the sources found in these surveys, lies at $`z=0.58`$ and, judging by the coldness of the dust and morphology of the galaxy, is forming stars at a more sedate rate. Whilst some sub-mm sources do now have spectroscopic redshifts that put most of them at $`z1`$ (Lilly et al, 1999; Chapman et al, 2003), it is reasonable to infer that some of the sub-mm sources in blank-field surveys are at lower redshifts than previously thought, and also have lower star formation rates. Whilst the star formation rates in these sources will be high, they still need to be carefully accounted for in current and future sub-mm surveys, to avoid overestimating the global history of star formation.

### 4.4 Implications for supernova cosmology

There is now strong, albeit indirect, evidence from observations of the cosmic microwave background (Jaffe et al, 2001; Spergel et al, 2003), and of large-scale structure (Peacock et al, 2001) that the total density of the Universe is dominated by ’dark energy’. To determine the nature of dark energy, it is necessary to track the evolution of dark energy with redshift directly, using a cosmological standard candle. Currently, the only method deployed to do this uses the luminosity distances of distant type Ia supernovae. The HZT and the SCP have performed surveys to find high redshift supernovae, and both groups claim that dark energy dominates the total density of the Universe, and that the expansion of the Universe is currently accelerating (Riess et al, 1998; Perlmutter et al, 1999). The HZT have further claimed (Riess et al, 2001) to have found an earlier epoch of decelerating expansion, based on an apparent relative brightening of SNeIa at $`z>1`$ and that accelerating expansion commenced at $`z1`$. This relative brightening supports a real, ’cosmological’ influence on the supernova magnitudes, as alternative explanations (e.g. ‘grey’ dust, progenitor evolution) do not readily predict this relative brightening at $`z>1`$.

These results are based upon the distant supernovae being $`\mathrm{\Delta }`$m = 0.40 dimmer (with a flat Universe prior) than if the expansion of the Universe was decelerating. There are alternative interpretations that can in principle account for this dimming, one of these being host galaxy extinction. The issue of extinction towards distant type Ia supernovae has been examined in detail by both the SCP and the HZT by measuring the reddening of the supernova light curves; the SCP by comparing the rest-frame colours of their low and high redshift samples, and the HZT by fitting to the light-curves via the MLCS method. The SCP find that the reddening of both their local and distant sample, after a small number of reddened supernovae are removed, has a small scatter with a mean value of close to zero. The HZT find that only a few of their distant supernovae are significantly reddened. An analysis of the colour excesses derived from the supernova light curves as a function of host galaxy morphology by the SCP (Sullivan et al, 2003) also found only a modest difference in extinction between the SNe in ellipticals and those SNe in spirals, with $`E(BV)=0.03`$ for the SNe in the ellipticals and $`E(BV)=0.07`$ for the SNe in spirals. Our result, of a mean value of $`A_V`$ in galaxies at $`z=0.5`$ that is $`25\%135\%`$ higher than locally, does not cast any doubt on the supernova teams results, as this extinction level at $`z=0.5`$ lies within $`1.3\sigma `$ of the local value, and is derived from the mean extinction in each galaxy rather than the line of sight extinction towards the supernovae. It does however highlight the need for caution in general in using supernovae as probes of the expanding Universe, as our derived mean extinction, $`A_V=0.56\pm 0.17`$, implies a rise that is at face value comparable to the dimming ascribed to dark energy. Therefore, our result emphasizes the need to accurately monitor the extinction towards distant supernovae if they are to be used in measuring the cosmological parameters.

## 5 Summary

We have presented deep sub-mm observations of seventeen type Ia supernova host galaxies at $`z0.5`$, from which we conclude the following:

1) The mean observed-frame sub-mm flux of the sample, excluding one bright object, is $`S_{850}=1.01`$mJy$`\pm 0.33`$mJy. Assuming that sub-mm flux scales linearly with optical extinction, then this implies a rise in optical extinctions in normal, inactive galaxies compared to locally of $`25\%135\%`$. Scaling from the local value of $`A_V`$ given by Rowan-Robinson (2003) for optically thin cirrus emission gives $`A_{V,z=0.5}=0.56\pm 0.17`$. This result is in good agreement with both chemical evolution models (Pei et al, 1999; Calzetti & Heckman, 1999) and with other, indirect, observational evidence (Lilly et al, 1996; Rowan-Robinson et al, 1997). The temperature and emissivity of the dust in the brightest sub-mm object in our sample are comparable to the temperature and emissivity of dust in local galaxies.

2) The discovery of a moderate redshift disk galaxy with a sub-milliimetre flux comparable to the sources found in blank-field sub-mm surveys suggests that some of these blank-field sources may not be high redshift starbursts, but lower redshift dusty disks. These sources must be carefully accounted for in current and future sub-mm surveys, to avoid overestimating the global history of star formation.

3) Our results, when combined with previous work, infer a level of extinction in galaxies at $`z=0.5`$ that could in principle produce a dimming that is comparable in size to the dimming ascribed to accelerated expansion, although the error on the value is large and the extinction level is the mean for the galaxy rather than the line of sight extinction towards the supernovae. This emphasizes the need to carefully monitor extinction levels towards distant supernovae if they are to be used to track the expansion rate history of the Universe.

We thank Margrethe Wold, David Branch and Carol Lonsdale for helpful discussion, and the referee for a very helpful report. The JCMT is operated by the Joint Astronomy Centre on behalf of the UK Particle Physics and Astronomy Research Council, The Netherlands Organization for Scientific Research and the Canadian National Research Council. DF was supported by NASA grant NAG 5-3370 and by the Jet Propulsion Laboratory, California Institute of Technology, under contract with NASA. MF and DLC were supported by PPARC. JA acknowledges support from the Science and Technology Foundation (FCT, Portugal) through the fellowship BPD-5535-2001.


text:
# Program summary

## Program summary

Title of the program: DPEMC, version 2.3

Computer: any computer with the FORTRAN 77 compiler under the UNIX or Linux operating systems.

Operating system: UNIX; Linux versions 7.x, 8.x

Programming language used: FORTRAN 77

High speed storage required: $`<`$ 25 MB

Keywords: Proton-(anti)proton collisions, diffraction, double pomeron exchange, double photon exchange.

Nature of the physical problem: Proton diffraction at hadron colliders can manifest itself in many forms, and a variety of models exist that attempt to describe it . This program implements some of the more significant ones, enabling the simulation of central particle production through color singlet exchange between interacting protons or antiprotons.

Method of solution: The Monte-Carlo method is used to simulate all elementary $`22`$ and $`21`$ processes available in HERWIG. The color singlet exchages implemented in DPEMC are implemented as functions reweighting the photon flux already present in HERWIG.

Restriction on the complexity of the problem: The program relying extensively on HERWIG, the limitations are the same as in .

Typical running time: Approximate times on a 800 MHz Pentium III : 5-20 minutes per 10000 unweighted events, depending on the process under consideration.

## 1 Introduction

Doubly diffractive Higgs boson production has been much debated during the last decade, in a variety of theoretical approaches and predictions . Now that the Tevatron is running, with both the D0 and CDF experiments well equipped for forward proton detection, there is good hope to clarify the situation in the mid-term. Meanwhile, potential upgrades of the LHC experiments ATLAS and CMS are being discussed, and a systematic evaluation of the DPE physics potential becomes necessary. The present work is meant as a step in this direction.

Available programs for diffraction simulation are POMPYT  and RAPGAP ; the factorized model of Cox and Forshaw, POMWIG ; and the Soft Color Interaction model of Enberg, Ingelman and Timneanu . POMWIG is implemented as an extension to HERWIG  and relies on the following branchings:

$$pp^{}+\mathrm{I}\mathrm{P},\mathrm{followed}\mathrm{by}\mathrm{I}\mathrm{P}q+X$$
(1)
where $`p`$ and $`p^{}`$ are the incoming and outgoing protons, the Pomeron is denoted $`\mathrm{I}\mathrm{P}`$, and $`q`$ is the parton (quark or gluon) involved in the hard subprocess. In this $`inclusive`$ picture, the final state can be made in lowest order of e.g. jets, photon pairs, lepton pairs, or a Higgs boson, accompanied by remnants from the $`\mathrm{I}\mathrm{PI}\mathrm{P}`$ collision and the outgoing protons.

Technically, the simulation of $`\gamma \gamma `$ interactions in $`ee`$ collisions, provided by HERWIG, is transformed into $`pp`$ induced $`\mathrm{I}\mathrm{PI}\mathrm{P}`$ interactions through a reweighting of the photon flux and a proper choice of the structure functions. The original implementation of this procedure, POMWIG, was based on a Regge fit to the Pomeron flux and Pomeron parton distribution functions provided by the H1 experiment .

Note that nothing prevents us from applying the same machinery to simulate $`exclusive`$ $`\mathrm{I}\mathrm{PI}\mathrm{P}`$ interactions, i.e. processes where the entire $`\mathrm{I}\mathrm{PI}\mathrm{P}`$ center-of-mass energy is used in the production of a heavy central system. Most often, the Pomeron is seen as a gluonic object, so that the branchings

$$pp^{}+\mathrm{I}\mathrm{P},\mathrm{followed}\mathrm{by}\mathrm{I}\mathrm{P}g+g_{NP}$$
(2)
result in a hard gluon-gluon interaction, with e.g. jet pairs, photon pairs or Higgs bosons in the final state. The color screening soft gluon, $`g_{NP}`$, couples to its counterpart from the opposite beam. Pomeron remnants are absent, and the protons again leave intact. Contrarily to the inclusive processes, which are being measured at the Tevatron (through dijet production ), the exclusive DPE has not yet been observed. This category of events is most appealing however, because Higgs boson production through this channel provides a very competitive measurement of its mass , and possible determination of its spin and parity, which is a major challenge in non-diffractive channels.

It is also simple to describe $`\gamma \gamma `$ interactions in the same framework : this is just going back to the original HERWIG where, to accommodate proton- or heavy ion-induced photon fluxes, we need to take into account the corresponding form factors. Compared to $`ee`$ induced $`\gamma \gamma `$ interactions, an additional complication comes from the fact that one needs to simultaneously forbid strong interactions between the nuclei. To achieve this, one generally works in impact parameter space and requires this impact parameter to be large enough for the nuclei not to geometrically overlap. These interactions have interest in their own right (they are a basic ingredient of the study of ultra-peripheral heavy-ion collisions, for example), but also provide a natural lower bound on the exclusive central production in $`pp`$ collisions. In the worst case, these channels thus provide a conservative estimate of the physics potential of DPE.

This flexible framework allows to implement a number of models alternative to that of Cox and Forshaw in an analogous way. We focus here on the exclusive model of Bialas and Landshoff (BL) , and on its inclusive extension by Boonekamp, Peschanski and Royon (BPR) . On the $`\gamma \gamma `$ side, we propose an implementation of the computations by Cahn and Jackson (CJ) , Drees, Ellis and Zeppenfeld (DEZ) , Papageorgiu , and Budnev et al .

We can unfortunately not pretend to be exhaustive, and the implementation of further models is left for future work. One notable example is the model of Khoze, Martin and Ryskin . We also escape discussing proton non-survival through soft parton rescattering; the cross-sections presented here should thus be understood as incomplete in that respect. Note however that the results of the inclusive DPE models can be rescaled to the measured DPE dijet cross-section, as in . On the other hand, in QED mediated diffraction, the proton collisions occur at much higher impact parameter, and the probability of parton rescattering is small .

The next section gives an overview of the models newly implemented (we leave most physical discussions aside and orient the reader to the references quoted above). Some results are then given in the third section, and the Appendix provides a guide to the package, and lists all settings and switches that control the output.

## 2 DPE models implemented in DPEMC

This section describes the new models and processes which can be studied using DPEMC. The model of Cox and Forshaw is described in  and of course still available.

### 2.1 Exclusive DPE in the Bialas-Landshoff model

The production cross section for a heavy system $`X`$ (such as a Higgs boson, a jet pair, or a photon pair) in exclusive DPE collisions:

$$p+pp+X+p,$$
(3)
as illustrated in Fig. 1a, takes the form :

$$d\sigma ^{(\mathrm{exc})}=𝑑\xi _1𝑑\xi _2F_{\mathrm{I}\mathrm{P}/p}(\xi _1)F_{\mathrm{I}\mathrm{P}/p}(\xi _2)𝑑\widehat{\sigma }(ggX).$$
(4)
The convolution in Eq. (4) is over the momentum fractions of the initial protons, $`\xi _1`$ and $`\xi _2`$, carried away by the Pomerons. The flux of the Pomeron inside the proton is parametrized as:

$$F_{\mathrm{I}\mathrm{P}/p}(\xi )=𝒩dv^2\frac{e^{\beta v^2}}{\xi ^{2\alpha _{\mathrm{I}\mathrm{P}}(v^2)1}},$$
(5)
where $`v^2`$ is related to the momentum transfer squared and the proton momentum loss through $`v^2=t/(1\xi )`$. Here the parameters $`𝒩^2=288/\pi ^5`$, $`\beta =4\mathrm{GeV}^2`$ and the Regge trajectory

$$\alpha _{\mathrm{I}\mathrm{P}}(t)=\alpha _{\mathrm{I}\mathrm{P}}(0)+\alpha _{\mathrm{I}\mathrm{P}}^{}t,$$
(6)
with intercept $`\alpha _{\mathrm{I}\mathrm{P}}(0)=1.08`$ and slope $`\alpha _{\mathrm{I}\mathrm{P}}^{}=0.25\mathrm{GeV}^2`$, are taken as in the original Bialas-Landshoff model, Ref. .

The partonic cross section $`d\widehat{\sigma }`$ in Eq. (4) describes the possible hard subprocesses $`ggX`$, of which a list is given in Table 1. In this context, the cross-sections satisfy color and angular momentum selection rules; namely the initial state is color singlet, and only the helicity configuration with $`J_z=0`$ and $`P=+1`$ contributes to $`\widehat{\sigma }`$ . These features have important consequences on the sensitivity to Higgs boson searches, since they imply important suppression of the backgrounds while being harmless to a naturally scalar and color-singlet signal.

Note that properly speaking, the Bialas-Landshoff expressions have been derived for Higgs bosons and massive quarks only. We propose here a generalization to gluon pair production, by using the corresponding color-singlet, $`J_z=0`$ subprocess cross-sections instead of the original ones.

### 2.2 Inclusive DPE

We implement the model of BPR as an alternative to that of Cox and Forshaw (CF) . Although different in spirit, both lead to similar expressions and final states. The production cross-section for a system $`X`$ ($`X`$ = Higgs, Dijet, Diphoton, Dilepton) in inclusive DPE collision (see Fig. 1b):

$$p+pp+X+p+\mathrm{Pomeron}\mathrm{remnants},$$
(7)
can be expressed as follows :

$$d\sigma ^{(\mathrm{inc})}=\underset{i,j}{}𝑑x_1𝑑x_2𝑑\xi _1𝑑\xi _2F_{\mathrm{I}\mathrm{P}/p}(\xi _1)F_{\mathrm{I}\mathrm{P}/p}(\xi _2)f_{i/\mathrm{I}\mathrm{P}}(x_1,\mu ^2)f_{j/\mathrm{I}\mathrm{P}}(x_2,\mu ^2)𝑑\widehat{\sigma }(ijX).$$
(8)
Besides the convolution over the proton momentum fractions, the integration runs over the Pomeron parton distribution functions. The Pomerons momentum fractions carried by the colliding partons are denoted $`x_1`$ and $`x_2`$. The proton induced Pomeron flux, $`F_{\mathrm{I}\mathrm{P}/p}(\xi _i)`$, is the same as in the exclusive case (the normalization and the parameters $`\beta `$, $`\alpha _{\mathrm{I}\mathrm{P}}(0)`$ and $`\alpha _{\mathrm{I}\mathrm{P}}^{}`$ are assumed identical to those used in ). The partonic content of the Pomeron is expressed in terms of distribution functions:

$$f_{i/\mathrm{I}\mathrm{P}}(x_i,\mu ^2)=x_iG_{i/\mathrm{I}\mathrm{P}}(x_i,\mu ^2),$$
(9)
where $`G_{i/\mathrm{I}\mathrm{P}}(x_i,\mu ^2)`$ are the “true densities” of the partons inside the Pomeron as measured by the H1 collaboration . They are evaluated at a scale $`\mu ^2`$ given by the mass of the central system. The integral of $`f_{i/\mathrm{I}\mathrm{P}}(x_i,\mu ^2)`$ is normalized to 1, so that in the limit $`f_{i/\mathrm{I}\mathrm{P}}(x_i,\mu ^2)\delta (x_i)`$ the exclusive expressions are recovered.

Note that the cross-sections obtained in this way are rescaled in  according to the double diffractive dijet cross-section measured in ; this is not done in the program.

Some of the possible hard subprocesses $`gg,q\overline{q}X`$ are listed in Table 2.

### 2.3 Two-photon reactions in $`pp`$ and heavy-ion collisions

The production cross section for a central system $`X`$ in heavy-ion collisions through double photon exchange (see Fig. 1c):

$$N+NN+X+N,$$
(10)
can be approximated by a factorized form:

$$d\sigma ^{(\gamma \gamma )}=𝑑\xi _1𝑑\xi _2F_{\gamma /\mathrm{N}}(\xi _1)F_{\gamma /\mathrm{N}}(\xi _2)𝑑\widehat{\sigma }(\gamma \gamma X).$$
(11)
The photon flux from an ion of electrical charge $`Z`$ and mass $`M_N`$ is parametrized as follows by DEZ :

$$F_{\gamma /N}(\xi )=\frac{\alpha Z^2}{\pi }\frac{1}{\xi }\left[e^{\xi ^2M_N^2/Q_0^2}+\left(1+\frac{\xi ^2M_N^2}{Q_0^2}\right)\mathrm{Ei}(\xi ^2M_N^2/Q_0^2)\right],$$
(12)
where the momentum transfer scale $`Q_060\mathrm{MeV}`$ and $`\mathrm{Ei}(x)=_x^{\mathrm{}}dte^t/t`$. Formula (12) is valid for heavy nuclei only. An alternative formulation, derived in impact parameter space, yields :

$$F_{\gamma /N}(\xi ,\xi _0)=\frac{2\alpha Z^2}{\pi \xi }\left[\frac{\xi }{\xi _0}K_0\left(\frac{\xi }{\xi _0}\right)K_1\left(\frac{\xi }{\xi _0}\right)\frac{1}{2}\left(\frac{\xi }{\xi _0}\right)^2\left(K_1^2\left(\frac{\xi }{\xi _0}\right)K_0^2\left(\frac{\xi }{\xi _0}\right)\right)\right].$$
(13)
In Eq. (13), $`K_{0,1}(\xi /\xi _0)`$ are the Bessel functions, while $`\xi _01/(RM_N)`$ with $`R`$ the radius of the ion. As discussed in the Introduction, this parameter appears because the impact parameter space integration is limited to the region where the colliding nuclei do not geometrically overlap. One can use $`R1.2A^{1/3}`$ for collisions of ions of atomic number $`A`$, or the proton radius to estimate the $`\gamma \gamma `$ cross-sections in proton-proton collisions. The partonic cross sections relevant in this context are listed in Table 3.

## 3 Results and Perspectives

All implemented fluxes are displayed in Figure 2. Apart from normalization, the CF and BL Pomeron fluxes differ mainly because of different choices of the Pomeron intercepts ($`\alpha _{\mathrm{I}\mathrm{P}}(0)=1.2`$ is taken by CF, while $`\alpha _{\mathrm{I}\mathrm{P}}(0)=1.08`$ is used by BL and BPR). The nucleus induced photon flux is illustrated for the proton ($`f_{\gamma /p}`$), Calcium ($`f_{\gamma /Ca}`$), and Lead ($`f_{\gamma /Pb}`$).

Figure 3 illustrates Higgs boson production cross-sections in various models, namely the inclusive models of CF and BPR, the exclusive cross-section by BL, and the photon-mediated cross-section by Papageorgiu. Figure 4 displays the “background” in the BPR model (the total and $`\mathrm{b}\overline{\mathrm{b}}`$ cross-sections are displayed) and in the exclusive BL model. The exclusive $`\mathrm{b}\overline{\mathrm{b}}`$ prediction is much suppressed as expected, and is more distributed towards small angles than its inclusive counterpart.

Considering the above examples and others that can be worked in an analogous way, it is possible using this program to perform complete studies of experimental sensitivity to diffractively produced Higgs bosons. These studies can be done following models of inclusive or exclusive double Pomeron exchange, and of double photon exchange in proton or heavy ion collisions. Supersymmetric enhancements of the quark and gluon couplings to the Higgs boson can also be included as in HERWIG.

The implementation of further models of hard diffraction, such as that of Ref. , will follow in subsequent versions of the program. Another important addition will be the inclusion of diffractive proton dissociation. Finally, since the beam particles currently need to be defined as electrons, it is impossible to simulate soft underlying events or multiple parton interactions. This shortcoming will also be handled in the near future.

In summary, the aim of this work is not to promote one model or the other, but to provide a tool that should allow to determine those that are compatible with the available and forthcoming data, and to prepare the analyses and constrain the predictions for the LHC.

## 4 Acknowledgments

We wish to thank B. Cox, J. Forshaw, V.A. Khoze, Ch. R. Peschanski, M. Rijssenbeek, C. Royon, M.H. Seymour and B.R. Webber for many helpful discussions. One of us (T.K.) is grateful to M. Rijssenbeek and CERN for support.

## Appendix A Manual

The source code for the program can be obtained from the authors or from Ref. , and includes an example main program. Some aspects of the implementation of DPEMC are described in Appendix B; here we discuss how to run the program. The processes are controlled using the following flags:

* IPROC,ID,IL select the hard process as in HERWIG;
* TYPEPR controls how the hard cross-section is to be evaluated; inclusively (TYPEPR=’INC’) or only the color-singlet $`J_Z=0`$ amplitude (TYPEPR=’EXC’);
* TYPINT switches between photon initiated processes (TYPINT=’QED’), and gluon/quark initiated processes (TYPINT=’QCD’);
* NSTRU controls the $`\mathrm{I}\mathrm{P}`$ or $`\gamma `$ flux; it can take values from 9 to 15, with meanings detailed in Table 4. For Pomeron fluxes, it controls the parton density functions to be used as well;
* when NSTRU=11, GAPSPR and CDFFAC set the absolute cross-section normalizations. GAPSPR controls the gap survival probability for exclusive DPE processes (TYPEPR=’EXC’); typical values are 0.1 at the Tevatron, and 0.03 at the LHC . CDFFAC controls the normalization of inclusive DPE processes (TYPEPR=’INC’), by comparison of the raw dijet computation to the measurement of ; the value 3.8 was found in .
* when appropriate, AION and ZION set the atomic and proton number of the colliding nuclei.

As in POMWIG, the beam particles need to be set as electrons in the main program, and the HERWIG facility generating $`\gamma \gamma `$ interactions is used and reweighted to the user’s needs, depending on the above settings.

Since NSTRU controls which flux is to be used in the cross-section evaluations, the introduction of TYPEPR and TYPINT may seem somewhat redundant. This is however needed for technical reasons. First, exclusive $`\mathrm{I}\mathrm{PI}\mathrm{P}`$ interactions can not be treated as such since the gluon initiating the hard process carries the full momentum of its parent “beam particle”, provoking errors in HERWIG. To circumvent this problem, we needed to modify a few HERWIG routines. For example, the routine evaluating $`\gamma \gamma f\overline{f}`$ has been modified so that it can compute the color-singlet, $`J_Z=0`$$`ggq\overline{q},gg`$ cross-section, depending on the values of TYPINT and TYPEPR.

For example, to simulate $`b\overline{b}`$ production in the exclusive Bialas-Landshoff model, one thus needs to request $`\gamma \gamma b\overline{b}`$ via IPROC=16005, set TYPINT=’QCD’, TYPEPR=’EXC’ so that the color-singlet, $`J_Z=0`$$`ggq\overline{q}`$ cross-section is actually computed, and NSTRU=11 requesting the Bialas-Landshoff flux. An exhaustive list of consistent flag settings is given in Table 5. These settings need to be followed strictly, the program is producing unpredictable results otherwise111We have tried to make the program detect unsound settings as much as possible..

It was also needed to provide some additional routines, for the computation of the $`\gamma \gamma H`$ and color-singlet, $`J_Z=0`$$`gg\gamma \gamma `$ cross-sections. For completeness, a computation of $`\gamma \gamma \gamma \gamma `$ is also provided. A list of the new IPROC values, and of ID settings specific to DPEMC is given in Table 7.

Finally, we provide an additional routine that corrects the event record after each event is finalized (remember that all events are recorded as $`ee`$ induced $`\gamma \gamma `$ interactions). This is needed for readability, but also to enable correct experimental simulations. When called, this routine resets the beam particles as protons, the secondary photon “beams” can be redefined as Pomerons or Reggeons, depending on the value of NSTRU, and the particles initiating the hard process, always stored as photons when TYPEPR=’EXC’, are redefined as gluons when TYPINT=’QCD’.

## Appendix B Structure of the program

We give here a list of the FORTRAN files that come with DPEMC. They have been tested to comply with HERWIG version 6.5. The first set of files contain the original HERWIG routines, which needed modifications for our purpose:

: the original source code of HERWIG; the routines HWSFUN, HWEGAM, HWEPRO and HWHQPM, superseded by the following files, are renamed to avoid conflicts;

: structure function evaluation;

: chooses the flux and generates Pomeron or photon momentum fractions and virtualities;

: decides on the hard process;

: calculates the cross-section for $`ggq\overline{q},gg`$ in the $`J_z=0`$ mode, for $`\gamma \gamma q\overline{q},W^+W^{},l^+l^{}`$ as in the original implementation of HERWIG, and for $`\gamma \gamma q\overline{q}`$ in the $`J_z=0`$ mode.

The following files contain the new routines needed to run DPEMC:

: implementation of the Pomeron and photon fluxes inside a proton or heavy ion;

: parton distribution functions inside the Pomeron;

: user defined Pomeron structure function;

: calculates the cross sections for $`ggH`$ and $`\gamma \gamma H`$;

: evaluates the cross sections for $`gg\gamma \gamma `$ in the case of exclusive diphoton production and for $`\gamma \gamma \gamma \gamma `$ in the heavy-ion collisions;

: calculates the normalization of parton distribution functions inside the Pomeron;

: corrects the event record.

In addition to the above mentioned files, a file vegas.f, containing the program VEGAS  for numerical integration, is also needed to obtain the normalization of Pomeron parton distribution functions. The example main programs provided with the package use HBOOK for histogramming, and thus need to be linked against the CERN program library .

An example main program is in example.f, and produces 10000 exclusive $`c\overline{c}`$ events in the Bialas-Landshoff model. It is compiled by executing gmake -f makeExample.

## Appendix C