text:
# Detection of dark matter concentrations in the field of Cl 1604+4304 from weak lensing analysis

## 1. INTRODUCTION

Weak shear fields of high-redshift galaxies are promising, efficient tools to investigate the mass distribution on cluster-supercluster scales (Kaiser & Squires 1993; Luppino & Kaiser 1997; Kaiser et al. 1998; Bartelmann & Schneider 2000) and provide a unique mean to detect dark mass concentrations (Schneider 1996; Erben et al. 2000). Wide-field weak-lensing surveys of projected mass overdensities can probe the statistical clustering properties and underlying cosmology (Bahcall et al. 1999).

Recent observations have revealed the existence of a supercluster at a high redshift of $`z0.9`$ (Lubin et al. 2000). This supercluster contains two massive galaxy clusters, Cl 1604+4304 at $`z=0.897`$ and Cl 1604+4321 at $`z=0.924`$. The two clusters are separated by $`17\mathrm{}`$ on the sky, corresponding to a projected separation of $`9h_{50}^1`$ Mpc. Cl1604+4304 is located at ($`\alpha _{\mathrm{J2000}}`$, $`\delta _{\mathrm{J2000}}`$) = ( $`16\mathrm{h}\mathrm{\hspace{0.17em}04}\mathrm{m}\mathrm{\hspace{0.17em}19.5}\mathrm{s}`$, $`+43\mathrm{°}\mathrm{\hspace{0.17em}04}\mathrm{}\mathrm{\hspace{0.17em}33}\stackrel{}{\mathrm{.}}9`$) and one of the optically-selected high-redshift cluster candidates studied by Oke, Postman, & Lubin (1998). Postman, Lubin, & Oke(1998) presented a detailed photometric and spectroscopic survey of the cluster and obtained a velocity dispersion of $`1226`$ km s-1 and a dynamical mass estimate of $`6.2\times 10^{15}h_{50}^1M_{}`$.

In this Letter, we present a weak lensing analysis of the Cl 1604+4304 field on the basis of deep images taken with the Hubble Space Telescope (HST)/Wide Field Planetary Camera 2 (WFPC2). We shall mainly describe the weak-lensing signal-to-noise ratio (S/N) analysis of the observed field; the details of the mass reconstruction procedure will be present elsewhere. Throughout this Letter, we adopt $`\mathrm{\Omega }_0=1`$ ,$`\mathrm{\Omega }_\mathrm{\Lambda }=0`$, and $`H_0=50h_{50}`$ km s-1 Mpc-1; $`1\mathrm{}`$ on the sky corresponds to $`0.52h_{50}^1`$ Mpc at the cluster redshift.

## 2. OBSERVATIONS AND DATA REDUCTION

The cluster Cl 1604+4304 field was observed in 1994 and 1995 using the WFPC2 camera with the F814W filter on board the HST (P.I.: J. Westphal; proposal 5234). We retrieved the calibrated data for Cl 1604+4304 field from the HST archive. Both the 1994 and 1995 observations have the same total exposure time of $`32`$ ksec, consisting of 16 single orbits of $`2`$ ksec. The 1995 pointing covers the central region of Cl 1604+4304 (Lubin et al. 1998), while the 1994 pointing is about $`20\mathrm{}`$ offset ($`\mathrm{\Delta }X0\stackrel{}{\mathrm{.}}1,\mathrm{\Delta }Y18\stackrel{}{\mathrm{.}}2`$) from the 1995 pointing. For this reason, Lubin et al. (1998) discarded the 1994 data from their analysis. We analyze both the 1994 and 1995 data separately. For each observation, the data were shifted and combined into the final frame to remove cosmic rays using the IRAF/STSDAS task CRREJ. The PC chip was discarded from our analysis because of its brighter isophotal limit, so that the final frame for each observation consists of three WFC chips. The side length of the WFPC2 field is about $`2\stackrel{}{\mathrm{.}}5`$ ( $`1.26h_{50}^1`$ Mpc at $`z=0.897`$).

We use the SExtractor package (Bertin & Arnouts 1996) for the object detection, photometry, and measurement of image shapes. We extract all objects with isophotal areas larger than 12 pixels ($`0\stackrel{}{\mathrm{.}}0996`$ pixel-1) above $`2\sigma `$ pixel-1 of the local sky level. We calculate the total magnitudes $`I_{814\mathrm{W}}`$ in STMAG system for each object using MAG\_BEST parameters in SExtractor. The corresponding detection thresholds for the 1994 and 1995 data are $`\mu _{814\mathrm{W}}=26.2`$ and $`26.5`$ mag arcsec-2, respectively. We measure the quadruple moments $`Q_{ij}=d^2\theta \theta _i\theta _jI(\stackrel{}{\theta })/d^2\theta I`$ of the surface brightness $`I(\stackrel{}{\theta })`$ for each object, where the center of light is chosen as the coordinate origin. From $`\{Q_{ij}\}`$, we construct the image ellipticity $`ϵ_i=(Q_{11}Q_{22},2Q_{12})/(\mathrm{Tr}(Q)+2(\mathrm{det}Q)^{1/2})`$ $`(i=1,2)`$. To construct catalogs of faint galaxies, we refer to the following SExtractor output parameters: We exclude objects with FLAGS (extraction flags) $`1`$ from our analysis. We make star/galaxy separation using the CLASS\_STAR parameter (stellarity index) and keep the objects with $`\mathrm{CLASS}\mathrm{\_}\mathrm{STAR}0.2`$ as faint galaxy candidates. Objects with FWHM\_IMAGE (FWHM profile from a Gaussian fit to the core) $`<4`$ pixels are excluded since image shapes of extremely small objects relative to the pixel scale may be affected by the anisotropic PSF. Finally, all objects with $`I_{814\mathrm{W}}(25.5,27.5)`$ are selected as background galaxy candidates. This detection and selection procedure leads to the final catalogs with total galaxy numbers of $`N_\mathrm{g}=`$ 200 and 241, corresponding number densities of $`n_\mathrm{g}=`$ $`45.8`$ and $`56.8`$ arcmin-2, for the 1994 and 1995 data, respectively.

## 3. WEAK LENSING ANALYSIS

### 3.1. Strategy

Our first aim is to obtain the distribution of weak-lensing S/N in the data field. Then, the high peaks in the S/N maps can be identified as clusters or mass overdensities. To do this, we make use of a variant of aperture mass statistics. The statistics rely on the fact that the shear $`\gamma =(\gamma _1,\gamma _2)`$ and the convergence $`\kappa :=\mathrm{\Sigma }/\mathrm{\Sigma }_{\mathrm{cr}}`$ are related to each other through $`\stackrel{}{}\kappa =\widehat{D}\gamma \stackrel{}{u}_\gamma `$, where $`\mathrm{\Sigma }`$ is the surface mass density of the deflector, $`\mathrm{\Sigma }_{\mathrm{cr}}=(c^2/4\pi G)D_\mathrm{s}/D_\mathrm{d}D_{\mathrm{ds}}`$ is the critical surface mass density, and $`\widehat{D}=\{\widehat{D}_{ij}\}`$ ($`i,j=1,2`$) is a differential operator defined by

$$\widehat{D}=\left(\begin{array}{cc}/\theta _1& /\theta _2\\ /\theta _2& /\theta _1\end{array}\right)$$
(1)
(Kaiser 1995). Since $`\stackrel{}{u}_\gamma `$ is the gradient of $`\kappa `$, operating $`\widehat{D}`$ further on $`\stackrel{}{u}_\gamma `$ yields

$$(\mathrm{}\kappa ,0)=\widehat{D}^2\gamma =\widehat{D}\stackrel{}{u}_\gamma =(\mathrm{div}\stackrel{}{u}_\gamma ,\mathrm{rot}\stackrel{}{u}_\gamma ).$$
(2)
In weak lensing limit ($`\kappa 1`$ and $`|\gamma |1`$), the expectation value of the image ellipticity, $`\mathrm{E}[ϵ(\stackrel{}{\theta })]`$, is the shear $`\gamma (\stackrel{}{\theta })`$. In practice, however, background sources have intrinsic ellipticities $`ϵ_{(\mathrm{s})}`$, so that weak lensing analysis involves the smoothing procedure to reduce the noise. We denote the smoothed fields by angular brackets $``$: e.g., $`\kappa =d^2\theta ^{}W(|\stackrel{}{\theta }\stackrel{}{\theta }^{}|;\vartheta )\kappa (\stackrel{}{\theta }^{})`$, where $`W(\theta ;\vartheta )`$ is a smooth, continuous window function with a characteristic scale of $`\vartheta `$. Because of the commutativity between smoothing and the mass reconstruction (Van Waerbeke 2000), the smoothed quantities $`\kappa `$ and $`\gamma `$ satisfies the same relations as those between $`\kappa `$ and $`\gamma `$: $`\stackrel{}{}\kappa =\widehat{D}\gamma \stackrel{}{u}_\gamma `$ ; $`(\mathrm{}\kappa ,0)=\widehat{D}^2\gamma =\widehat{D}\stackrel{}{u}_\gamma =(\mathrm{div}\stackrel{}{u}_\gamma ,\mathrm{rot}\stackrel{}{u}_\gamma )`$. The Laplacian of $`\kappa `$, i.e. $`\mathrm{div}\stackrel{}{u}_\gamma `$, is just the convergence convolved with the compensated filter function $`\mathrm{}W`$, and hence equivalent to the aperture mass. Defining $`\stackrel{}{u}_ϵ:=\widehat{D}ϵ`$, we see that the observable $`\mathrm{div}\stackrel{}{u}_ϵ`$ traces the distribution of $`\mathrm{\Sigma }`$ in the limit of weak lensing, as pointed out by Luppino & Kaiser (1997). On the other hand, $`\mathrm{rot}\stackrel{}{u}_ϵ`$ measures the ‘pure’ noise. A discretized estimator for $`\widehat{D}\stackrel{}{u}_\gamma `$ is given by $`\widehat{D}\stackrel{}{u}_ϵ=(\mathrm{div}\stackrel{}{u}_ϵ,\mathrm{rot}\stackrel{}{u}_ϵ)=n_\mathrm{g}^1_{m=1}^{N_\mathrm{g}}p(|\stackrel{}{\theta }\stackrel{}{\theta }_m|;\vartheta )[ϵ_\mathrm{t}(\stackrel{}{\theta }_m;\stackrel{}{\theta }),ϵ_\mathrm{r}(\stackrel{}{\theta }_m;\stackrel{}{\theta })]`$, where $`ϵ_\mathrm{t}`$ and $`ϵ_\mathrm{r}`$ are the tangential and the radial components of the image ellipticity $`ϵ=ϵ_1+iϵ_2`$ defined by $`ϵ_\mathrm{t}(\stackrel{}{\theta };\stackrel{}{\theta }_0):=\mathrm{}[ϵ(\stackrel{}{\theta })e^{2i\varphi }]`$ and $`ϵ_\mathrm{r}(\stackrel{}{\theta };\stackrel{}{\theta }_0):=\mathrm{}[ϵ(\stackrel{}{\theta })e^{2i\varphi }]`$ with $`\varphi =\mathrm{Arg}(\stackrel{}{\theta }\stackrel{}{\theta }_0)`$ and $`p(\theta ;\vartheta )=W^{\prime \prime }(\theta ;\vartheta )W^{}(\theta ;\vartheta )/\theta `$. The noise properties of $`\widehat{D}^2ϵ`$ due to the intrinsic source ellipticities are contained in the covariant matrix $`\sigma _{ij}^2:=\mathrm{E}[\widehat{D}^2ϵ^{(\mathrm{s})}\widehat{D}^2ϵ^{(\mathrm{s})}]_{ij}`$:

$$\sigma _{ij}^2(\stackrel{}{\theta })=\delta _{ij}\frac{\sigma _ϵ^2}{2n_\mathrm{g}^2}\underset{m=1}{\overset{N_\mathrm{g}}{}}p^2(|\stackrel{}{\theta }\stackrel{}{\theta _m}|;\vartheta )\delta _{ij}\sigma ^2(\stackrel{}{\theta }),$$
(3)
where $`\sigma _ϵ`$ is the dispersion of the intrinsic source ellipticities and we have assumed $`\mathrm{E}[ϵ^{(\mathrm{s})}(\stackrel{}{\theta }_m)ϵ^{(\mathrm{s})}(\stackrel{}{\theta }_n)]_{ij}=\delta _{ij}\delta _{mn}\sigma _ϵ^2/2`$. The Kronecker delta $`\delta _{ij}`$ in equation (3) ensures that the dispersions of $`\mathrm{div}\stackrel{}{u}_ϵ`$ and $`\mathrm{rot}\stackrel{}{u}_ϵ`$ are the same and that the two fields are statistically uncorrelated. The local weak-lensing S/N $`\nu (\stackrel{}{\theta })`$ at position $`\stackrel{}{\theta }`$ is then defined by

$$\nu (\stackrel{}{\theta }):=\frac{\mathrm{div}\stackrel{}{u}_ϵ(\stackrel{}{\theta })}{\sigma (\stackrel{}{\theta })}=\frac{\sqrt{2}}{\sigma _ϵ}\frac{_mp(|\stackrel{}{\theta }\stackrel{}{\theta }_m|;\vartheta )ϵ_\mathrm{t}(\stackrel{}{\theta }_m;\stackrel{}{\theta })}{\sqrt{_mp^2(|\stackrel{}{\theta }\stackrel{}{\theta _m}|;\vartheta )}}.$$
(4)
The resulting formula (eq. ) is equivalent to the one derived by Schneider (1996) using aperture mass statistics. In the present Letter, we use a Gaussian window function of the form $`W_\mathrm{G}(\theta ;\vartheta )=\mathrm{exp}(\theta ^2/\vartheta ^2)/\pi \vartheta ^2`$, in which case $`p(\theta ;\vartheta )=4(\theta /\vartheta )^2\mathrm{exp}(\theta ^2/\vartheta ^2)/\pi \vartheta ^4p_\mathrm{G}(\theta ;\vartheta )`$; $`p_\mathrm{G}(\theta ;\vartheta )`$ has its maximum at $`\theta =\vartheta `$ and falls off rapidly at $`\theta >\vartheta `$. These statistics can be applied to the strong-lensing regime (e.g., cluster central region), because the contribution of image ellipticities to the aperture mass comes mainly from galaxies within an annulus at radius $`\vartheta `$ and thus we can avoid the strong lensing regime. For each independent HST observation, we perform a local S/N analysis using equation (4).

To obtain the projected mass distribution, we perform a mass reconstruction to the HST/WFC field. Taking into account the high redshift ($`z=0.897`$) of Cl 1604+4304 and small field-of-view ($`2\stackrel{}{\mathrm{.}}5`$ on a side) of the HST/WFC field, we adopt a non-linear finite-field inversion method developed by Seitz & Schneider (1997), which takes account of the source redshift distribution. Since little is known about the redshift distribution of field galaxies, we assume a source redshift distribution of the form $`p_z(z)=\beta z^2\mathrm{exp}[(z/z_0)^\beta ]/\mathrm{\Gamma }(3/\beta )z_0^3`$ (Brainerd, Blandford, & Smaili 1996), in which case the mean redshift $`z`$ is given by $`z=z_0\mathrm{\Gamma }(4/\beta )/\mathrm{\Gamma }(3/\beta )`$. In the present Letter, we consider only the case $`(z,\beta )=(1.0,1.0)`$.

Once a smoothed ellipticity field $`ϵ(\stackrel{}{\theta })`$ is obtained from the observed image ellipticities, a convergence map can be obtained through the integral equation

$$\kappa _{\mathrm{}}(\stackrel{}{\theta })\overline{\kappa }_{\mathrm{}}=_𝒰d^2\theta ^{}\stackrel{}{H}(\stackrel{}{\theta },\stackrel{}{\theta }^{})\stackrel{}{u}_\gamma _{\mathrm{}}(\stackrel{}{\theta }^{}),$$
(5)
where quantities with $`\mathrm{}`$-subscript represent the values for sources at infinite redshift, $`\overline{\kappa }_{\mathrm{}}`$ is the unknown constant which represents the average of $`\kappa _{\mathrm{}}`$ within the data field $`𝒰`$, $`\stackrel{}{H}`$ is the kernel which is the gradient of the scalar field that satisfies the Neumann boundary problem (see Seitz & Schneider 1996), and $`\stackrel{}{u}_\gamma _{\mathrm{}}=\widehat{D}\gamma _{\mathrm{}}`$. In general, the shear is not direct observable, so that the integral equation (5) is non-linear and solved iteratively: A mass reconstruction scheme applied to the irregular HST/WFC field is outlined in Seitz et al. (1996) and Umetsu, Tada, & Futamase (1999).

In order to determine the constant $`\overline{\kappa }_{\mathrm{}}`$, we employ a dynamical mass estimate of CL 1604+4304 obtained by Postman et al. (1998). Postman et al. (1998) estimated the projected mass inside the circular aperture of radius $`R_{500}=500h_{50}^1`$ kpc at the dynamical center to be $`M_{500}=0.842\times 10^{15}h_{50}^1M_{}`$, corresponding to the mean convergence within the circular aperture of $`\overline{\kappa }_{\mathrm{},500}=M_{500}/\pi R_{500}^2\mathrm{\Sigma }_{\mathrm{cr},\mathrm{}}=0.811`$. However, owing to the small field-of-view of the HST/WFPC2, the observed data fields do not cover the whole aperture area even for the 1995 pointing. The effective aperture area covered by the 1995 data field is about $`60\%`$ of the circular aperture area. We assume that the mean convergence $`\overline{\kappa }_{\mathrm{},\mathrm{eff}}`$ inside the effective aperture is the same as the mean convergence $`\overline{\kappa }_{\mathrm{},500}`$ inside the circular aperture: $`\overline{\kappa }_{\mathrm{},\mathrm{eff}}\overline{\kappa }_{\mathrm{},500}=0.811`$. We, thus, solve the integral equation (5) iteratively varying the value of $`\overline{\kappa }_{\mathrm{}}`$ in such a way that the average $`\overline{\kappa }_{\mathrm{},\mathrm{eff}}`$ of the reconstructed convergence within the effective aperture area equals $`\overline{\kappa }_{\mathrm{},500}=0.811`$.

### 3.2. Analysis of Cl 1604+4304

Since the 1995 pointing covers the central cluster region where a dynamical analysis was performed by Postman et al. (1998), we first analyzed the 1995 data. We calculated the fields $`\widehat{D}\stackrel{}{u}_ϵ`$ with $`\vartheta =0\stackrel{}{\mathrm{.}}4`$ ($`150h_{50}^1`$ kpc) on $`24\times 24`$ grid points $`\stackrel{}{\theta }_{ij}`$ (except on the PC chip), from which we obtained the distribution of the weak-lensing S/N within the WFC field. Here we have measured the dispersion $`\sigma _ϵ`$ from the rotation field $`\mathrm{rot}\stackrel{}{u}_ϵ(\stackrel{}{\theta }_{ij})`$ in the following way: We calculated the variance $`\sigma _{\mathrm{rot}}^2`$ of the noise field $`\mathrm{rot}\stackrel{}{u}_ϵ(\stackrel{}{\theta }_{ij})`$ by averaging over the grid points. Equating $`\sigma _{\mathrm{rot}}^2`$ with the average of $`\sigma ^2(\stackrel{}{\theta })`$ over the grid points $`\stackrel{}{\theta }_{ij}`$, we obtained an estimate of $`\sigma _ϵ=0.21`$.

Figure 1 displays the resulting weak-lensing S/N map obtained for the 1995 data. We find two S/N peaks above a threshold of $`\nu _{\mathrm{th}}=3`$. The first maximum has a peak height of $`\nu =4.5`$ and its location is about $`1\stackrel{}{\mathrm{.}}7`$ ($`850h_{50}^1`$ kpc) southwest of the the optical cluster center (marked with cross) determined from the dynamical analysis of Postman et al. (1998). The peak hight of the second maximum is $`\nu =3.6`$ and located in the vicinity of the dynamical cluster center. This offset is about $`0\stackrel{}{\mathrm{.}}15`$ ($`75h_{50}^1`$ kpc) and comparable to the grid spacing of $`0\stackrel{}{\mathrm{.}}1`$ ($`50h_{50}^1`$ kpc), so that the second peak is consistent with the dynamical cluster center and associated with the cluster member galaxies (Lubin et al. 1998). In this field, there are a total of 12 cluster members whose redshifts are spectroscopically confirmed by Keck observations (Lubin et al. 1998). However, no cluster members are observed at around the first S/N peak (see Fig. 4 of Lubin et al. 1998) because objects located on the upper edge of the WF4 chip (upper-right WF chip) are not contained within the sampled region of spectroscopic observations (Postman et al. 1998).

Although the detection of the mass concentration cannot be confirmed spectroscopically from the available data, we have checked the stability and reliability of S/N maps. We find that the main features, i.e. the locations and heights of the peaks, are insensitive to the choice of the weight function or filtering scale. It shoud be noted that these main features are stable even if objects with large ellipticities ($`|ϵ|0.4`$) are removed from our galaxy catalog, in which case the peak heights of the first and second maxima are $`3.7`$ and $`3.0`$, respectively. There still remains a possibility that the detected S/N peaks might be noise peaks due to the intrinsic source ellipticities or the anisotropic PSF. Van Waerbeke (2000) has shown from the simulated data that the probability distribution of the noise peaks in mass maps deviates from Gaussian and the probability of finding high noise peaks increases owing to the boundary effects. Since the field-of-view of the HST/WFC is small and we adopt a constant filtering scale, the obtained S/N maps may be affected by boundary effects. We can investigate such possibilities using independent deep observation of this field. In Fig. 2, we show the weak-lensing S/N map obtained from the 1994 data. Here we have used the same filtering scale and grid spacing as for the 1995 data. Comparing Fig. 2 with Fig. 1, we see that the first maximum in the 1994 S/N map corresponds to that in the 1995 S/N map. The first maximum in the 1994 S/N map has a peak height of $`\nu =4.2`$, whose location coincides well with that of the first maximum in the 1995 S/N map ($`\alpha _{\mathrm{J2000}}=16\mathrm{h}\mathrm{\hspace{0.17em}04}\mathrm{m}\mathrm{\hspace{0.17em}15.7}\mathrm{s}`$ and $`\delta _{\mathrm{J2000}}=+43\mathrm{°}\mathrm{\hspace{0.17em}03}\mathrm{}\mathrm{\hspace{0.17em}47}\stackrel{}{\mathrm{.}}2`$). The identification of the first S/N maximum from two-different pointing observations indicates that the S/N maximum does not originate from the anisotropic PSF. Moreover, since the first maximum in the 1994 S/N map is located at the center of the WF4 chip, the S/N measure at the first peak cannot be affected by boundary effects. These results confirm the detection of ‘dark’ mass concentration associated with the first S/N maximum.

We have reconstructed the projected mass distribution within the 1995 data field. We calculated the smoothed ellipticity $`ϵ(\stackrel{}{\theta })`$ on $`20\times 20`$ grid points (except on the PC chip) using $`W_\mathrm{G}(\theta ;\vartheta )`$ with $`\vartheta =0\stackrel{}{\mathrm{.}}35`$ ($`176h_{50}^1`$ kpc). In Fig. 3, we show the reconstructed $`\kappa _{\mathrm{}}`$-map. In mass reconstruction, we implicitly assumed that the dark mass concentration has the same redshift $`z=0.897`$ as Cl 1604+4304. We see from Figs. 1 and 3 that the local maxima in the S/N maps correspond to the mass peaks and that the weak-lensing S/N traces the projected mass, as expected. The $`\kappa _{\mathrm{}}`$-map shows a super-critical feature around the dark mass concentration ($`\kappa _{\mathrm{}}=1.19`$) as well as in the cluster central region. The projected mass within the 1995 data field is obtained as $`1.0\times 10^{15}h_{50}^1M_{}`$. From the mass map, the enclosed mass in a $`280h_{50}^1`$ kpc radius centered on the dark mass concentration is calculated to be $`2.4\times 10^{14}h_{50}^1M_{}`$: For the mass-to-light ratio, a lower limit is obtained as $`M/L_\mathrm{B}500h_{50}`$ in solar units from an upperlimit on the dark clump’s surface brightness.

## 4. DISCUSSION AND CONCLUSIONS

We have performed a weak-lensing analysis on the deep images of a region around Cl 1604+4304 taken with the HST/WFPC2. We detected two significant maxima in the resulting weak-lensing S/N map: the second peak associated with the dynamical cluster center and the first peak located about $`1\stackrel{}{\mathrm{.}}7`$ south of the dynamical center having no optical counter parts. The identification of high S/N peak from two-different pointing observations provides a robust evidence of dark mass concentration. Assuming that the dark mass clump has the same redshift as Cl 1604+4304, we reconstructed the projected mass distribution in the 1995 data field, from which we have estimated the mass within a circular aperture of radius $`280h_{50}^1`$ kpc centered on the dark clump to be $`2.4\times 10^{14}h_{50}^1M_{}`$. The obtained results indicate that this dark clump could be a compact, massive bound system associated with the supercluster at $`z0.9`$; the obtained mass will be overestimate if the redshift of the dark clump is much less than $`0.897`$.

Our results demonstrate that the weak-lensing S/N statistics are powerful and efficient tools even for high-redshift cluster surveys ($`z1`$). For a further confirmation of the dark mass concentration, deep wide-field observations covering the entire supercluster field ($`20\mathrm{}`$) are required. Such observations can probe the mass distribution in the supercluster field extending over $`10h_{50}^1`$ Mpc.

We are very grateful to Toru Yamada, Masahiro Takada, and Masaru Kajisawa for helpful discussions. We thank the anonymous referee for useful comments. Part of this work is based on the observation with the NASA/ESA Hubble space Telescope, obtained from the data archive at the Space Telescope Science Institute (STScI). STScI is operated by the Association of Universities for Research in Astronomy, Inc. under the NASA contract NAS 5-26555.


format: text

text:
# Observational constraints upon quintessence models arising from moduli fields

## I Introduction

A cosmological constant $`\mathrm{\Lambda }`$ has been considered as the missing energy of the universe for a long time. Recently, a varying vacuum energy or “quintessence”  has become popular as an alternative candidate. Since more parameters are involved in such models, they are more flexible in explaining some of the problems left by the cosmological constant model. One such problems is the “cosmic coincidence problem” : the missing energy and the matter energy densities decrease at different rates as the universe expands, so it seems purely coincident that they are comparable today with the concordant values from several observational tests of $`\mathrm{\Omega }_{\varphi 0}\frac{2}{3}`$ and $`\mathrm{\Omega }_{m0}\frac{1}{3}`$ . Several candidates for quintessence fields have been proposed. Typically quintessence models possess attractor solutions with common evolutionary properties for a wide range of initial conditions. For example, in some cases the scalar field energy maintains constant ratio to the matter density at late time , where as in other models it tracks the dominant matter component in some other sense .

In this paper, I will study quintessences arising from string moduli. In string or Kaluza-Klein type models the moduli fields associated with the geometry of the extra dimensions may have effective potentials which depend exponentially on the moduli fields, due to the curvature of the internal spaces, or alternatively through the interaction of the moduli with the form fields on the internal spaces (see, e.g.,  and references therein). Single exponential potentials of the form

$$V=V_0e^{\lambda \kappa \varphi }$$
(1)
which give rise to a scaling solution have been well-studied in the literatures (see, e.g.,  and references therein). In this paper I will concentrate on double expontial potentials of the form

$$V=V_0\mathrm{exp}(Ae^{\sqrt{2}\kappa \varphi })$$
(2)
which typically arise from supersymmetry breaking via gaugino condensation . Furthermore, it has also recently been argued that if the superpotential of the moduli field in string theory is T-duality invariant, then it cannot be approximated by a single exponential function, but must depend on double exponentials .

Binetruy , and later de la Macorra , have argued that scalar fields with double exponential potentials cannot act as quintessences for two reasons. Firstly, such models have a global attractor solution leading to a matter-dominated universe at late times. Secondly, near the attractor the equation of state is positive, $`w_\varphi >0`$ and approaching zero as $`t\mathrm{}`$. This would appear to contradict the latest observational results from type IA supernovae that a quintessence is dominating over matter and the universe is entering a phase of accelerated expansion , similar to inflation. The purpose of this paper is to show that for parameter values away from the attractor, there can exist models which are consistent with the observational tests, although some tuning is required for the universe not to have reached the attractor at present. This paper also updates the observational constraints on the single exponential potential model, which have been given by Frieman and Waga .

The organization of this paper is as followed: in section II, I will present a phase-space analysis of the double exponential potential model. In section III, I will discuss the numerical integration of the evolution equations for both the single and double exponential potential models, and obtain values for $`\mathrm{\Omega }_{\varphi 0}`$ and $`H_0t_0`$. In section IV, I will constrain both models using the light-curve calibration luminosity distances of type IA supernovae (see  and references therein) and the gravitational lensing statistics of high luminosity quasars by intervening galaxies (see  and references therein). Similar constraints on other quintessence models have been presented in .

## II Phase space analysis for double exponential potential model

We begin by considering a universe which consists of a scalar field with the potential (2) and a barotropic fluid with equation of state $`P_\gamma =(\gamma 1)\rho _\gamma `$, $`0\gamma 2`$. For a spatially-flat Friedmann-Robertson-Walker (FRW) universe, the governing equations are given by

$`\dot{H}`$ $`=`$ $`{\displaystyle \frac{\kappa ^2}{2}}\left(\rho _\gamma +P_\gamma +\dot{\varphi }^2\right),`$ (3)
$`\dot{\rho _\gamma }`$ $`=`$ $`3H(\rho _\gamma +P_\gamma ),`$ (4)
$`\ddot{\varphi }`$ $`=`$ $`3H\dot{\varphi }{\displaystyle \frac{dV}{d\varphi }},`$ (5)
subject to the Friedmann constraint

$$H^2=\frac{\kappa ^2}{3}\left(\rho _\gamma +\frac{1}{2}\dot{\varphi }^2+V\right),$$
(6)
where an overdot denotes the ordinary differentiation with respect to the time $`t`$. In the above equations, $`\kappa ^28\pi G`$, and $`H\dot{a}/a`$ is the Hubble parameter.

We follow the treatment of Copeland, Liddle and Wands  by defining the variables

$$x=\frac{\kappa \dot{\varphi }}{\sqrt{6}H},y=\frac{\kappa \sqrt{V}}{\sqrt{3}H},$$
(7)
and introducing a third variable

$$\lambda (\varphi )=\frac{dV/d\varphi }{\kappa V}=\sqrt{2}Ae^{\sqrt{2}\kappa \varphi }.$$
(8)
In the case of the single exponential potential (1), $`\lambda (\varphi )`$ reduces to a constant. In terms of these variables the evolution equations (3) – (6) become

$`x^{}`$ $`=`$ $`3x+\lambda \sqrt{{\displaystyle \frac{3}{2}}}y^2+{\displaystyle \frac{3}{2}}x\left[2x^2+\gamma (1x^2y^2)\right],`$ (9)
$`y^{}`$ $`=`$ $`\lambda \sqrt{{\displaystyle \frac{3}{2}}}xy+{\displaystyle \frac{3}{2}}y\left[2x^2+\gamma (1x^2y^2)\right],`$ (10)
$`\lambda ^{}`$ $`=`$ $`2\sqrt{3}x\lambda ,`$ (11)
where a prime denotes a derivative with respect to the logarithm of the scale factor, $`N\mathrm{ln}(a)`$. The system is bounded by the cylinder $`x^2+y^2=1`$. Assuming that $`A`$ is positive, the system is confined to $`\lambda 0`$. The system is symmetric under the reflection $`yy`$.

As many properties of the solutions in the double exponential potential model can be related to the single exponential potential model, let us review the single exponential potential model as has been studied in the literature ( and references therein). Depending on the values $`\lambda `$ and $`\gamma `$, there are up to five critical points in the $`x`$$`y`$ plane, two of which are related to late times. For $`\lambda ^2<3\gamma `$, the late times attractor is a scalar field dominated solution $`\mathrm{\Omega }_\varphi \kappa ^2\rho _\varphi /3H^21`$ with $`\gamma _\varphi (\rho _\varphi +P_\varphi )/\rho _\varphi \frac{\lambda ^2}{3}<\gamma `$. For $`\lambda ^2>3\gamma `$, the late times attractor is the scaling solution  with $`\gamma _\varphi =\gamma `$ and $`\mathrm{\Omega }_\varphi \frac{3\gamma }{\lambda ^2}`$, a constant. Fig. 1 shows the late times attractors and how they evolve as $`\lambda ^2`$ increases from zero. As $`\lambda ^2`$ increases the attractor evolves from a stable node into a stable spiral at $`\lambda ^2=\frac{24\gamma ^2}{(9\gamma 2)}`$. The attractor approaches the origin as $`\lambda ^2\mathrm{}`$.

In the double exponential potential model, we can identity critical points at finite $`\lambda `$ from the evolution equations (9) – (11). They are three discrete points at $`(x_c,y_c,\lambda _c)=`$ $`(\pm 1,0,0)`$, $`(0,1,0)`$, and a 1-parameter family at $`(0,0,\lambda _c)`$ for $`\lambda _c`$ arbitrary. The stability of the critical points can be obtained by linearizing (9) – (11) about the these points and solving for the eigenvalues of small perturbations. The critical points and the eigenvalues are listed in Table I.

The scalar field kinetic energy dominated solution $`(1,0,0)`$ is an unstable node, while $`(1,0,0)`$ is a saddle point since trajectories are repelled in the $`x`$ and $`y`$ directions and are attracted in the $`\lambda `$ direction. The scalar field dominated solution $`(0,1,0)`$ attracts a two-dimensional bunch of trajectories but is degenerate in the $`\lambda `$ direction. For the barotropic fluid dominated solution $`(0,0,\lambda _c),0\lambda _c<\mathrm{}`$, trajectories are attracted in the $`x`$ direction, repelled in the $`y`$ direction and degenerate in the $`\lambda `$ direction.

To examine the critical points at infinite $`\lambda `$, we use the Poincaré sphere method by transforming

$$\lambda =\frac{1}{ϵ}.$$
(12)
On the surface $`ϵ=0`$ we find

$`{\displaystyle \frac{dx}{d\tau }}=\sqrt{{\displaystyle \frac{3}{2}}}y^2`$ (13)
$`{\displaystyle \frac{dy}{d\tau }}=\sqrt{{\displaystyle \frac{3}{2}}}xy`$ (14)
where $`\tau `$ is a new time coordinate defined by $`d\tau =\lambda dN`$. Therefore, there are infinitely many critical points lying on the $`ϵ=0`$ surface with $`y_c=0`$ and arbitrary $`x_c[1,1]`$. The eigenvalues for the matrix of perturbations are $`0`$, $`\sqrt{\frac{3}{2}}x_c`$, $`0`$. Therefore the critical points are degenerate in both the $`x`$ and $`\lambda `$ direction. They repel or attract trajectories in the $`y`$ direction depending on $`1x_c<0`$ or $`0<x_c1`$ respectively. For $`x_c=0`$ all the eigenvalues vanish, and higher-order perturbations show that it is a stable attractor.

We will numerically integrate the evolution equations (9) – (11) forwards and backwards from some initial point $`(x_i,y_i,\lambda _i)`$, in the case of a matter-dominated universe $`\gamma =1`$. In Fig. 3 and 3, we project the trajectories onto the $`x`$$`y`$ plane. When integrating backwords from the initial points $`(x_i,y_i,\lambda _i)`$, the trajectories in Fig. 3 approach the critical point $`(1,0,0)`$, and the trajectories in Fig. 3 approach the critical points $`(x_c,0,\lambda _c)`$, $`x_c[1,0)`$, $`\lambda \mathrm{}`$.

The scalar field is rolling down the potential at sufficiently late times so $`\varphi `$ and hence $`\lambda (\varphi )`$ is increasing with time. As $`\lambda (\varphi )`$ increases with time, the behaviour of the solutions effectively approximates that of single exponential potential models, but with the parameter $`\lambda `$ evolving in the direction shown in Fig. 1. For sufficiently small $`\lambda _i`$ the scalar field-dominated solution is dynamically important and the trajectories get very close to the circular boundary. As $`\lambda `$ increases, the solutions approach scaling solutions ($`x=y`$ for $`\gamma =1`$) before spiralling towards the origin.

## III Numerical integration

In order to determine $`\mathrm{\Omega }_{\varphi 0}`$ and $`H_0t_0`$, and the luminosity distance – redshift relation in the next section, we ought to integrate the evolution equations (3) – (6) numerically. In this section, we will discuss the numerical integration for both the single and double exponential potential models in a manner similar to a previous study on the model with a pseudo Nambu-Goldstone boson (PNGB) potential .

### A Single exponential potential

Here we update the numerical integration result for the single exponential potential by Frieman and Waga . We use an alternate set of dimensionless variables:

$$u=\frac{\kappa \dot{\varphi }}{H_0},v=\mathrm{\Omega }_{m0}^{1/3}(1+z),w=\lambda \kappa \varphi \mathrm{ln}\frac{\kappa ^2V_0}{H_0^2},$$
(15)
the variables $`u`$ and $`v`$ are the same as were used in . In terms of these variables, the field equations become:

$`u^{}`$ $`=`$ $`3{\displaystyle \frac{H}{H_0}}u+\lambda e^w,`$ (16)
$`v^{}`$ $`=`$ $`{\displaystyle \frac{H}{H_0}}v,`$ (17)
$`w^{}`$ $`=`$ $`\lambda u,`$ (18)
where a prime denotes a derivative with respect to the dimensionless time parameter $`H_0t`$. The Hubble parameter is defined implicitly according to

$$\frac{H}{H_0}\left(v^3+\frac{1}{6}u^2+\frac{1}{3}e^w\right)^{\frac{1}{2}}.$$
(19)
We shall begin the integration at last-scattering. We choose $`H_0t_i=0`$. At this initial stage, $`H/H_01`$ and the field evolution (16) is overdamped by the expansion of the universe, driving $`u`$ to zero. We therefore choose $`u_i=0`$. At last-scattering, $`z1100`$ and $`\mathrm{\Omega }_{m0}^{1/3}1`$ for $`\mathrm{\Omega }_{m0}1`$, hence $`v_i=1100`$ is a reasonable value. Note that results of the integration do not change significantly if $`H_0t_i`$, $`u_i`$ and $`v_i`$ are altered to values within the same order of magnitude. The remaining parameters are $`w_i`$ and $`\lambda `$, which are determined by the physical origin of the model.

The integration proceeds until a value $`H/H_0=1`$ is reached. In Figs. 4 and 5, the contours of $`\mathrm{\Omega }_{\varphi 0}=1v_0^3`$ and $`H_0t_0`$ are displayed in the $`w_i`$$`\lambda `$ parameter space. In view of recent estimates of the ages of globular clusters , a mean value of 12.8 Gyr for the age of the universe is arrive at. With $`h0.7`$ this would require $`H_0t_00.9`$. A concordant value of $`\mathrm{\Omega }_{\varphi 0}=0.67\pm 0.05`$ is given by several observational tests . By comparing these two contour plots, we see that parameters in the region of small $`\lambda `$ are favoured. This region shall represent models where the solutions approach the attractor, either the scalar field dominated solution ($`\lambda ^2<3\gamma `$) or the scaling solution ($`\lambda ^2>3\gamma `$), only recently. Note that a seperate bound of $`\lambda ^2>20`$ has been obtained from nucleosynthesis by assuming that the late times attractor is the scaling solution, and the attractor is approched within a few expansion times of the end of inflation .

Another useful quantity which can be determined from the numerical integration is the deceleration parameter

$$q_0=\frac{3}{2}v_0^3+\frac{1}{2}u_0^21.$$
(20)
$`q_0<0`$ corresponds to a universe whose expansion is accelerating at the present epoch.

### B Double exponential potential

With the double exponential potential model, we again use the variables $`u`$ and $`v`$ as defined in (15), and redefine $`w`$ as

$$w=\sqrt{2}\kappa \varphi +\mathrm{ln}A.$$
(21)
In terms of these variables the field equations become

$`u^{}`$ $`=`$ $`3{\displaystyle \frac{H}{H_0}}u+\sqrt{2}V_1e^w\mathrm{exp}\left(e^w\right),`$ (22)
$`v^{}`$ $`=`$ $`{\displaystyle \frac{H}{H_0}}v,`$ (23)
$`w^{}`$ $`=`$ $`\sqrt{2}u,`$ (24)
where we have replaced $`V_0`$ with a dimensionless parameter $`V_1=\frac{\kappa ^2V_0}{H_0^2}`$, and the Hubble parameter is

$$\frac{H}{H_0}\left[v^3+\frac{1}{6}u^2+\frac{1}{3}V_1\mathrm{exp}\left(e^w\right)\right]^{\frac{1}{2}}.$$
(25)
We choose the same $`H_0t_i`$, $`u_i`$, and $`v_i`$ as above. The parameters $`w_i`$ and $`V_1`$ determine the phyisical origin of the model. The contours of $`\mathrm{\Omega }_{\varphi 0}`$ and $`H_0t_0`$ are displayed in the $`w_i`$$`V_1`$ parameter space in Fig. 7 and 6. By comparing these two contour plots, the interesting region would seem to be the bottom left-hand corner of the parameter space. The interesting region corresponds to the scalar fields being still nearly frozen to their initial states, or having become dynamical and starting to roll down the potential slope only recently.

## IV Observational Constraints

A constraint on the single exponential potential model has been obtained by Frieman and Waga  using five type IA supernovae (Sne IA). In this section, we will constrain both the single and double exponential potential models using a larger set of Sne IA, together with a separate constaint using the gravitational lensing statistics of high luminosity quasars.

Empirical calibration of the light curve – luminosity relationship of Sne IA provides absolute magnitudes that can be used as distance indicators. The Supernova Cosmology Project leaded by Perlmutter et al.  and the High Redshift Supernovae Search Team leaded by Riess et al.  are two different groups have been working extensively in searching for high redshift Sne IA. I will use the larger available data set of the two, namely the 60 Sne IA published by Perlmutter et al. . The luminosity distance – redshift relation so obtained provides a good observational constraint on quintessence models provided that there is no intrinsic evolution of the peak luminosities of the Sn IA sources. The issue of how the constraints are altered in the presence of such evolution has been discussed in .

Gravitational lensing of distant quasars due to galaxies along the line of sight provides another relatively sensitive constraint  on the quintessence models. The statistics of abundances of multiply imaged quasars and observed separations of the images to the source can be used to estimate the distances to the quasars. We follow the calculation as described by Waga and Miceli , who used a total of 862 ($`z>1`$) high luminosity quasars plus 5 lenses from seven major optical survey.

We display the results of both tests as 68.3% and 95.4% joint credible regions of the parameter spaces described in the previous section. Figs. 8 and 9 display results for the parameter spaces of the single and double exponential potential models respectively. For both models, the interesting regions that give rise to expected values of $`\mathrm{\Omega }_{\varphi 0}0.67`$ and $`H_0t_00.9`$ are well within the 68.3% confident region of both the Sne IA and the gravitational lensing statistics tests. The regions to the left of the $`q_0=0`$ contours correspond to models that give rise to accelerated expansion.

## V Discussion

We have performed a phase-space analysis on the double exponential potential model, its properties can be understood in terms of a single exponential potential model with a varying coefficient $`\lambda `$. The analysis shows that unlike the single exponential potential model which has two late times attractor depending on $`\lambda `$, the double exponential potential model has only one global attractor which leads to a matter dominated universe at late times. However, it is always possible for the model to be dominated by the scalar field at present and dominated by matter in the future. The problem with this model is we need to fine-tune the parameters in order not to have reached the attractor at the present epoch, thereby circumventing the objections of Binetruy  and de. la. Macorra . However, all other quintessence models appear to also require a degree of fine-tuning .

We studied the observation constraints for both the single and double exponential potential models using a more updated type IA supernovae data and the gravitational lensing statistics. The results show that there are regions in the parameter spaces for which the models are consistent with the observations at the same time giving appropriate values for the scalar field energy density and the age of the universe at present.

## Acknowledgments

I would like to thank Chris Kochanek for supplying me with the gravitational lensing data, Nelson Nunes, Ioav Waga, and David Wiltshire for helpful discussions about aspects of the paper.


text:
# Energy Spectra, Altitude Profiles and Charge Ratios of Atmospheric Muons

## I Introduction

Measurements of the flux of atmospheric neutrinos by large underground detectors have consistently disagreed with theoretical predictions, a fact that has been interpreted in terms of possible neutrino oscillations. The most compelling evidence thus far comes from the SuperKamiokande experiment . The discrepancies between experiment and theory are well beyond the statistical uncertainties in the measurements. The correct interpretation of the effect requires a detailed understanding of neutrino production in the atmosphere.

Model predictions  of the absolute flux of neutrinos produced in air showers are uncertain due to a number of difficulties: one problem is the normalization of the primary cosmic-ray flux, for which measurements vary by $`\pm `$15% or more. Because of this, different assumptions are made by different authors. For instance, the model of Honda et al. , used as the starting point in the analyses of the SuperKamiokande collaboration, assumes a primary cosmic ray flux normalization that is in excess of the current measurements by balloon experiments , by as much as $``$20% in the energy range beyond about 5 GeV, which is relevant to the production of sub- to multi-GeV atmospheric neutrinos. In contrast, the model of the Bartol group  uses a primary spectrum in good agreement with the most recent measurements. Different assumptions are also made about the details of particle production in atmospheric interactions, and about the treatment of geomagnetic effects (see Ref. ). Thus, while the authors of these calculations expect that the predictions of absolute neutrino fluxes might have a total uncertainty of about 18% , it seems fortuitous that the predicted neutrino rates have been found to agree with each other within about 10%. It should also be noted that in the neutrino oscillation analysis of the SuperKamiokande collaboration, an additional normalization factor of 1.16 in the cosmic ray intensity is introduced to obtain the best oscillation fit to the measured muon and electron neutrino rates as a function of zenith angle.

Neutrino production in the atmosphere is closely coupled to muon production, as both types of particles are produced together in pion and kaon decays, and as some of the muons decay to produce further neutrinos. Monte Carlo simulations of neutrino production naturally predict the spectra of muons as a function of atmospheric depth. Therefore, a detailed experimental test of the predictions is possible through measurements of the intensity of muons at different depths in the atmosphere, for instance with a balloon-borne particle detector during atmospheric ascent. If such a detector includes a magnet spectrometer, separate measurements of the $`\mu ^{}`$ and $`\mu ^+`$ fluxes are possible. This is relatively easy for negative muons, as negatively charged particles other than electrons are rare in air showers (and electrons are easily rejected), while non-interacting protons generate a large background for positive muons above $``$1 GeV. Such measurements have been made by a number of instruments, including the MASS , HEAT , CAPRICE , and IMAX  magnet spectrometers and older, less sensitive balloon payloads.

Here we describe measurements made with the High Energy Antimatter Telescope (HEAT), a balloon instrument optimized for the study of high-energy cosmic-ray electrons and positrons. HEAT uses several complementary particle identification techniques, which are also well suited to the identification of muons during atmospheric ascent of the balloon. This instrument was flown twice, in 1994 and 1995. Preliminary results from the first flight on relative abundances of $`\mu ^+`$ and $`\mu ^{}`$ were reported previously . The present paper describes measurements made during the second balloon flight: the relative abundances of $`\mu ^+`$ and $`\mu ^{}`$ at momenta between 0.3 and 0.9 GeV/c, and the differential spectra of $`\mu ^{}`$ at momenta between 0.3 and 50 GeV/c, at atmospheric depths between 13 and 960 g/cm2. We compare the results with other measurements and with calculations made with a modified version of the TARGET Monte Carlo algorithm , developed by the Bartol group to predict atmospheric neutrino rates.

## II Muon Identification and Backgrounds

The HEAT instrument is described in detail elsewhere . It combines a superconducting magnet spectrometer (using a drift-tube tracking hodoscope), a time-of-flight (TOF) system, a transition-radiation detector (TRD) and an electromagnetic calorimeter (EMC). For its second flight, it was launched from Lynn Lake, Manitoba, Canada, at a vertical geomagnetic cutoff rigidity of well below 1 GV. The flight took place on 23 August 1995, approximately at minimum solar activity. The ascent from an atmospheric overburden at the ground of 960 g/cm2 to a float altitude at 3 g/cm2 required 3.1 hours, during which charged atmospheric secondary particles were detected and recorded. Near the end of ascent, at about 13 g/cm2, the instrument trigger configuration was changed to commence preferential measurements of electrons and positrons, affecting the ability to measure absolute muon intensities at shallower depths. Therefore we report measurements of the muon charge ratio at all atmospheric depths, but measurements of absolute muon intensities only at depths greater than 13 g/cm2.

In identifying muon events, TRD information is not used. The TOF system measures the velocity of the particle $`\beta =v/c`$ with a resolution of $`\sigma _\beta =0.15`$, permitting complete rejection of upward-going particles; in addition, the amount of light generated in the scintillation counters of the TOF system is used to measure the magnitude of the particle’s electric charge $`Ze`$ with a resolution of $`\sigma _Z=0.11`$, permitting the unambiguous identification of singly-charged particles. The magnet spectrometer measures the sign of the particle’s charge from the direction of the deflection in a magnetic field of about 1 T, with a field integral over the particle’s trajectory of $`B𝑑l4.2`$ kGm. The magnetic rigidity $`R=pc/Ze`$ of the particle is determined from the amount of deflection; at rigidities between 0.3 and 0.9 GV, the resolution achieved is $`\sigma _R=0.08`$ to 0.11 GV. The EMC records the pattern of energy deposits in 10 scintillation counters, each topped by a 0.9 radiation length-thick lead sheet. In each layer, the energy deposited is measured in units of the energy loss by a vertical minimum-ionizing particle (MIP). A shower sum is obtained by adding the signals from the 10 scintillators.

The selection of muon events proceeds in three steps. First, a high-quality spectrometer track and measurement of the rigidity $`R`$ are required, as described in Ref. . Also required are down-going and singly-charged TOF characteristics. Second, as shown in Figure 1, the distribution of reconstructed particle rigidity $`R`$ is studied as a function of velocity $`\beta `$, and compared with ideal curves $`R=(m/Z)\beta /\sqrt{(1\beta ^2)}`$ where $`m`$ is the mass and $`Z`$ the charge of the particle. Clearly identifiable in the figure are populations due to He (and d), p, $`\mu ^+`$ and $`\mu ^{}`$ events. $`\pi ^\pm `$ and some $`K^\pm `$ events cannot be distinguished from the $`\mu ^\pm `$ populations, and are a small background to the muon signal. We select $`\mu `$-like events by requiring: 1) $`\beta 0.85`$ and 2) $`0.3\left|R\right|0.9`$ GV, for the low-energy $`\mu ^+/\mu ^{}`$ ratio, and 2) $`R0.3`$ GV for the $`\mu ^{}`$ energy spectra. For $`R0.9`$ GV, $`\mu ^+`$ events become indistinguishable from non-showering hadrons (mostly protons), so that $`\mu ^+`$ spectra are not measured with this instrument. Third, we study the behavior of the shower sum measured by the EMC as a function of particle rigidity $`R`$, as shown in Figure 2. Ideal curves for electrons and positrons are obtained by assuming a simple linear relationship between EMC sum (proportional to energy) and $`R`$. For heavier particles, a calculation of EMC sum is made by integrating Bethe-Bloch energy losses within the EMC scintillators. Populations of events due to e±, p and $`\mu ^\pm /\pi ^\pm /K^\pm `$ are identifiable. By augmenting the selection criteria of Figure 1 with the requirement that particles not shower in the EMC (EMC sum $``$ 15), e± events with $`\left|R\right|>0.3`$ GV are rejected. The low-rigidity proton events that range out in the EMC and would appear to contaminate the $`\mu ^+`$ population are rejected by the $`\beta 0.85`$ requirement of Figure 1.

With these criteria, we achieve essentially complete rejection of electron events, but there remains a small background to the muon signal due to pions and kaons. We estimate from Monte Carlo simulations of air showers based on CERN’s GEANT-FLUKA algorithms  that the $`\pi ^\pm `$ flux at a depth of 13 g/cm2 is only 2% that of $`\mu ^\pm `$, in agreement with another calculation  – with $`K^\pm `$ fluxes at a much lower level – and that this further decreases with increasing atmospheric depth. Moreover, only pions that do not interact can be mistaken for muons, which occurs 39% of the time, so that the background to the muon measurement due to atmospheric pions is only about 0.7-0.9% near float, decreasing to less than 0.4% at depths greater than 300 g/cm2. Such a small background is not corrected for here. Occasionally, cosmic-ray interactions in the instrument result in $`\pi ^\pm `$ production, at an even more modest level than the atmospheric pion background. GEANT-based simulations indicate that only 0.04% of proton-induced events yield a misidentification as a $`\mu ^{}`$. As the payload slowly rotated throughout the flight (as determined with a solar sensor attached outside the gondola), and did not align itself with the Earth’s magnetic field, possible geomagnetic East-West asymmetries are averaged out. Furthermore, such asymmetries in the primary proton flux are only expected at momenta near the geomagnetic cutoff, which for the Lynn Lake flight is well below the energies of interest here.

## III Results and Discussion

### A Muon Charge Ratio

The number of low-energy $`\mu ^+`$ and $`\mu ^{}`$ events detected as a function of atmospheric depth and the $`\mu ^+/\mu ^{}`$ ratio are summarized in Table I. The $`\mu ^+/\mu ^{}`$ ratio is shown as a function of atmospheric depth in Figure 3 (labeled HEAT 95), together with the measured ratio from the first HEAT flight  (labeled HEAT 94) and other recent measurements . As noted on the figure, the various analyses have used different magnetic rigidity ranges, and moreover, the measurements were made at different solar epochs and different geomagnetic rigidity cutoffs Rcutoff, so that the data are not truly directly comparable. The two HEAT measurements of the muon ratio are the most statistically significant, and are essentially consistent with each other and with other measurements, within the appreciable errors. Based on the HEAT measurements, no clear correlation of the muon ratio with geomagnetic cutoff rigidity is observed. The depth dependence of the charge ratio seems to be essentially flat, although a slight decrease from about 1.3-1.4 at high altitudes (3-50 g/cm2) down to about 1.1 at the ground cannot be excluded. We note that the charge ratio measured by both HEAT 95 and CAPRICE 94 at small depths appears anomalously high. This effect is not understood at present. Also shown on Figure 3 are calculations with the TARGET algorithm  (widely used for neutrino flux calculations), for conditions of solar minimum and maximum activity. Both HEAT flights occurred under essentially solar-minimum conditions. The TARGET calculations have been made for average primary fluxes at solar minimum and maximum, which may not exactly represent the actual spectrum at the time of the flight. The calculations are for a location with no geomagnetic cutoff, and are intended for comparison with the HEAT 95 data only. The agreement between the HEAT measurements and the solar-minimum TARGET calculation appears to be fairly good. Note that the fluctuations in the simulated distributions at shallow atmospheric depths are statistical, and indicative of the small number of muons having been produced in the air showers at the highest altitudes.

### B Energy Spectra of Negative Muons

The absolute intensity of $`\mu ^{}`$, in a rigidity interval $`\mathrm{\Delta }R`$ with an average rigidity $`\overline{R}`$, at atmospheric depth $`d`$ is obtained with:

$$j_\mu (d,\overline{R})=\frac{N_\mu }{\mathrm{\Delta }tϵ_tϵ_lϵ_\theta ϵ_{dt}\mathrm{\Delta }R(\mathrm{\Omega }A)ϵ_{scan}ϵ_{acc}},$$

where $`N_\mu `$ is the number of $`\mu ^{}`$ events recorded at $`(d,\overline{R})`$, $`\mathrm{\Delta }t`$ is the time spent at depth $`d`$, $`ϵ_t`$ is the live time fraction, $`ϵ_l`$ is an event-transmission loss correction, $`ϵ_\theta `$ is a correction to account for the fact that the muon flux is increasingly less isotropic deeper in the atmosphere, $`ϵ_{dt}`$ is the efficiency of the basic event cleanliness criteria applied to the drift-tube hodoscope track, $`(\mathrm{\Omega }A)`$ is the geometrical factor, $`ϵ_{scan}`$ is a “scanning efficiency” correction (described below), and $`ϵ_{acc}`$ is the muon acceptance efficiency. The acceptance of the instrument decreases rapidly for particles incident at a zenith angle greater than about 25, so that particle intensities reported here are essentially for vertical incidence.

Both $`\overline{R}`$ and $`\mathrm{\Delta }R`$ are weighted to account for the details of the energy spectrum, according to:

$$\overline{R}=\frac{_{R_i}^{R_j}Rf(R)𝑑R}{_{R_i}^{R_j}f(R)𝑑R}\mathrm{\Delta }R=\frac{_{R_i}^{R_j}f(R)𝑑R}{f(\overline{R})}$$

where $`f(R)R^\alpha `$ is the rigidity power-law spectrum, with spectral index $`\alpha `$ varying between $`0.56`$ and 3.5 depending on both $`R`$ and $`d`$ ($`\alpha `$ is experimentally determined from the spectra before any of the normalization corrections are applied). $`\mathrm{\Delta }t`$ is measured with an on-board clock, $`ϵ_t`$ is determined using on-board scalers which count clock cycles while the instrument is available for a trigger or busy processing an event, and $`ϵ_l`$ is determined by careful accounting of event numbers generated on board compared to events successfully transmitted. $`ϵ_\theta `$ is calculated using a standard prescription , where the zenith dependence of the muon flux is taken to be $`\mathrm{cos}^{n(d)}\theta `$, with the exponent $`n`$ a function of atmospheric depth $`d`$; we have used $`n(d)=(d/1030\mathrm{g}/\mathrm{cm}^2)\times n_{sealevel}`$ with $`n_{sealevel}=2`$. $`ϵ_{dt}`$ is obtained by a careful accounting of the number of events recorded compared to the number of events with a successful minimal track reconstruction. $`(\mathrm{\Omega }A)`$ and $`ϵ_{acc}`$ are determined with the aid of a GEANT-based simulation of the response of the HEAT instrument. $`ϵ_{scan}`$ is a correction factor introduced based on the visual scanning of several hundred events to account for residual differences between the reconstruction efficiency of real events compared to that of simulated ones, and is found to be $`ϵ_{scan}=(0.9\pm 0.1)`$. The various parameters described above are given in Tables II and III.

The final $`\mu ^{}`$ intensities as a function of momentum for various atmospheric depths are shown in Figure 4 and given in Table IV. Also shown in the figure are the measurements of the MASS  and CAPRICE  experiments. The HEAT sample of 10327 $`\mu ^{}`$ events collected during ascent is to be compared with the MASS samples of 2893 events (1989 flight) and 4471 events (1991 flight) and the CAPRICE sample of 4627 events. Although the 1989 MASS measurements were also made in Northern Canada (from Prince Albert, Saskatchewan), the flight occurred at a different solar epoch (1989), at the time of a significant Forbush decrease. The 1991 MASS measurements were made from Fort Sumner, New Mexico. The CAPRICE data were collected at Lynn Lake, in 1994, and so are more directly comparable with our measurements. The general level of agreement between the data sets should be noted.

### C Comparison with Model Calculations

#### 1 One-Dimensional TARGET Algorithm

In Figure 5, we compare the HEAT measurements reported here with predictions of the TARGET algorithm , for conditions of solar minimum and maximum activity, shown as solid and dotted curves, respectively. (The solar-minimum curves are the ones of interest here, but the solar maximum curves are shown as well to illustrate the extent of the effect of the solar cycle on muon production.) These curves are obtained with the standard TARGET algorithm, which simulates vertically incident cosmic rays, and which follows the development of the air shower in one dimension only. No corrections for geomagnetic effects are made. This is the algorithm developed by the Bartol group and used in predicting underground neutrino rates from atmospheric sources . In Figure 6, we show the $`\mu ^{}`$ growth curves for different momentum intervals, also compared with the 1D TARGET-based predictions for solar minimum conditions (solid curves). The calculations were not made for the highest momentum bin.

There is general similarity between the experimental and simulated distributions, with some notable differences. For instance, the predictions are significantly in excess of the measurements below 4 GeV/c at atmospheric depths between 13 and 250 g/cm2. The ratio of simulated to measured intensity varies from $`1.2\pm 0.2`$ near shower maximum at 200 g/cm2 to $`1.7\pm 0.3`$ at depths between 13 and 140 g/cm2. At depths beyond shower maximum, or at momenta greater than 4 GeV/c, the simulations agree very well with our measurements. A similar trend was found by the CAPRICE collaboration : simulations predict more $`\mu ^{}`$ events than they measure below about 1 GeV/c, but they find that the ratio of simulated to measured intensity is greatest at shower maximum, with a value of $`1.8\pm 0.1`$.

#### 2 Three-Dimensional TARGET Algorithm

In an attempt to understand the origin of the discrepancy between the predicted and measured muon intensities, a new version of the TARGET algorithm was produced in collaboration with T. Gaisser and T. Stanev of the Bartol Research Institute. In this, three dimensional air shower development effects are taken into account, and the primary cosmic ray arrival direction is sampled isotropically, rather than assuming vertical incidence. Geomagnetic effects are not yet included in the calculations. Figure 5 also shows the $`\mu ^{}`$ momentum spectra at various atmospheric depths obtained with the 3-dimensional TARGET algorithm (dashed curves), for solar minimum conditions. The 3D calculations are in substantially better agreement with the data than the 1D calculations. In Figure 6, the measured $`\mu ^{}`$ growth curves are also compared with 3D TARGET predictions (dashed curves). Here again the 3D predictions are a more adequate representation of the data.

Figure 7 shows distributions of the ratios of predicted to measured $`\mu ^{}`$ intensities, for 1D TARGET calculations (top panel) and 3D TARGET calculations (bottom panel), respectively. These are cumulative distributions for all momentum and atmospheric depth bins. Each ratio is weighted by the square of the error on the ratio derived from the experimental error on the intensity. The distribution for 1D calculations has a mean 1.13, indicating an average overprediction of 13%, whereas for 3D calculations the mean is 1.07, a slightly better agreement. The main improvement however is in the reduced RMS variance of 0.17 for the distribution for 3D calculations compared to 0.27 for the 1D calculations. Thus, the more realistic calculations that take into account 3D air shower development and primary zenith arrival direction constitute a clear improvement in the representation of muon production.

#### 3 Neutrino Production

The 1D and 3D TARGET algorithms were used to predict ($`\nu _\mu +\overline{\nu }_\mu `$) and ($`\nu _\mathrm{e}+\overline{\nu }_\mathrm{e}`$) intensities at different atmospheric depths and in different momentum bins. The calculations are made for solar minimum conditions, for no geomagnetic rigidity cutoff, and for primary cosmic rays arriving within 30 of the zenith. The resulting neutrino growth curves in different momentum intervals are shown in Figure 8. The calculations are made only up to 32 GeV/c. Although there are differences between the 1D and 3D predictions at momenta less than about 1 GeV/c, these differences are most important at mid-to-high altitudes. The neutrino intensities at the ground level, which are the ones of relevance to the underground neutrino studies, are summarized in Table V, for 1D and 3D calculations. The neutrino intensities at ground level appear not to be altered much by the 3D effects. A similar conclusion was also reached in a study by Battistoni et al. , where detailed calculations of atmospheric muon and neutrino production are made in one and three dimensions.

## IV Conclusions

We have made statistically significant measurements of air shower muons as a function of atmospheric depth. We report the muon charge ratio $`\mu ^+/\mu ^{}`$ in the momentum range 0.3-0.9 GeV/c and the momentum spectra of $`\mu ^{}`$ in the range 0.3-50 GeV/c, at atmospheric depths from 13 to 960 g/cm2. The charge ratio is essentially constant with altitude within errors, with a possible decrease from 1.3-1.4 at high altitudes to 1.1 at the ground. A comparison of our measured $`\mu ^{}`$ momentum distributions with model calculations indicates significant discrepancies with the predictions of the standard one-dimensional TARGET algorithm: our measured fluxes are lower than the calculated ones at shallow depths before about shower maximum. Calculations of the muon intensities with a new version of the TARGET algorithm, accounting for three-dimensional air shower development, lead to a substantially improved agreement with our data. A detailed representation of atmospheric secondary production thus benefits from the more realistic simulations. The average excess of about 7% of the 3D calculations over our measured intensities is comparable to possible systematic effects in our experiment. Thus, within this uncertainty, the 3D TARGET algorithm generates atmospheric secondary particle intensities which are in agreement with the measurements.

The three-dimensional air shower development effects do not appear to impact significantly the atmospheric neutrino rates at the ground, but merely the pattern of neutrino production altitudes. Thus, we estimate that the neutrino intensities predicted by the 1D version of the TARGET algorithm are also accurate to about 7%. This is to be compared with the accuracy of 14-18% first estimated for such calculations . Even though the Honda et al.  model uses different assumptions about the primary cosmic ray flux and about the atmospheric interaction characteristics, it predicts neutrino intensities on the ground which agree with the TARGET predictions and with our data at the level of 7-10%. Therefore, one might put into question the additional normalization factor of 1.16 of the primary cosmic ray spectrum that is introduced in the SuperKamiokande neutrino oscillation analysis.

###### Acknowledgements.

We are grateful to T. Stanev and T. K. Gaisser for helpful discussions, for sharing with us their TARGET Monte Carlo algorithm and for contributing to the effort to introduce 3D effects in the algorithm. We thank the NSBF balloon crews that have supported the HEAT flights. This work was supported by NASA grants NAG5-5059, NAG5-5069, NAG5-5070, NAGW-5058, NAGW-1995, NAGW-2000 and NAGW-4737, and by financial assistance from our universities.
.

.


text:
# Far Infrared Study of IRAS 00494+5617 & IRAS 05327-0457

## 1 Introduction

A long term programme of studying the distribution of cold dust ( down to $``$ 15 K) in and around Galactic star forming regions, is in progress using the 1-m TIFR balloon-borne far-infrared (FIR) telescope. Under this programme, high angular resolution ($``$ 1′) mapping is carried out simultaneously in two trans-IRAS FIR wavebands (Ghosh et al., 1996; Mookerjea et al., 1998, 1999; Verma et al., 1999). In this paper we present FIR mapping of the sources IRAS 00494+5617 and 05327-0457 in wavebands centered around 143 & 185 $`\mu `$m. These sources share the similarity that each is heated by a luminous external source in addition to one or more possible embedded sources. The simultaneity of observations in the two wavebands with nearly identical beams is useful in deriving reliable temperature distributions of the interstellar dust in these regions.

The source IRAS 00494+5617 is part of the western fragment of the molecular cloud NGC281, located in the Perseus arm at a distance of 2.2 kpc (Cesaroni, Felli & Walmsley, 1999). A compact star cluster containing the multiple star HD 5005 as the brightest star excites an ionization front at the edge of NGC 281W. HD 5005 is an O6 V type star and the cluster is at an angular distance of $``$ 5′ from IRAS 00494+5617. Observational evidences in the form of detection of C34S (3$``$2) (Megeath & Wilson, 1997), NH3 (Henning et al., 1994), 22 GHz H2O masers (Tofani et al., 1995), CO emission (Carpenter, Snell & Schloerb, 1990) and molecular outflows (Snell, Dickman & Huang, 1990) together with FIR and millimeter continuum emission, strongly suggest ongoing star formation activity in this region. Existing literature on this source suggests that the ongoing star formation activity has been induced by the compression due to the propagation of the ionization front (energized by HD 5005) into the molecular cloud (Elmegreen & Lada, 1978; Megeath & Wilson, 1997). The gas emission features from this source are reasonably well studied. There are not as many high angular resolution observations of the emission from the dust component. Longward of IRAS wavelengths, the only available observations are the high resolution maps at 1.3 mm (Henning et al., 1994).

The source IRAS 05327-0457 is associated with the nebulosity NGC 1977, located near the northern star in the sword of Orion. It is bounded on the south by the northern end of the Orion molecular cloud (Kutner, Evans & Tucker, 1976) and is at a distance of 450 pc. The exciting star for the HII region is the B1 V star HD 37018, also known as 42 Ori and is at an angular distance of $``$ 5$`\stackrel{}{\mathrm{.}}`$9 from IRAS 05327-0457. Earlier observations of the NGC 1977 region include FIR ($``$ 160 $`\mu `$m), near-infrared and radio continuum mapping together with extensive molecular line observations (Makinen et al., 1985; Kutner et al., 1985). High resolution molecular line maps (Kutner et al., 1985) indicate an increase in radial velocity from south to north implying an expansion of the HII region into the molecular cloud, with a velocity of a few km/sec.

In the present study the structural details and dust temperature distributions obtained from FIR observations have been used to quantify the relative contributions of the external and internal sources towards heating these regions. In sections 2 & 3, the observations and the results are described. The discussion including the radiation transfer modeling of these sources is presented in section 4.

## 2 Observations

### 2.1 The 143 & 185 $`\mu `$m maps

The Galactic star forming regions associated with the sources IRAS 00494+5617 and 05327-0457 were observed using a two band far-infrared photometer system at the Cassegrain focus of the TIFR 1-m (f/8) balloon-borne telescope. The FIR telescope was flown from the TIFR Balloon Facility, Hyderabad, in central India (Latitude = 17$`\stackrel{}{\mathrm{.}}`$47 N , Longitude = 78$`\stackrel{}{\mathrm{.}}`$57 E ) on 1995 November 12. Details of the telescope and the observational procedure have been given by Ghosh et al. (1988). Additional information specific to this balloon flight has been presented by Mookerjea et al. (1999). The photometer consists of 12 composite silicon bolometers, each having a FOV of 1$`\stackrel{}{\mathrm{.}}`$6 and arranged in a 3$`\times `$2 array for each band. The sky was chopped at 10 Hz by wobbling the secondary mirror. The chopper throw was 4$`\stackrel{}{\mathrm{.}}`$2 along the cross elevation axis. The same region of the sky was viewed simultaneously in two bands. The effective wavelengths for the two bands are 143 $`\mu `$m and 185 $`\mu `$m for a greybody spectrum with a temperature of 36 K and $`\lambda ^1`$ emissivity. All flux densities presented in this paper also use the same assumptions regarding temperature and emissivity. Saturn was observed for absolute flux calibration as well as for the determination of the instrumental Point Spread Function (PSF) including the effect of sky chopping. The absolute flux calibration was done following the method described by Ghosh et al. (1988).

The regions (refer to the Figures Far Infrared Study of IRAS 00494+5617 & IRAS 05327-0457 and Far Infrared Study of IRAS 00494+5617 & IRAS 05327-0457) around IRAS 00494+5617 and 05327-0457 were mapped by scanning the sky in cross-elevation with steps in elevation at the end of each scan line. The chopped FIR signals were gridded in a two dimensional sky matrix (elevation & cross elevation) with a cell size of 0$`\stackrel{}{\mathrm{.}}`$3$`\times `$ 0$`\stackrel{}{\mathrm{.}}`$3 . The observed chopped signal matrix was deconvolved using an indigenously developed procedure based on the Maximum Entropy Method (MEM) similar to that of Gull & Daniell (1978) (see Ghosh et al. (1988) for details). The FWHM sizes of the deconvolved maps of the point-like source (Saturn) are 1$`\stackrel{}{\mathrm{.}}`$6$`\times `$1$`\stackrel{}{\mathrm{.}}`$9 and 1$`\stackrel{}{\mathrm{.}}`$6$`\times `$1$`\stackrel{}{\mathrm{.}}`$8 in the 143 and 185 $`\mu `$m bands respectively. An optical photometer at the Cassegrain focal plane was used to improve the absolute positional accuracy of the telescope to $``$0$`\stackrel{}{\mathrm{.}}`$5 (Mookerjea et al., 1999).

### 2.2 The HIRES processed IRAS maps

To supplement our balloon-borne observations, we have used the IRAS survey data for all four (12, 25, 60 and 100 $`\mu `$m ) bands for the regions of the sky around IRAS 00494+5617 and 05327-0457. These data were HIRES processed (Aumann, Fowler & Melnyk, 1990) in the Infrared Processing and Analysis Center (IPAC111IPAC is funded by NASA as part of the part of the IRAS extended mission program under contract to JPL., Caltech) for improving the angular resolutions of the raw maps. These maps have been used to obtain flux densities and angular sizes in the four IRAS bands.

## 3 Results

### 3.1 IRAS 00494+5617

Figure Far Infrared Study of IRAS 00494+5617 & IRAS 05327-0457 shows the intensity maps of the region around IRAS 00494+5617 in the two TIFR and the four IRAS bands. Table 1 presents the coordinates and flux densities (in a circle of 5′ dia) at all the six wavelengths. The main peak at 185 $`\mu `$m is shifted by $``$ 1′ compared to the peaks detected at all other wavelengths. The brightest source at 12 $`\mu `$m (IRAS 00492+5614) is too faint in the FIR maps and is not the main source of interest.

The 143 and 185 $`\mu `$m maps show more structural details of the region as compared the 100 $`\mu `$m HIRES map. These maps show extended features (diffuse component) along the east-west direction at low contour levels. The presence of a diffuse component is also noticed in all the HIRES maps. The FIR maps presented here have the shape of a peanut, similar to the shape at 1.3 mm observed by Henning et al. (1994).

In addition to the diffuse component, the 143 $`\mu `$m map shows a secondary peak (S2) towards the west; this is not seen in the maps of other bands. The position and flux density (after correcting for the diffuse emission) of S2 at 143 $`\mu `$m are also presented in Table 1. The limits in the other bands refer to the flux densities measured in a circle (2′ dia) centered at the position of S2 at 143 $`\mu `$m. Although the HIRES maps do not detect S2, the HIRAS map at 60 $`\mu `$m processed with the Groningen IRAS Software Telescope (Henning et al., 1994) with better angular resolution indicates the presence of the same.

Figure Far Infrared Study of IRAS 00494+5617 & IRAS 05327-0457 presents the contour maps of the temperature (T(143/185)) and optical depth ($`\tau _{150}`$) for this region, generated using the method described by Mookerjea et al. (1999). For these maps a dust emissivity, $`ϵ_\lambda `$ $``$ $`\lambda ^1`$ has been assumed. The temperatures determined are accurate to within $`\pm `$ 2 K between 20 and 50 K and within $``$ $`\pm `$ 5 K between 50 and 75 K. The structural details in these maps are highly reliable due to (i) the simultaneity of observations and (ii) conservative processing e.g., 3 pixel by 3 pixel smoothing of the flux densities done prior to the determination of the temperature. The temperature and optical depth maps have been restricted to the regions where the intensities in both the wavebands are substantially higher ($`>`$ 10 times) than the measured noise level. From the temperature map it is seen that there are two regions of enhancement close to the boundary, the highest temperature being close to the boundary facing HD 5005. Further, the temperature distribution seems to be featureless over most of the central region and shows a minimum near the center. In contrast, the peak of the $`\tau _{150}`$ map coincides with the peak at 185 $`\mu `$m and the contours decrease smoothly outwards, as in the intensity map.

### 3.2 IRAS 05327-0457

Figure Far Infrared Study of IRAS 00494+5617 & IRAS 05327-0457 presents the intensity maps of the region surrounding IRAS 05327-0457 in the six wavebands considered in this paper. The positions of the global peak in the 143 & 185 $`\mu `$m maps agree with the IRAS Point Source Catalog (PSC) coordinates within the achieved absolute positional accuracy in the maps. The source is extended along a direction approximately perpendicular to the line joining the cloud to the star 42 Ori. The basic features of the intensity distributions at 143 and 185 $`\mu `$m are similar to the 60 and 100 $`\mu `$m HIRES maps as well as the 158 $`\mu `$m \[C II\] maps by Howe et al. (1991). The emission due to 42 Ori (located to the north of IRAS 05327-0457) is clearly seen in the 12 and 25 $`\mu `$m maps. Table 2 presents the position and flux densities in a circle of 5′ dia centered on the global peak in each of the six wavebands.

In both the 143 and 185 $`\mu `$m maps, a fainter source (P2) is detected towards the south of IRAS 05327-0457. Positionally this matches very well with the coldest sub-millimeter source (CSO 10) detected at 350 $`\mu `$m by Lis et al. (1998). This is associated with the Orion Molecular Cloud -3 and has also been detected at 1.3 mm by Chini et al. (1997). The dust temperature, T(143/185) for P2 has been estimated to be $``$ 19 K assuming $`ϵ_\lambda \lambda ^2`$, the choice of the dust emissivity index being guided by Lis et al. (1998). This low temperature source is not detected in any of the HIRES maps. Table 2 presents the position and flux densities of P2.

Figure Far Infrared Study of IRAS 00494+5617 & IRAS 05327-0457 shows the maps of the dust temperature (T(143/185)) and optical depth ($`\tau _{150}`$) for IRAS 05327-0457, generated using same methods as for IRAS 00494+5617. In the $`\tau _{150}`$ map, the peak has a value of 0.02 and is shifted, north-east towards the direction of the ionizing star. The optical depth decreases monotonically towards the south-west. Higher temperatures are seen towards the south-eastern and western edges; rest of the region is seen to have smoothly varying temperatures, between 25 and 35 K. There is no increase in the dust temperature near the edge facing 42 Ori. The T(143/185) and the $`\tau _{150}`$ map together suggest the presence of (i) embedded sources of heating at the positions of temperature enhancements and (ii) a high density shell facing 42 Ori.

## 4 Discussion

### 4.1 Radiation transfer scheme & Models explored

We have explored the major sources of dust heating based on radiation transfer modeling of the sources. The radiation transfer equations have been solved assuming a 2-point boundary condition for a spherically symmetric cloud of dust and gas. Based on the boundary conditions, we have constructed 3 types of models : $`modelA`$ – the cloud is heated by centrally embedded sources and by an external radiation field due to the average Galactic Interstellar Radiation Field (ISRF) and nearby stars; $`modelB`$ \- the cloud is heated by internal sources and ISRF only; $`modelC`$ – there is no embedded source, the cloud is heated only externally by the ISRF and radiation due to nearby stars. The contribution of the nearby star in the models $`A`$ and $`C`$ has been calculated in the following way: the geometrically diluted stellar radiation intercepted by the cloud surface is estimated and it is then smeared out uniformly over the entire surface of the cloud. This calculated contribution could be a slight overestimate since the absorption in the intervening medium has been neglected. For modeling the observed Spectral Energy Distribution (SED) two types of dust have been considered. The dust properties have been taken from two sources viz., Draine & Lee (1984) (hereafter DL) and Mezger, Mathis & Panagia, (1982) (hereafter MMP). The dust compositions include mixtures of astronomical silicate and graphite. For the DL type of dust grains, size averaged properties are used taking a size distribution of $`n(a)daa^{3.5}da`$ (Mathis, Rumpl & Nordsieck, 1977), with $`a`$ ranging between 0.01 and 0.25 $`\mu `$m. The average Galactic ISRF used for this problem has been taken from Mathis, Mezger & Panagia (1983). The cloud is parameterized by the following physical quantities: Rmax, the outer size of the cloud, Rmin, the radius of the inner dust cavity, $`\tau _{100}`$, the total radial optical depth at 100 $`\mu `$m and the radial dust (hence gas) density distribution ($`r^\beta `$, $`\beta `$=0, -1, -2). The contributions of ISRF and the external source (if any) are kept fixed, while the luminosity of the embedded source, the dust composition, $`\tau _{100}`$ value and radial distribution of the dust density and the physical sizes of the cloud are varied to obtain a good match to observations. The gas to dust ratio has been assumed to be 100:1 by mass. Using this scheme a best fit model matching the observed spectral energy distributions (SED), radial profiles at selected wavelengths and the radio continuum flux are obtained. These radiation transfer calculations have been done using a modified version of the code CSDUST3 (Egan, Leung & Spagna, 1988), which is capable of treating gas and dust in a self-consistent manner (Mookerjea et al., 1999; Mookerjea & Ghosh, 1999).

In the following subsections we explore the tenability of the models $`A`$, $`B`$ and $`C`$ in the light of available observational constraints for individual sources.

### 4.2 IRAS 00494+5617

#### 4.2.1  The main source

The observed SED for the source IRAS 00494+5617 has been constructed using the flux densities presented in Table 1 along with flux densities at 1.3 mm (Henning et al., 1994) and 2.9 mm (Walker, Adams & Lada, 1990). The total observed luminosity is 2.4$`\times `$104 $`L_{}`$ (Carpenter, Snell & Schloerb, 1990) and no radio continuum emission is detected at 8.4 GHz (Tofani et al., 1995) from this cloud. T(143/185) map shows enhanced temperature towards the boundary closer to HD 5005. In order to explain these, the embedded source in radiation transfer models $`A`$ and $`B`$ is assumed to be a cluster of stars. In models $`A`$ and $`C`$ an intercepted luminosity of 6.9$`\times `$103$`L_{}`$ from HD 5005 has been considered. For each type of model it is found that the best fit is obtained for the case of uniform density and grains of DL type. Figure Far Infrared Study of IRAS 00494+5617 & IRAS 05327-0457 shows the observed and the best fit SEDs. Model $`A`$ fits the SED at all wavelengths very well, while model $`B`$ shows some mid-infrared excess. Model $`C`$ fails to reproduce any of the observed flux densities. The predicted radial profiles (as obtained by convolving the model profiles with instrumental PSFs) for the models $`A`$ and $`B`$ are almost similar. Table 3 compares the FWHMs predicted by model $`A`$ with those observed at 25, 60, 100, 143 and 185 $`\mu `$m. The observed FWHMs presented here refer to the FWHMs for the source along its major and the minor axes. We find that model $`A`$ not only fits the observed SED better than model $`B`$ but also achieves an overall consistency by supporting the role of HD 5005 in heating the cloud and fitting the observed FWHMs reasonably well. Table 4 presents the parameters of the best fit model $`A`$ for this source alongwith the uncertainties in the $`\tau _{100}`$ and the embedded luminosity. These uncertainties have been estimated using only the SED fit and keeping the other parameters fixed. From this radiation transfer modeling and the T(143/185) map we conclude that the source IRAS 00494+5617 is primarily heated by one or more embedded sources with non-negligible contribution from HD 5005.

We compare he FWHM predicted from the best fit spherically symmetric model with the observed FWHM along the minor axis. For IRAS 00494+5617 we find that at 100, 143 and 185 $`\mu `$m model predictions for FWHMs match the observed FWHMs (along minor axis) fairly well. However at 25 $`\mu `$m (which traces much hotter dust) the observed FWHM is substantially larger than the model prediction. This trend is also seen at 60 $`\mu `$m, though not to the same extent as at 25 $`\mu `$m. This can be explained by invoking spatially distributed heating sources embedded in clumps. It may be noted that our dust temperature maps also shows the presence of two temperature enhancements. This clumpiness is also substantiated by high resolution C18O observations by Megeath & Wilson (1997), which have revealed that this western fragment of NGC 281 is composed of three clumps. The best fit model $`A`$ obtained above corresponds to a gas density of 2.2$`\times `$104 cm-3 ( for a gas to dust density of 100:1 by mass) and a total mass of 1890$`M_{}`$. Megeath & Wilson (1997) obtain 1080$`M_{}`$ as the lower limit for the total mass in the three clumps. The residual mass (810 $`M_{}`$) as predicted by the model can be present in the form of the inter-clump medium. The projected radii (taking the distance to the source to be 2.9 kpc) of the clumps were given by Megeath & Wilson (1997) assuming the clumps to be spherical. We assume that all three clumps are spherical and are confined within the boundary of the spherical model $`A`$ above. After scaling the clump sizes to our assumed distance of 2.2 kpc, we find that the volume filling factor is 4.5 %. Using the total mass from our model $`A`$ and assuming an interclump medium of uniform density we obtain a value of 4.5$`\times `$103 cm-3 for the number density of the interclump medium. We also calculate the clump-interclump density contrast to be $``$ 80. The resultant clump density of 3$`\times `$105 cm-3 is consistent with the detection of C34S. Although the estimates presented here are based on crude approximations, the values for filling factors of clumps and clump-interclump density contrast obtained here can be called very typical, while comparing with the measured values of such parameters for other clumped sources (Howe et al., 1991).

#### 4.2.2 Diffuse emission

The 143 $`\mu `$m map of this region shows a considerable amount of diffuse emission particularly towards the west. This emission is to a certain extent positionally coincident and structurally similar to that in the 185 $`\mu `$m map. We have looked into the possible energy sources for this emission. The flux densities (in an aperture of 2′ dia centered at 0h 49m 09s +56° 17′ 30″ (1950)) at 143 and 185 $`\mu `$m are 236 and 265 Jy respectively. The bolometric luminosity of the diffuse emission over a diameter of 2′ has been estimated to be $``$ 870 $`L_{}`$. This luminosity cannot be explained due to heating by either the the average ISRF or the star cluster HD 5005. We explore the possibility of radiation actually leaking from the main source (around 00494+5617) and heating the dust in the neighbourhood. Under this assumption we obtain that about 8% of the total luminosity of the embedded source would have to be intercepted by the region where diffuse emission is detected. The separation between the region of diffuse emission and IRAS 00494+5617 is such that the solid angle covered by the diffuse region is adequate to intercept the requisite amount of radiation. This evidence of luminosity leakage further supports the presence of clumpiness discussed in Section 4.2.1.

### 4.3 IRAS 05327-0457

The observed SED for the region around IRAS 05327-0457 has been constructed using the flux densities presented in Table 2 alongwith the IRAS Low Resolution spectrum (LRS) (Chen, Gao & Xiong, 1995). VLA observations of this source by Kutner et al. (1985) have not detected any radio continuum emission at 5 GHz and from the detection limit predict the embedded source to be a ZAMS B2.5 or later type star. Based on their observations Kutner et al. (1985) have suggested that 42 Ori plays the most significant role in fueling the FIR emission from this source. However the T(143/185) map presented here shows no temperature enhancement towards the boundary illuminated by the star 42 Ori. It rather shows an increase from the boundary towards the position of IRAS 05327-0457.

We have applied the radiation transfer models $`A`$, $`B`$ and $`C`$ (see Section 4.1) to IRAS 05327-0457. The embedded source for models of type $`A`$ and $`B`$ is taken to be a star of type B2. In models of types $`A`$ and $`C`$ the intercepted contribution of 600$`L_{}`$ from 42 Ori is considered. Figure Far Infrared Study of IRAS 00494+5617 & IRAS 05327-0457 presents the SEDs predicted from the models $`A`$, $`B`$ and $`C`$ in comparison with the observations. The model outputs shown are the best fits for each type. As in the case of IRAS 00494+5617 uniform density gives the best fit for all three models. Further the MMP type grains lead to better fit than the DL type grains. While the predicted SED from model $`A`$ matches the observations at all wavelengths reasonably well, model $`B`$ predicts fluxes less than the observed values longward of 60 $`\mu `$m. Predictions from model $`C`$ do not the fit the observations at all. The predicted FWHMs (post convolution) are similar for models of type $`A`$ and $`B`$. Table 3 presents the comparison of FWHM sizes from model $`A`$ at different wavelengths with observations. From model $`A`$ the predicted radio continuum emission at 5 GHz is $``$14 mJy which is well above the sensitivity of the VLA observations by Kutner et al. (1985). The reasons for this non-detection could be : (i) the embedded source is a cluster of stars (with masses less than that of a B2 star) with combined luminosity equal to that of a B2 star. This is supported by the detection of multiple stars in the near-infrared (Makinen et al., 1985). They have detected 7 stars having a total luminosity of 2.3$`\times `$103$`L_{}`$ (model $`A`$ total luminosity = 2.8$`\times `$103$`L_{}`$). Our temperature map also indicates presence of two sources, again close to the boundary. (ii) Since the radio continuum emission has a strong n$`{}_{}{}^{2}{}_{e}{}^{}`$ dependence, a smaller gas to dust ratio could explain the reduced emission. (iii) The radio emission could largely be attenuated by free-free absorption. The third possibility proposed would arise, if the gas density distribution very close to the star (r$`<<`$Rmin) is r-2, a typical situation for ionized wind.

Of the models presented here, model $`A`$ explains the observed SED, the radial profile, the T(143/185) map and the role of 42 Ori in heating the cloud most satisfactorily. Table 4 presents the parameters of this model and the estimated uncertainties in the embedded luminosity and $`\tau _{100}`$. From radiation transfer modeling we conclude that it is absolutely necessary to have an embedded source to explain the total emission from IRAS 05327-0457. Clearly 42 Ori cannot be the only heating source since the luminosity intercepted by IRAS 05327-0457 is $``$ 20 % of the observed luminosity. Further we have verified that the temperature profile from radiation transfer calculation using a slab geometry and 42 Ori as the source, does not fit observations.

The observed angular sizes for IRAS 05327-0457 at 100 $`\mu `$m and more prominently at 60 $`\mu `$m are larger than the model values. As in IRAS 00494+5617 this is again consistent with a distribution of embedded hot sources as suggested by the T(143/185) map (Figure Far Infrared Study of IRAS 00494+5617 & IRAS 05327-0457) as well. The apparent reduction in the observed FWHMs at 143 and 185 $`\mu `$m as compared to the FWHM at 100 $`\mu `$m could be due to the effect of sky chopping.

## 5 Summary

This paper presents simultaneous far-infrared mapping observations at 143 and 185 $`\mu `$m of the regions around IRAS 00494+5617 and IRAS 05327-0457. Both sources have been well resolved at these wavelengths. The cold dust source CSO 10 located in OMC-3 (Lis et al., 1998) has also been detected in our map of IRAS 05327-0457. HIRES processed IRAS maps in all 4 bands have been used for comparison. Reliable dust temperature (T(143/185)) and optical depth ($`\tau _{150}`$) maps for these sources are presented. Radiation transfer models have been constructed to explain the observed SEDs, radial profile and the influence of external heating agencies for both these sources. These models provide useful physical parameters for the star forming clouds associated with these sources. For both IRAS 00494+5617 and 05327-0457 it has been demonstrated that although HD 5005 and 42 Ori contribute to their heating respectively, it is necessary to have dominant contribution from embedded sources to explain the observed SEDs. Comparison of predicted radial profiles from models with observations suggest presence of a distribution of embedded clumps in both sources. This is corroborated by spatially resolved peaks in the maps of dust temperature. A simplistic model for the clumpiness of IRAS 00494+5617 has been proposed based on parameters derived from radiation transfer models and previous molecular line observations.

It is a pleasure to thank all members of the Infrared Astronomy Group of T.I.F.R. for their support during laboratory tests and balloon flight campaigns. All members of the Balloon Group and Control Instrumentation Group of the TIFR Balloon Facility, Hyderabad, are thanked for their technical support during the flight. We thank IPAC, Caltech, for providing us the HIRES processed IRAS products. We thank the anonymous referee for the suggestions which have improved the quality of the paper.