text:
# Introductory Remarks on The Extragalactic Infrared Background and its Cosmological Implications

## 1. Fundamental Assumptions

It is useful to begin by recalling the basic assumptions and observations that lead us to the concept of the extragalactic cosmic background radiation, as opposed to radiation surface brightness that may be some highly variable function of position and direction.

Deep counts of objects detected at a broad range of wavelengths, from gamma ray sources to radio galaxies, are close to isotropic across the sky. It is an excellent bet therefore that the integrated radiation from observed sources plus those too faint to be detectable as individual objects also is quite to isotropic. This allows us to think of the local extragalactic radiation background as a function of one variable, the radiation energy density $`u_\nu `$ per interval of frequency $`\nu `$. The fluctuations around the mean as a function of position in the sky are important too, as a measure of large-scale structure, but $`u_\nu `$ is the center of attention in these Proceedings.

The argument for large-scale homogeneity — against a universe with a radial density gradient and us at the center — is less direct but I think persuasive; my review of the considerations is in Peebles (1993). If we are persuaded then we conclude that within our Hubble length space is filled with a near uniform sea of radiation with spectral energy density $`u_\nu `$: the cosmic extragalactic background radiation.

If the propagation of the radiation is described by a metric theory then it satisfies the Liouville or brightness theorem. If the metric describes a homogeneous isotropic spacetime then the geometry is fixed by the expansion factor $`a(t)`$, a function of the proper world time $`t`$ alone, together with the radius of curvature $`a(t)R`$ of sections of constant time, where the comoving radius $`R`$ is a constant. In this spacetime the radiation energy density $`u(t)=𝑑\nu u_\nu `$ integrated over frequency at time $`t`$ is an integral over the history of production and absorption of radiation,

$$u(t)=^t𝑑t^{}j(t^{})(a(t^{})/a(t))^4.$$
(1)
At tine $`t^{}`$ the net rate of production of radiation (emitted minus absorbed) per unit proper volume is $`j(t^{})`$, and $`j(t^{})(a(t^{})/a(t))^3`$ is the rate of production of energy per comoving volume normalized to the time $`t`$ of observation. The remaining factor in the integrand, $`a(t^{})/a(t)=(1+z)^1`$, where $`z`$ is the redshift at the epoch $`t^{}`$ observed at time $`t`$, represents energy lost due to the cosmological redshift.

If spacetime were static, $`a`$ independent of time, equation (1) says $`j`$ could not have been constant: there would have to have been a characteristic time at which star formation commenced. The point, associated with the name Olbers, is not often mentioned now; an edifying discussion is to be found in Bondi (1960). In the classical steady state cosmology (which also is well described by Bondi) the universe is expanding, $`ae^{Ht}`$, where $`H`$ is Hubble’s constant. This makes the integral converge even when $`j`$ is constant, stars forming at a fixed mean rate per physical volume back to the indefinitely remote past. But we know now this is not a viable picture: Cowie and Lilly describe in these Proceedings observations of galaxies and an intergalactic medium at high redshift that are distinctly different from what is observed nearby; the more youthful appearance of objects at high redshift agrees with the interpretation that they are seen closer to the time when the structure we see started forming. In the general relativistic Friedmann-Lemaître model with a classical stress-energy tensor that satisfies $`\rho +3p>0`$ the integral in equation (1) has to have a lower limit, at the singular start of expansion at $`a(t)=0`$. In the eternal inflation scenario (Linde 1990) this unsatisfactory situation is relieved by the return to a steady state philosophy: the lower limit to the integral extends back along our world line to the remote past.

## 2. Cosmic Energy Densities

Let us consider now the interpretation of the radiation background under the standard relativistic cosmology. Evolution after inflation — or whatever produced the initial conditions for the present state of our expanding universe — was accompanied by exchanges of energy among different forms. An accounting of the integrated results of the transactions at the present epoch offers a measure of cosmic evolution, and in particular it informs our interpretation of the infrared background.

The estimates in Table 1 are expressed in units of the Einstein-de Sitter value, $`\rho _c=3H_o^2/(8\pi G)`$, at Hubble constant $`H_o=70\pm 10`$ km s-1 Mpc-1. That is, these numbers are contributions to the cosmological density parameter. The first set of numbers, labeled primeval, are thought to have been fixed by physical processes operating in the early universe, well before stars and galaxies started forming; the second set are estimates of the effects of the formation and evolution of structure on scales ranging from clusters of galaxies down to star remnants.

The accounting in this table accepts the evidence for a Friedmann-Lemaître model that is close to cosmologically flat, the stress-energy tensor being dominated by a term that acts like Einstein’s cosmological constant, $`\mathrm{\Lambda }`$. The next most important term appears to be some form of nonbaryonic dark matter. The baryon density in the third line agrees with the theory of the origin of the light elements in the early universe, with the fluctuation spectrum of the 3 K thermal background radiation — within reasonable-looking uncertainties (eg. Hu et al. 2000) — and with the observational constraints on the baryon budget (Fukugita, Hogan, & Peebles 1998). The baryon entry seems secure to 30% or so, a truly remarkable advance. It is a measure of the state of our subject that the two largest entries are conjectural. The evidence for low pressure dark matter at about the density indicated in the table is compelling if we accept general relativity theory (and hence the inverse square law for gravity); the evidence for $`\mathrm{\Lambda }`$ or its near operational equivalent is strong if we accept the adiabatic CDM model for structure formation. Work in progress promises to establish more tests of general relativity theory applied on the scale of the Hubble length, of the cosmology, and of the theory of structure formation. The results certainly will be searched for potential insights into the enigmatic leading entries in the table.

The entry for relict neutrinos that broke thermal equilibrium with the background radiation at $`z10^{10}`$ assumes the tau neutrino mass is no smaller than 0.03 eV, from atmospheric neutrinos (Super-Kamiokande Collaboration 2000), and no larger than about 1 eV, to avoid undue effect on structure formation (Klypin et al. 1993). The accepted provenance of the 3 K thermal radiation and its related neutrinos — entropy originating near the end of inflation — is conjectural. The observed peak in the angular fluctuation spectrum of this radiation, at the length scale set by decoupling at $`z1000`$, is good evidence this radiation is primeval, present well before the observed stars and galaxies could have started forming.

The meaning of the last of the primeval entries is illustrated by the comparison of two Friedmann-Lemaître models, both containing only cold dark matter particles with the same particle mass, and with the same comoving radius of curvature measured relative to the mean distance between particles. In one model the mass distribution is close to homogeneous; in the other model primeval curvature fluctuations have placed most of the dark matter in gravitationally bound nonrelativistic halos. At a given mean particle number density the mean mass density is smaller in the latter model by the amount of the mean halo gravitational binding energy: the sum of the kinetic energy of proper motion of the particles relative to the general expansion and the gravitational potential energy relative to a homogeneous mass distribution. Thus when the particle number densities are the same the physical radii of curvature are the same in the two models, and the expansion rate, from the relativistic expression $`(\dot{a}/a)^2=8\pi G\rho /3\pm 1/(aR)^2`$, is lower in the model with halos. This gravitational binding energy in the model with halos is primeval in the sense that it is a consequence of the gravitational growth of structure out of given initial conditions; there is no energy transaction when the halos collapse to form bound systems. We have good evidence that this gravitational growth of structure is responsible for the large-scale structure of our universe; most dramatic is the consistency of the angular fluctuation spectrum of the 3 K background with the simple adiabatic CDM model for the gravitational instability picture (eg. Hu et al. 2000 and references therein). In this CDM model radiation pressure suppresses the growth of density fluctuations on scales less than the matter-radiation Jeans length at decoupling, $`10`$ Mpc, lowering the binding energy, and the dissipation of the small-scale pressure waves adds a little entropy. The number in the table for the primeval gravitational binding energy assumes the rms peculiar velocity is comparable to that of the Local Group relative to the 3 K background, $`600`$ km s-1, in matter at the sum of the mass densities in the second and third entries.

Many have commented on the small differences of values of successive entries in this first group, which may in part reflect the anthropic consideration that these parameters can be adjusted so we could not have been here to measure them, may in part result from physical relations to be discovered, and may in part be pure coincidence. Similar remarks apply to the second group of entries, of course.

The first entry in this second group is based on the contribution to the density parameter by the mass bound in compact nuclei of galaxies (from the review by Fabian 2000). If these objects are relativistic black holes, mass may flow into them without the emission of energy in radiation or jets. For our purpose this accretion without emission represents an exchange of small particles for large ones, both gravitationally bound to the massive halos of galaxies, and it has no effect on the accounting in Table 1. The same applies to the kinetic energy of a jet that is dissipated within the galaxy. The factor $`ϵ`$ in the table is the mean fraction of energy in radiation and jets that has left the gravitationally bound halos within which these compact objects live.

Most of the baryons are thought to be in intragroup and intracluster plasma, distributed in about the same way as the dark matter. The gravitational binding energy per unit mass consequently is about the same, and, as we have noted, in the adiabatic CDM model this primeval binding energy appeared without transfer of energy to some other form. The dominant energy exchange in this component of the baryons may be the thermal bremsstrahlung and emission line X-ray radiation from intracluster plasma; less important transactions include the energy exchanges due to galactic winds and the energy exchanges between baryons and dark matter via fluctuations in the gravitational potential.

The formation of baryon-rich spheroids and disks of galaxies requires about the same dissipative production of binding energy per unit mass as does star formation. Here too one can think of subdominant energy transactions: magnetic fields produced by galaxy dynamos and then blown out with the galactic winds, and escaping cosmic rays. The largest transaction likely was the contribution to the cosmic background radiation.

The entry for the binding energy in atomic nuclei heavier than helium is from Fukugita (2000); it uses the baryon budget of Fukugita, Hogan, & Peebles (1998) at $`H_o=70`$ km s-1 Mpc-1, and assumes the heavy element abundance is $`0.4`$ times Solar in the intracluster plasma and $`0.1`$ times Solar in the intragroup plasma. The production of helium in stars assumes the ratio of production of helium and heavier elements is $`1<\mathrm{\Delta }Y/\mathrm{\Delta }Z<4`$.

The energies in the optical/near ir and far ir/sub-mm extragalactic radiation backgrounds are taken from Pei, Fall, & Hauser (1999). These measurements, and their interpretations, are the subject of these Proceedings; a few comments suggested by the entries Table 1 will be noted here.

The optical extragalactic background at $`\lambda `$5000 Å likely is dominated by starlight from low redshift, because the spectra of most galaxies decrease toward the ultraviolet. If so this background is not sensitive to the parameters of the cosmological model, but it is of considerable importance to cosmology as a test for stars in the extreme outer halos of galaxies or in objects with effective radii or surface brightnesses too small to be included in galaxy counts (Arp 1965; Peebles 1971). The measurement of the extragalactic contribution to the light of the night sky has a long history (eg. Dube, Wickes, & Wilkinson 1979; Bernstein 2000). It shows we have not grossly underestimated the density of starlight, but further advances in the measurements will be followed with interest.

The estimate of the energy released by star formation is less than the upper bound on the background radiation between the optical/near ir and far ir/sub-mm peaks, not an observationally promising situation.

The energy released during dissipative contraction to the central baryon-dominated parts of the galaxies could make an interesting contribution to the X-ray background. The contribution by spheroids depends on the mode of formation, whether by dissipative contraction of plasma followed by star formation or by near dissipationless merging of dense star clusters. Wyse (1999) reviews the observational considerations on which might be closer to what happened. If the latter, and if the star cluster building blocks had small individual escape velocities, the radiation accompanying the dissipative assembly of the star clusters could be lost in the uv/optical background. Disks likely formed by dissipative settling, and might be a significant source for the soft X-ray background.

The energy density in the X-ray background is three orders of magnitude down from the estimate of the energy density in black holes in galactic nuclei. Fabian (2000), Hasinger (2000) and others note that even with the $`1+z`$ redshift dimming (eq. ), and a low radiation efficiency, $`ϵ`$, there is an ample budget for significant contributions by AGNs to the optical and sub-mm backgrounds.

About half the extragalactic B-band background radiation comes from the disks of spiral galaxies. The dust in disks absorbs a substantial fraction of the starlight, so the contribution to the energy density in the sub-mm radiation background from starlight absorbed and reradiated by dust would be expected to be comparable to the contribution to the optical/near infrared radiation by unabsorbed starlight. This is in line with the background measurements discussed by Hauser and others at this meeting. The entry for the binding energy in atomic nuclei is comparable to the measured energy in the infrared background. And with Pettini’s (1999) calibration the time integral of the observed star formation history agrees with the present mass density in stars. These three results could be taken to suggest we are getting close to a reconciliation of our accounts. There are some complications to consider, however, as discussed next.

## 3. Star Formation Histories

Michael Fall and Piero Madau discuss in these Proceedings an important observational and conceptual advance, the establishment of a global star formation history. This history can be compared to the record of the evolution of heavy element abundances, and to the resulting extragalactic infrared background. I will present reasons for thinking complications in the history may require a generalization of the star formation history to a function of two variables, world time and environment. The latter might have just three values: normal $`LL_{}`$ galaxy bulges, galaxy disks, and everywhere else. Even with this simple second parameter the assembly of useful observational constraints would be a messy project, but it may be necessary for concordance in the interpretation of a more sophisticated version of Table 1.

To begin, suppose we choose cosmological parameters (and the range of values now under popular discussion is not that broad) so as to fix the world time $`t`$ as a function of redshift. Then the Madau plot (Pei & Fall 1995; Madau et al. 1996) — the observed star formation rate $`dM_{}/dt`$ per unit comoving volume as a function of redshift — could be replotted as the product $`tdM_{}/dt`$ as a function of $`\mathrm{log}t`$ (or, what is equivalent, $`z(dM_{}/dt)(dt/dz)`$ against $`\mathrm{log}z`$). The motivation is the same as for the representation of the energy density of the extragalactic background as $`\nu u_\nu =\lambda u_\lambda `$: in a semi-logarithmic plot the area under the curve is the contribution to the integral over the independent variable per logarithmic interval of the variable. That is, this way to represent the star formation history gives us a picture of when there were significant contributions to the net mass in stars. One would see that about half the observed star formation was relatively recent, at redshift $`z<1`$. Where are these young stars? The answer may well be complicated; I will mention two simple possibilities.

First, the star formation observed at $`z<1`$ might be concentrated somewhere other than in the normal textbook galaxies. Maybe this has something to do with the enigmatic fading faint blue galaxy population; we don’t know where they ended up either. In this scenario most stars may have formed at $`z<1`$ and produced most of the metals and most of the extragalactic infrared radiation. But by assumption the extragalactic infrared background would have no direct relation to the origin and evolution of the populations of normal galaxies.

A second scenario is that most of the stars that formed at $`z<1`$ end up in normal $`LL_{}`$ galaxies. Then the evidence as I understand it is that these stars would have to be subdominant additions to the predominantly older star populations. We might recall that this evidence includes these four points. a. About two thirds of the star mass in the high surface density $`L>L_{}`$ galaxies is in spheroids. The evidence I have been hearing is that most of these spheroid star populations are a good deal older than $`z1`$. b. Galaxy counts as a function of redshift at $`z<1`$ are consistent with near passive evolution of the star populations already present at $`z=1`$. c. It seems to be agreed that the present population of normal $`L>L_{}`$ galaxies by and large were in place, with near familiar morphologies, and close to the present-day comoving number densities, at $`z=1`$. d. The normal-looking $`LL_{}`$ spirals and ellipticals observed at $`z1`$ are close to the Tully-Fisher and fundamental plane dynamical relations, after correction for passive evolution, consistent with the idea that these are full-grown textbook galaxies. I think we must conclude that in this scenario the observed star formation history does not include the main event: the bulk of the star formation would have to be off scale to larger redshift, or maybe present at the redshift range of the Madau plot but heavily obscured by dust and in concentrations small enough to have avoided too many SCUBA detections.

We have a constraint on this hypothetical “main event” from the production of heavy elements. Pei, Fall, & Hauser (1999) estimate that the observed star formation rate integrated from high redshift to $`z2.5`$ is about enough to account for the observed accumulation of mass in heavy elements at this redshift, while Pettini (1999) and Pagel (1999) argue that the observed mass in heavy elements at $`z2.5`$ may actually be significantly less than what might have been expected from the seen amount of star formation up to this epoch. Under the former estimate the “main event” requires a special conjecture: we have to suppose the heavy elements associated with the substantial early star formation are hidden, perhaps in dense clouds. Under the latter estimate heavy elements have to be hidden at $`z2.5`$, so it would not be such a great stretch to imagine the heavy elements associated with the “mail event” are hidden too. Advances in the subtle analysis of the constraints on the star formation history, including the relation between star formation and heavy element production, will be fascinating to follow.

## 4. Concluding Remarks

It has been widely and properly noted that we have significant observational evidence galaxies as we know them formed recently, at $`z<1`$: the energy in the infrared extragalactic background, without large redshift dimming, is comparable to the binding energy of the heavy elements in normal $`L>L_{}`$ galaxies, and the time integral of the seen star formation rate, which implies most stars are young, agrees with the seen mass in stars. But I have argued other lines of evidence indicate the cohort of stars and their heavy elements that formed at $`z<1`$ had to have ended up somewhere other than the high surface brightness parts of normal galaxies, or else are subdominant additions to the normal galaxies. True coincidences happen, of course: we have a sensible case for similar contributions to the sub-mm background by two quite different sources — the radiation from stars and from AGNs — with quite different fractions processed through dust.

The lesson is that untangling the relations among the extragalactic radiation background, the histories of star formation in different environments, and the accumulations of stars and heavy elements, is likely to be a complicated operation. This is hardly surprising, considering that we live in a complicated universe. But the suite of observational evidence in hand is rich, rapidly growing, and fascinating to see sorted out.

### Acknowledgments.

I am grateful to Masataka Fukugita for his help in constructing Table 1, to David Hogg for his help with the astronomy, and to Michael Fall and Piero Madau for their helpful comments on the paper. This work was supported in part by the US National Science Foundation.

## References

Arp, H. C. 1965, ApJ, 142, 402

Bernstein, R. A. 2000, to be published

Bondi, H. 1960, Cosmology (Cambridge University Press), 2nd edition

Dube, R. R., Wickes, W. C., & Wilkinson, D. T. 1979, ApJ, 232, 333.

Fabian, A. C. 2000, astro-ph/0001178

Fukugita, M. 2000, private communication

Fukugita, M., Hogan, C. J., & Peebles, P. J. E. 1998, ApJ, 503, 518

Hasinger, G. 2000, astro-ph/0001360

Hu, W., Fukugita, M., Zaldarriaga, M., & Tegmark, M. 2000, astro-ph/0006436

Klypin, A., Holtzman, J., Primack, J., & Regös, E. 1993, ApJ, 416, 1

Linde, A. D. 1990, Particle Physics and Inflationary Cosmology (Harwood Academic)

Madau, P, Ferguson, H. C., Dickinson, M. E., Giavalisco, M., Steidel, C. C., & Fruchter, A. 1996, MNRAS, 283, 1388

Pagel, B. E. J. 1999, in Galaxies and the Young Universe II, ed. H. Hippelein, astro-ph/9911204

Peebles, P. J. E. 1971, Physical Cosmology (Princeton University Press), p. 59

Peebles, P. J. E. 1993, Principles of Physical Cosmology (Princeton University Press)

Pei, Y. C. & Fall. S. M. 1995, ApJ, 454, 69

Pei, Y. C., Fall, S. M., & Hauser, M. G. 1999, ApJ, 522, 604

Pettini, M. 1999, in Chemical Evolution from Zero to High Redshift, eds. J. Walsh & M. Rosa, astro-ph/9902173

Super-Kamiokande Collaboration 2000, astro-ph/0009001

Wyse, R. F. G. 1999, in The Formation of Galactic Bulges, eds. C. M. Carollo, H. C. Ferguson, & R. F. G. Wyse (Cambridge University Press)


format: text

text:
# On the Astrometric Behavior of Binary Microlensing Events

## 1 Introduction

Searches for Galactic dark matter by monitoring light variations of stars located in the Magellanic Clouds and Galactic bulge caused by gravitational microlensing have been or are being carried out by several groups (MACHO: Alcock et al. 2000; EROS: Afonso et al. 1999; OGLE: Udalski, Kubiak, & Szymański 1997). The searches have been very successful and the total number of candidate events now reaches up to 1,000 (P. Popowski 2000, private communication). Among these events, a notable fraction of events are believed to be caused by binary lenses (Udalski et al. 1994; Bennett et al. 1995; Alard, Mao, & Guibert 1995; Rhie & Bennett 1996; Pratt et al. 1996; Albrow et al. 1996; Dominik & Hirshfeld 1996; Alcock et al. 1999a, 1999b). Detecting binary lensing events is important because one can obtain valuable information about Galactic binary lens population such as the distributions of binary mass ratios and separations.

Due to the difference in the lens system geometry, the light curve of a binary lensing event deviates from the single-peak symmetric one of an event caused by a single point-mass lens. The most important difference in the geometry of a binary lens system from that of a single lens system is the formation of caustics. When a source star crosses the caustics, an extra pair of images appear (or disappear) and its flux is greatly amplified, resulting characteristic sharp spikes in the light curve. If the source does not cross the caustics, on the other hand, the light curve is less perturbed and in most cases has no distinctive features such as the spikes of the caustic crossing event. Di Stefano & Perna (1997) pointed out that the detection efficiency of non-caustic crossing binary lensing events is significantly lower than that of caustic crossing events due to the absence of distinctive features in their light curves combined with small deviations. This was well demonstrated by the recent MACHO result where non-caustic crossing events comprise only $`25\%`$ of the total detected binary lensing events despite their predicted higher frequency of occurrence (Alcock et al. 1999b). Another problem in investigating binaries from microlensing is that even if an event is detected, it is difficult to determine the solution of the binary lens parameters, especially for non-caustic crossing events, due to the severe degeneracy problem (Mao & Di Stefano 1995). The level of degeneracy is so high for most gently perturbed light curves of non-caustic crossing events that many of these events allow fits to models of the binary companion ranging from equal-mass binary to planets (Di Stefano & Perna 1997). Although the level of degeneracy for caustic crossing events is relatively low due to the characteristic spike feature in their light curves, they cannot be totally free from the degeneracy problem (Dominik 1999; Albrow et al. 1999). Therefore, for the acquisition of accurate information about Galactic binaries from microlensing, it will be essential to devise a new method that can effectively detect binary lensing including non-caustic crossing events and can resolve their degeneracy problem.

In addition to photometry, microlensing can also be detected and observed astrometrically. When a source star is lensed by a single point-mass lens, its image is split into two. The typical separation between the two images for a Galactic event caused by a stellar mass lens is on the order of a milli-arcsec. This separation is too small for direct observation of the individual images even with the Hubble Space Telescope. However, if microlensing is observed by using several planned high-precision interferometers from space-based platform, e.g. the Space Interferometry Mission (SIM, http://sim.jpl.nasa.gov), and ground-based large telescopes, e.g. the Keck and the Very Large Telescope, one can measure the astrometric displacements in the source star image light centroid, $`𝜹𝜽_\mathrm{c}`$. Once the centroid shifts are measured and their trajectory is constructed, one can determine the lens proper motion, from which one can better constrain the lens mass compared to the photometrically obtained Einstein ring radius crossing time $`t_\mathrm{E}`$. Astrometric microlensing observation was first proposed nearly simultaneously by Miyamoto & Yoshii (1995), Høg, Novikov, & Polnarev (1995), and Walker (1995) and further analysis considering actual performances of specific instrument has been done by Boden, Shao, & Van Buren (1998). A modified idea of monitoring a nearby high proper motion star working as a lens instead of monitoring a source star for the accurate determination of the lens mass was proposed by Miralda-Escudé (1996) and Paczyński (1998). For a single lensing event, the light centroid traces out an elliptical trajectory during the event (Walker 1995; Jeong, Han, & Park 1999). If an event is caused by a binary, on the other hand, the resulting trajectory of the light centroid (astrometric trajectory) deviates from the elliptical one of the single lensing event (Chang & Han 1999), which is analogous to the deviations in the photometric light curve.

In this paper, we investigate the properties of binary lensing events when they are astrometrically observed by using high precision interferometers. For this, we construct vector field maps of excess centroid shifts, $`\mathrm{\Delta }𝜹𝜽_c(\xi ,\eta )`$, which represent the deviations of the binary lensing centroid shifts from those of a single lensing event as a function of source position $`𝜻=(\xi ,\eta )`$. From the analysis of the maps, we find that the excess centroid shifts are substantial in a considerably large area around caustics. In addition, they have characteristic sizes and directions depending strongly on the source positions with respect to the caustics and the resulting astrometric trajectories have distinctive features which can be distinguished from the deviations caused by other reasons. We classify the types of the deviations and investigate where they occur.

## 2 Basics of Binary Lensing

When lengths are normalized to the combined Einstein ring radius, the lens equation in complex notations for a binary lens system is represented by

$$\zeta =z+\frac{m_1}{\overline{z}_1\overline{z}}+\frac{m_2}{\overline{z}_2\overline{z}},$$
$`(1)`$
where $`m_1`$ and $`m_2=1m_1`$ are the mass fractions of the individual lenses, $`z_1`$ and $`z_2`$ are the positions of the lenses, $`\zeta =\xi +i\eta `$ and $`z=x+iy`$ are the positions of the source and images, and $`\overline{z}`$ denotes the complex conjugate of $`z`$ (Witt 1990). The combined Einstein ring radius is related to the total mass and the location of binary the lens by

$$r_\mathrm{E}=\left(\frac{4GM}{c^2}\frac{D_{ol}D_{ls}}{D_{os}}\right)^{1/2},$$
$`(2)`$
where $`D_{ol}`$, $`D_{ls}`$, and $`D_{os}`$ represent the separations between the observer, lens, and source star, respectively. For a binary lensing event, the number of images is either three or five depending on the source position. The amplifications of the individual images are given by the Jacobian of the transformation (1) evaluated at the image position, i.e.

$$A_i=\left(\frac{1}{|\mathrm{det}J|}\right)_{z=z_i};\mathrm{det}J=1\frac{\zeta }{\overline{z}}\frac{\overline{\zeta }}{\overline{z}}.$$
$`(3)`$
The caustic refers to the source position on which the amplification of a point source event becomes infinity, i.e. $`\mathrm{det}J=0`$. The total amplification is given by the sum of the amplifications of the individual images, i.e. $`A=_iA_i`$. Since the position of the image centroid equals to the amplification weighted mean position of the individual images, the source image centroid shift with respect to the unlensed source position is given by

$$𝜹𝜽_\mathrm{c}=\frac{\underset{i}{}A_i𝒛_i}{_iA_i}𝜻,$$
$`(4)`$
where the positions of the individual images $`𝒛_i`$ are obtained by numerically solving the lens equation (1).

## 3 Properties of Various Types of Binary Lensing Events

The size of caustics depends strongly on the binary separation $`\mathrm{}`$ (normalized by the combined angular Einstein ring radius $`\theta _\mathrm{E}=r_\mathrm{E}/D_{ol}`$), causing strong dependency of the caustic crossing probability on the binary separation. Caustic crossing is optimized when the binary sepration is equivalent to the angular Einstein ring radius. Han (1999) estimated that when the binary separation is in the range $`0.7\mathrm{}1.8`$, the caustic crossing probability is greater than $`50\%.`$111He estimated the probability under the definition of a binary lensing event as ‘a close lens-source encounter within the combined Einstein ring radius with its center at the center of mass of the binary lens system’. However, as the separation further increases or decreases, the caustic crossing probability decreases rapidly. Therefore, most non-caustic crossing events are caused by binaries with separations either substantially smaller (narrow binaries) or larger (wide binaries) than $`\theta _\mathrm{E}`$. In this work, we focus on the astrometric properties of events caused by these types of binary lenses because they are much more common.

### 3.1 Narrow Binary Lensing Events

If an event is caused by a narrow binary, the lens system forms 3 separated closed caustics. One of them is located near the center of mass (the central caustic) and has a diamond shape with 4 cusps. The other two are located away from the center of mass (the outer caustics) and each of them has a wedge-like shape with 3 cusps. The outer caustics are symmetric with respect to the binary axis and are located on the heavier lens side with respect to the center of mass. Compared to the central caustic, the outer caustics are much smaller (see Figure 1). If the binary separation is very small, both the central and outer caustics become very small and the geometry of the lens system mimics that of a single lensing event with a mass equal to the total mass of the binary and located at the center of mass of the binary (Gaudi & Gould 1997). For a narrow binary lensing event, we, therefore, define the excess centroid shift by

$$\mathrm{\Delta }𝜹𝜽_\mathrm{c}=𝜹𝜽_{\mathrm{c},\mathrm{b}}𝜹𝜽_{\mathrm{c},\mathrm{s}},$$
$`(5)`$
where $`𝜹𝜽_{\mathrm{c},\mathrm{b}}`$ and $`𝜹𝜽_{\mathrm{c},\mathrm{s}}`$ represent the centroid shifts of the binary and the corresponding single lensing event, respectively.

In Figure 1, we present the vector field map of the excess centroid shifts, $`\mathrm{\Delta }𝜹𝜽_\mathrm{c}(\xi ,\eta )`$, for an example narrow binary lens system with $`\mathrm{}=0.5`$ and $`q=0.2`$. The positions in the map are chosen so that the center of mass of the binary is at the origin. Both lenses are located on the $`\xi `$ axis and the heavier lens is to the right. The circle drawn with a long-dashed line represents the combined Einstein ring and the closed figures drawn by thick solid line are caustics. Also drawn on the map are the contours of the absolute value of the excess centroid shift, $`\mathrm{\Delta }\delta \theta _\mathrm{c}`$. The contours are drawn at the levels of $`\mathrm{\Delta }\delta \theta _\mathrm{c}=0.02\theta _\mathrm{E}`$ (dotted line) and $`0.2\theta _\mathrm{E}`$ (solid line). The angular Einstein ring radius of a Galactic bulge event caused by a stellar mass lens with $`M0.3M_{}`$ located at the half way between the observer and the source, i.e. $`D_{ol}/D_{os}=0.5`$, is $`\theta _\mathrm{E}300`$ $`\mu `$-arcsec. Therefore, these thresholds correspond to $`\mathrm{\Delta }\delta \theta _\mathrm{c}6`$ $`\mu `$-arcsec and $`\mathrm{\Delta }\delta \theta _\mathrm{c}60`$ $`\mu `$-arcsec, which are respectively related to the astrometric precisions of the SIM (Unwin et al. 1997) and the interferometers on the large ground-based telescopes (Colavita et al. 1998; Mariotti et al. 1998). To better show $`\mathrm{\Delta }𝜹𝜽_\mathrm{c}(\xi ,\eta )`$ near the caustic regions, we expand map of the regions and presented in Figure 2: upper panel for the central caustic region and lower panel for the outer caustic region. In Figure 3, we also present the light curves (upper 6 panels) and the astrometric trajectories (lower 6 panels) for several example non-caustic crossing events whose source trajectories are marked in Figure 1a (straight lines with arrows). Note that the number in each panel of Fig. 3 matches with the number of the corresponding source trajectory in Fig. 1.

From the analysis of $`\mathrm{\Delta }𝜹𝜽_\mathrm{c}(\xi ,\eta )`$ in various regions outside the caustics222 Since our major interest is the astrometric properties of non-caustic crossing events, in this section we limit our discussion only about excess centroid shifts in the regions outside of caustics. The properties of $`\mathrm{\Delta }𝜹𝜽_\mathrm{c}`$ in the region inside of caustics will be discussed in § 3.3  and the resulting light curves and astrometric trajectories, one finds the following patterns of the excess centroid shifts. First, as already known by Chang & Han (1999), the amount of excess centroid shifts is substantial in a considerably large area around caustics, c.f. the excess amplification map of Gaudi & Gould (1997), implying that astrometric binary lensing observation is an efficient method in detecting non-caustic crossing events caused by narrow binary lenses.333 It would be useful to estimate the average amount of deviation expected from a typical Galactic event. For this estimation, it is required to know not only the physical distributions of Galactic binaries and their mass function but also the distributions of their mass ratios and separations. However, our knowledge about these quantities are very poor and thus we just compare the expected deviation for events caused by an example binary lens system.  The excess centroid shifts are especially large in the regions along the binary axis and the lines connecting the central and outer caustics. In the region along the binary axis, the excess centroid shifts diverges from the on-binary-axis cusps of the central caustic. Then, if a source passes through this region (e.g. the source trajectories 4 and 5), the resulting astrometric trajectory becomes convex compared to the elliptical one of the corresponding single lensing event. By contrast, in the regions along the lines connecting the central and outer caustics, the excess centroid shifts converge towards the off-binary-axis cusps of the central caustic. Then, the astrometric trajectory of an event with a source trajectory passing through this region (e.g. the source trajectories 1, 2, and 3) will become concave with respect to that of the corresponding single lensing event or even twisted. Due to these various types of deviations, the astrometric trajectories of non-caustic crossing binary lensing events have distinctive feaures. By contrast, the light curves in most cases have no strong characteristic features.

### 3.2 Wide Binary Lensing Events

If an event is caused by a wide binary, on the other hand, the lens system forms 2 separate closed caustics, which are located close to the individual lenses along the binary axis. Each caustic has a wedge-like shape with 4 cusps (see Figure 4). For an event caused by these systems, the observed light curve and the astrometric trajectory appear to be the superpositions of those resulting from two events in which the component lenses behave as if they are two independent single lenses (Di Stefano & Mao 1996; Di Stefano & Scalzo 1999; An, Han, & Park 2000). Therefore, we define the excess centroid shifts of a wide binary lensing event by

$$\mathrm{\Delta }𝜹𝜽_\mathrm{c}=𝜹𝜽_{\mathrm{c},\mathrm{b}}𝜹𝜽_{\mathrm{c},\mathrm{s},\mathrm{comb}},$$
$`(6)`$
where the combined centroid shift $`𝜹𝜽_{\mathrm{c},\mathrm{s},\mathrm{comb}}`$ represents the vector sum of the centroid shifts caused by the individual single lenses.

In Figure 4, we present excess centroid shift map for an example wide binary lens system with $`\mathrm{}=2.0`$ and $`q=0.5`$. In the figure, the two circles in the map drawn with a long-dashed line represent the Einstein rings of the individual single lenses. Besides that there are two Einstein rings, the positions and notations are same as those in Fig. 1. Maps in the regions around caustics are presented in Figure 5. In Figure 6, we also present the light curves and astrometric trajectories for several example events resulting from the source trajectories marked in Fig. 4.

From the figures, one finds that $`\mathrm{\Delta }𝜹𝜽_\mathrm{c}`$ is also large for the wide binary lens system in a considerably large area around caustics and the resulting astrometric trajectories have distinctive features. Especially large excess centroid shifts occur in the region along the binary axis. One also finds that the general patterns of $`\mathrm{\Delta }𝜹𝜽_\mathrm{c}`$ in the region around each caustic is similar to those in the region around the central caustic of the narrow binary lens case. That is, the excess centroid shifts diverges from the on-bianry-axis cusps and converges towards the off-binary-axis cusps. One major difference is that in the region between the right-side on-binary-axis cusp of the left caustic and left-side on-binary-axis cusp of the right caustic, the combination of the two diverging excess centroid shift vectors from the individual cusps makes $`\mathrm{\Delta }𝜹𝜽_\mathrm{c}`$ directed almost vertically from the binary axis. As a result, the centroid shift trajectory of an event passing through this region (e.g. part of the trajectory 3 around the moment of the source star’s closest approach to the center of mass) becomes convex compared to the astrometric trajectory of the combined single lensing events. Another interesting deviations occur when the source passes closely both lenses. For this case, the astrometric trajectory will have double bumps or twisted loops (e.g. the trajectories 1, 2, and 3), which characterize the event is caused by a wide binary.

### 3.3 Caustic Crossing Evenst

When binary separation becomes equvalent to $`\theta _\mathrm{E}`$, the separated caustics merge together and form a single large closed caustic encompassing the center of mass of the binary. The caustic has 6 cusps and two of them are located on the binary axis (see Figure 7). Since the size of the caustic is maximized at around this binary separation, events caused by these binaries have high chance of caustic crossings.

To see the astrometric properties of caustic crossing events, we present in Figure 7 the vector field map of a binary lens system with $`\mathrm{}=1.0`$ and $`q=0.5`$. We also present in Figure 8 two example astrometric trajectories (lower two panels) whose source trajectories are marked in the map. Although this binary lens system cannot be well approximated either by a single lens nor by the superposition of the two independent single lenses, we compute the excess centroid shifts by using equation (5). From the figures, one finds the general pattern of the deviation vector in the region inside of the caustic is that $`\mathrm{\Delta }𝜹𝜽_\mathrm{c}`$ converges towards the on-binary-axis cusps and diverges from the off-binary-axis cusps. Note that this pattern is exactly opposite to that of the excess centroid shifts in the region outside of the caustic. As a result, the two adjacent excess centroid shift vectors just inside and outside of the caustic line have opposite directions. When a source crosses the caustic, then, the deviations in the astrometric trajectory becomes very big not only because the absolute value of the centroid shift is very large near the caustics but also the direction of $`\mathrm{\Delta }𝜹𝜽_\mathrm{c}`$ is reversed during the caustic crossing. We note that these patterns of $`\mathrm{\Delta }𝜹𝜽_\mathrm{c}`$ in the region very close to and inside of caustics apply also to caustic regions of the narow and wide binary lens systems (see the lower panels of Figure 1 and 3). The only exception is the outer caustic of the narrow binary lens system, in which the direction of $`\mathrm{\Delta }𝜹𝜽_\mathrm{c}`$ is not reversed during caustic crossings (see the lower right panel of Fig. 2).

## 4 Confusion with Other Types of Deviations

In previous section, we illustrate various forms of deviations in the astrometric trajectories caused by the lens binarity. Besides the lens binarity, however, the trajectory can also be distorted due to various other reasons. First, the astrometric trajectory can be distorted by the change of the observer’s location during the event caused by the orbital motion of the Earth around the Sun: parallax effect (Paczyński 1998; Han 2000). Second, if a star not participating in the lensing process is located too close to the source star (not the lens) to be resolved even with a high resolution interferometer, such as a close companion of the source star, the astrometric trajectory can also be distorted: blending effect (Han & Kim 1999). Third, if the size of the source star is not negligible compared to the Einstein ring radius, the astrometric trajectory is distorted: finite source effect (Mao & Witt 998).

However, the distortions caused by these effects have different forms from those affected by the lens binarity. The trajectory affected by the blending effect has characteristic linear component caused by the displacement towards the blended star \[see the example trajectories in Figure 4 of Han & Kim (1999)\], and thus can be distinguished from the trajectories affected by the lens binarity. The deviation caused by the parallax effect is periodic with a known period of a year, allowing one to notice the effect. The devations caused by the finite source effect is somewhat confusing because it can produce concave and twisting distortions similar to those of binary lensing events. However, the resulting astrometric trajectories of events affetced by the finite source effect are exactly symmetric with respect to the semi-minor axis of the unperturbed elliptical trajectory. On the other hand, the astrometric trajectories of binary lensing events are asymmetric in general, and thus one can distinguish the two types of astrometric deviations.

## 5 Summary

The findings about the astrometric properties of binary lensing events from the analysis of the vector field maps of excess centroid shifts for various types of binary lens systems are summarized as follows.

1. The lens binarity causes large deviations in the astrometric trajectory from the elliptical one of a single lensing event. The excess centroid shifts are substantial even in the regions considerably away from the caustics, implying that the distortions will be significant even for a significant fraction of non-caustic crossing binary lensing events.
2. Depending on the locations with respect to lens caustics, the excess centroid shifts have very characteristic sizes and directions. As a source passes through these regions, the resulting astrometric trajectories have distinctive features, which are categorized by concave, convex, and twisting distortions.
3. The binary-lens induced deviations in the astrometric trajectory can be distinguished from those caused by other effects, such as the parallax, blending, and finite source effects.

Therefore, astrometric microlensing observations will be an efficient method to detect events caused by binary lenses and to constrain their parameters with accuracy.

Acknowledgements: This work was supported by grant KRF-99-041-D00442 of the Korea Research Foundation.

Figure 1: Vector field map of the excess centroid shifts for an example narrow binary lens system with a separation $`\mathrm{}=0.5`$ and a mass ratio $`q=0.2`$. The positions are chosen so that the center of mass is at the origin. Both lenses are located on the $`\xi `$ axis and the heavier lens is to the right. The circle drawn with a long-dashed line represents the combined Einstein ring. The three closed figures drawn by thick line are caustics. Also drawn are the contours of the absolute values of the excess centroid shift, $`\mathrm{\Delta }\delta \theta _\mathrm{c}`$. The contours are drawn at the levels of $`\mathrm{\Delta }\delta \theta _\mathrm{c}=0.02\theta _\mathrm{E}`$ (dotted line) and $`0.2\theta _\mathrm{E}`$ (solid line). The straight lines with arrows represent the source trajectories whose resulting light curves and astrometric trajectories are presented in Fig. 3.

Figure 2: Vector field map of the excess centroid shifts in the region near the caustics of the same lens system as in Fig. 1. The upper and the lower panels are for the central and the outer caustic regions, respectively

Figure 3: Light curves (upper 6 panels) and astrometric trajectories (lower 6 panels) for several example non-caustic crossing narrow binary lensing events whose source trajectories are marked in Fig 1. The number in each panel matches with the number of the corresponding source trajectory in Fig. 1. The dotted curve in each of the lower set of panels represents the astrometric trajectory of a single lensing event with a mass equal to the total mass of the binary and located at the center of mass of the binary.

Figure 4: Vector field map of the excess centroid shifts for an example wide binary lens system with a separation $`\mathrm{}=2.0`$ and a mass ratio $`q=0.5`$. The two circles drawn by a long-dashed line represent the Einstein rings of the individual single lenses. The positions and notations are same as those in Fig. 1.

Figure 5: Vector field map of the excess centroid shifts in the region near the caustics of the same lens system as in Fig. 4.

Figure 6: Light curves (upper 6 panels) and astrometric trajectories (lower 6 panels) for several example non-caustic crossing wide binary lensing events whose source trajectories are marked in Fig. 4. The notations are same as those in Fig. 3.

Figure 7: Vector field map of the excess centroid shifts for a binary lens system with $`\mathrm{}=1.0`$ (upper panel) and $`q=0.5`$. The notations are same as those in Fig. 1.

Figure 8: The astrometric trajectories of caustic crossing events whose source trajectories are marked in Fig 7. The dotted curves are the corresponding single lensing events.


text:
# Exploring the gravitationally lensed system HE 1104$-$1805: VLT spectroscopy of the lens at z=0.729 Based on observations collected during ESO Program 65.O-0566(A)

## 1 Introduction

HE~1104$-$1805 was discovered as part of the Hamburg/ESO Quasar Survey and was first identified as a gravitational lens candidate by Wisotzki et al. (Wizot93 (1993)). It consists of two lensed images of a radio-quiet quasar at $`z`$= 2.319 that are separated by $``$3.2″. The lensing galaxy was discovered from ground based near-IR (Courbin, Lidman & Magain, Courbin98a (1998)) and HST optical observations (Remy et al. Rem1998 (1998); hereafter R98) and is 1″ from the brighter quasar image (component A).

From a spectrophotometric monitoring program that lasted several years, Wisotzki et al. (Wizot98 (1998); hereafter W98) measured a time delay of 0.73 years between the two components. To convert this into an estimate of the Hubble Constant, one needs to determine the geometry of the system (astrometry and lens and source redshifts) and to model the mass distribution of the lens. Several mass models for the lens have been published (W98; R98; Courbin et al. Courbin2000a (2000); hereafter C2000; Lehar et al. Lehar00 (2000)) and precise astrometry from HST images is available. Since the source redshift is known, the remaining unknown is the redshift of the lens.

In their discovery paper, Wisotzki et al. (Wizot93 (1993)) noted that the continuum of component A was considerably harder than that of component B, and that the broad emission lines of the two components were identical, only scaled by a factor of 2.8. This was still the case during the spectrophotometric monitoring carried out by Wisotzki et al. (Wizot98 (1998)) and in IR spectra covering the 1 to 2.5 $`\mu `$m range (C2000). This has been interpreted as microlensing of the continuum emitting region of the quasar in component A, which is twice as close to the lensing galaxy as component B. Presumably, the broad line region is too large to be affected. It may be possible to use this information to gain insight into the geometries and sizes of the continuum emitting and broad line regions; however, detailed predictions require the lens redshift.

From high resolution spectroscopic observations (Smette et al. Smette95 (1995); Lopez et al. Lopez99a (1999)), several metallic absorption line systems have been detected: $`z=0.517,0.728,1.280,1.320,1.662,1.860,2.220`$ and $`z=2.300`$. The systems at $`z=0.728`$, $`z=1.320`$ and $`z=1.662`$ contain lines that are mostly detected in component A, the component that is closest to the lens. At one time or another, each system has been individually proposed to be the lens.

Despite its importance, the redshift of the lensing galaxy has proved elusive. Apart from the many unpublished attempts to measure it, the redshift of the lens has been estimated by indirect means. From IR and optical photometry, C2000 gave a redshift in the range $`z=0.8`$ to $`z=1.2`$ while from the time delay and a model of the lens, W98 estimated $`z=0.79`$. From the position of the lens on the fundamental plane, Kochanek et al. (Koe20 (2000)) derived $`z=0.77`$. Although model dependent, the two latter estimates prove to be very close to the truth.

In this paper we describe the successful measurement of the lens redshift, $`z=0.729`$, using the spectral deconvolution method described by Courbin et al. (Courbin2000b (2000)).

## 2 Observations-Reductions

The observations were taken with FORS2 on Kueyen (VLT/UT2) at the Cerro Paranal Observatory on April 1st, 2000, which also happened to be the first scheduled observing night for this telescope, and consisted of three 1080 second exposures with the G200I grism and the high resolution collimator. This gives a scale of 0.1″ per pixel in the spatial direction and a scale of approximately 4Å per pixel in the spectral direction. The observing conditions were clear and the external seeing varied between 0.5″ and 0.9″.

We used the movable slits in the FORS2 focal plane to select targets and set the slit width. One slit was placed on HE~1104$-$1805 and aligned along the two quasar components, three slits were placed on field stars that were used to determine the PSF required for the deconvolution (see below) and to correct for telluric absorption. Used in this mode, the slits of FORS2 have a fixed length. This is less flexible than punched or laser designed masks, but it was adequate for our experiment.

The slit width was set to 1″ , which is larger than the FWHM of the images. For the deconvolution to work well, it is important that the targets are not too heavily occulted by the slit edges, since this can lead to differential occultation between the PSF stars and the quasar images and hence to a less accurate deconvolution. This can be minimized by carefully centroiding the stars in the slits and by having slits that are significantly wider than the PSF. Additionally, preparing the mask from an image that was taken with the same instrument, as was done in this case, also minimises differential alignment errors.

Instrumental signatures were removed in the standard manner using IRAF. The bias was removed by subtracting a constant from each frame and the pixel-to-pixel sensitivities were normalised with screen flats. The arc frames were used to correct for slit distortion and to set the wavelength scale. All spectra were rebinned to the same linear wavelength scale. The position of the night sky lines and the absorption features in component A were used to check the accuracy of the wavelength calibration, which was found to be around $`3\AA `$.

## 3 Spectral deconvolution: deblending the lens and source

The lensing galaxy is four to five magnitudes fainter than component A of HE~1104$-$1805 and only one arc second away. Even with excellent seeing, the spectrum of the lens is heavily contaminated by that of the quasar and needs to be extracted with sophisticated techniques.

The technique implemented by Courbin et al. (Courbin2000b (2000)) spatially deconvolves spectra and decomposes them into point-sources and extended components. It is therefore very well suited to the present problem: extracting the spectrum of a faint extended lens galaxy close to a bright quasar spectrum. As is the case with image deconvolution, a reference PSF is necessary. Out of the three PSF stars that were available, two were unfortunately too faint. In order to build the PSF spectrum, we therefore used only one PSF star (which was slightly fainter than the quasar itself) in combination with the spectrum of the bright quasar image itself. Only half of the spatial profile of the quasar was actually used, the one unaffected by light contamination from the lensing galaxy. The deconvolved spectrum has a pixel size of 0.05″ while the undeconvolved data have a pixel size of 0.1″. The deconvolved spectrum is displayed in Fig. 1.

As the PSF stars were observed at exactly the same time and under the same conditions as the lensed quasar, they could also be used to remove telluric features. The spectra of the PSF stars were combined and then divided directly into the spectra of the lens and the two components of the quasar. Before division, the continuum of the PSF stars was normalised with a high-order function and obvious absorption features, such as $`\mathrm{H}\alpha `$ and the Calcium triplet at 8600 Å, were removed. As we do not know the spectral type of the PSF stars, weaker features may be present. Ideally, one should choose hot, featureless stars to make such a correction. Despite this, the correction works rather well.

The 1-D spectrum of the lens, extracted directly from the 2-D deconvolved image (extended component only), is shown in Fig. 2. An unsmoothed version and a smoothed version shifted by 0.03 units on the vertical scale are shown in the lower half of the figure. Also shown, but reduced by a factor of 20, is the spectrum of component A. The Calcium H and K absorption lines of the lensing galaxy are clearly detected and the G-band is marginally detected. The redshift was measured by cross correlating the lens spectrum with the elliptical template from the Kinney-Calzetti spectral atlas (Kinney et al. Kinney96 (1994)), which has been plotted (third curve from the bottom) with a shift of 0.05 units in Fig. 2. For this we used the IRAF task RVSAO v2.0 (Kurtz and Mink 1998). The redshift is $`z=0.729\pm 0.001`$, where the error incorporates the internal error reported by the cross-correlation routine and the error in the wavelength calibration. An alternative measure of the accuracy is given by the r-statistic of Tonry and Davis (Tonry79 (1979)), which is the signal-to-noise ratio of the main peak in the cross correlation. We find that $`r=4.1`$, which indicates that the redshift is reliable (Kurtz & Mink Kurtz98 (2000)).

The redshift is in good agreement with the redshift of the multiple absorption line system at a redshift of $`z=0.728`$ (Smette et al. Smette95 (1995); Lopez, Lopez99b (1999)). This absorption line system is made up of many sub-components spread over a range of redshifts ($`\mathrm{\Delta }z0.001`$). Although it is not possible to say if the absorption line system and the lens are the same object, the two systems are probably dynamically associated.

As the quality of the PSF was not optimal, the quasar spectra from the deconvolved images are noisier than what one would calculate from the photon noise only. Therefore, the extraction of the quasar spectra was not done from the deconvolved spectrum but with narrow apertures from the undeconvolved spectrum. Since the galaxy is much fainter than either component of the quasar and since it is not spatially coincident, the contamination of the quasar spectra by the lensing galaxy is negligible.

The spectra of the two components of the quasar are shown in Fig. 3. In this figure we plot the data before and after the telluric features have been removed. The A and B atmospheric bands are removed very well and the FeII absorption lines at $`z=1.66`$ become visible. Additionally, the MgII emission line, which is almost obscured by the strong telluric feature at 9300Å, is recovered.

As in W98 and C2000, we also plot the difference between the spectrum of component A and a scaled version of the spectrum of component B, that is $`f_\lambda (A)c.f_\lambda (B)`$, where, $`c`$, is set to 2.8. The difference is plotted in Fig. 4, where we plot it before and after telluric correction. Apart from a slight excess at 6390Å and the metallic absorption lines in component A, the broad lines of the quasar disappear and only the excess continuum from component A remains. The strongest absorbers are now clear and are marked with the vertical dotted lines. The wavelengths of the main absorption lines and the redshifts of the corresponding absorption line systems are given in Smette et al. (Smette95 (1995)).

## 4 Discussion

The aim of our observations were modest: to measure the redshift of the lensing galaxy in HE~1104$-$1805, one of the ever growing list of lenses with known time delays. This was possible with only one hour of VLT time during average weather conditions. The prospects of measuring the lens redshifts in other systems are therefore very good, at least for double systems, where aligning the slit simultaneously on the quasar images and on the lens is straightforward.

With the time delay and the lens and source redshifts known and several mass models published, it is relatively straight forward to derive the Hubble constant. W98 present two models to describe the lens: a singular isothermal sphere with external shear and a singular isothermal ellipsoid without external shear. For the former (the model presented by R98 is very similar), one derives $`H_0=52`$ km/s/Mpc and for the latter, one derives $`H_0=59`$ km/s/Mpc. We have set $`(\mathrm{\Omega }_M,\mathrm{\Omega }_\mathrm{\Lambda })=(0.3,0.7)`$. C2000 also present two models: a singular ellipsoid without external shear and a two component model which consists of a singular isothermal ellipsoid to represent the lensing galaxy and a more extended component representing a galaxy cluster centered on the lensing galaxy. The first model leads to $`H_0=62`$ km/s/Mpc, but the second leads to values for the Hubble Constant which are a factor of 1.5 to 2 lower.

HE~1104$-$1805 is somewhat unusual in that the lensing galaxy is closest to the brighter component. It has been pointed out by several authors that the lens probably does not consist of a single galaxy (W98, R98, C2000; Lehar et al. Lehar00 (2000)). Additional extended mass is required to match all the observational constraints: the position of the two images relative to the lens, the magnification ratio of the images, and the lens orientation and ellipticity. Deep images in the near-IR may be able to reveal the distribution of this additional mass.

As the lens appears to be associated with the absorption system seen in the quasar spectrum at $`z=0.728`$, one could use the velocity spread in this system to constrain the lens velocity dispersion. However, this assumes that the absorption line system and the lens are the same object or, at the very least, in the same potential well. Without further supporting evidence that links the absorption line system to lens, it is premature to use the velocity spread as a constraint on the mass models.

The lens redshift estimates of W98 and Kochanek et al. (Koe20 (2000)) have proven to be accurate, whereas the estimate of C2000 from IR and optical photometry is less so. The optical photometry is in good agreement with the VLT spectrum; however, if the lens is an elliptical galaxy, then the IR points are over-estimated by 0.3 magnitudes. The cause of this discrepancy requires further investigation.

###### Acknowledgements.

It is a pleasure to thank Thomas Szeifert for his expert support at the telescope. Frédéric Courbin acknowledges financial support from the Chilean grant FONDECYT/3990024. Additional support from the European Southern Observatory, through CNRS/CONICYT grant 8730 “Mirages gravitationnels avec le VLT: Distribution de matière noire et contraintes cosmologiques”as well as through the TMR LENSNET network (ERBFMRXCT97-0172), is also gratefully acknowledged.


text:
# Simulation accuracy of long range muon propagation in medium: analysis of error sources

## I INTRODUCTION

In order to obtain muon spectrum and value for integral muon flux at large depths where existing and planned neutrino telescopes are located (one to several kilometers of water equivalent) one has to perform several steps starting with interactions of primary cosmic rays in atmosphere. The last step of this procedure is transportation of muons from sea level down to detector location. It can be done both with analytical or semianalytical technique and Monte Carlo (MC) simulation which directly takes into account that muon energy loss represents an essentially stochastic process. MC codes for muon transportation have to use some simplifications to optimize equilibrium between reasonable computation time, desirable accuracy and neccessary statistics.

The usual approximation when calculating the muon propagation through a medium is the following: muon interactions with comparatively large energy transfers, i.e., when fraction of energy lost $`v=\mathrm{\Delta }E/E`$ exceeds some value $`v_{cut}`$ 10-3–10-2, are taken into account by direct simulation of $`vv_{cut}`$ for each single interaction according to shape of differential cross sections $`d\sigma (E,v)/dv`$ (these interactions lead to “stochastic” energy loss or SEL) while the part of interaction with relatively small $`v`$ is treated by the approximate concept of “continuum” energy loss (CEL), i.e., using the function $`[dE(E)/dx]_{CEL}`$ which is obtained from formula

$`\left[{\displaystyle \frac{dE}{dx}}(E)\right]_{CEL}`$ $`=`$ $`{\displaystyle \frac{N_A}{A}}E\left({\displaystyle \underset{v_{min}^\gamma }{\overset{v_{cut}}{}}}{\displaystyle \frac{d\sigma ^\gamma (E,v)}{dv}}v𝑑v\right)+{\displaystyle \frac{N_A}{A}}E\left({\displaystyle \underset{v_{min}^{e^+e^{}}}{\overset{v_{cut}}{}}}{\displaystyle \frac{d\sigma ^{e^+e^{}}(E,v)}{dv}}v𝑑v\right)`$ (1)
$`+`$ $`{\displaystyle \frac{N_A}{A}}E\left({\displaystyle \underset{v_{min}^{pn}}{\overset{v_{cut}}{}}}{\displaystyle \frac{d\sigma ^{pn}(E,v)}{dv}}v𝑑v\right)+\left[{\displaystyle \frac{dE}{dx}}(E)\right]_{BB}{\displaystyle \frac{N_A}{A}}E\left({\displaystyle \underset{v_{cut}}{\overset{v_{max}^\delta }{}}}{\displaystyle \frac{d\sigma ^\delta (E,v)}{dv}}v𝑑v\right).`$ (2)
Here $`N_A`$ is the Avogadro number; $`A`$ is the atomic weight; indexes $`\gamma `$, $`e^+e^{}`$, $`pn`$ and $`\delta `$ correspond to bremsstrahlung, direct $`e^+e^{}`$-pair production, photonuclear interaction and knock-on electron production, respectively; $`v_{min}`$ and $`v_{max}`$ are the minimum and the maximum kinematically allowed fraction of energy lost for corresponding process; term $`[dE(E)/dx]_{BB}`$ represents Bethe-Bloch formula with correction for density effect. Notice that actually one has to compute CEL, as well as all total cross sections separately for each kind of atoms given material consists of and then add them to each other with different weights but here and below in the text we omit this detail and give only general expressions. One is forced to decompose energy loss because simulation of all interactions with $`vv_{min}`$ would result in infinite computation time due to steep dependence of muon cross sections on $`v`$ (they decrease with $`v`$ at least as $`d\sigma (E,v)/dvv^1`$ and for some processes are not finite at $`v`$ 0). Number of interactions to be simulated per unit of muon path grows, roughly, as $`N_{int}v_{cut}^1`$ along with computation time. On the other hand, setting $`v_{cut}`$ to large value may affect simulation accuracy. Thus, the question is how large value of $`v_{cut}`$ may be chosen to keep result within desirable accuracy?

The second problem can be formulated as follows: if it is necessary to include ionization in SEL, at all? Small energy transfers strongly dominate at knock-on electron production ($`d\sigma ^\delta (E,v)/dvv^2`$), so this process is almost non-stochastic and it seems to be reasonable to omit the last term in expression for CEL (Eq. (2)) and exclude knock-on electrons from simulation procedure when simulating SEL, i.e., to treat ionization completely as “continuous” process which saves computation time noticeable. How much does it affect the result of simulation?

Influence of these factors on simulated result was discussed in literature (see, e.g., Refs. ) but, in our opinion, more detailed analysis on which further improvement of muon propagation MC algorithms could be based is still lacking. In presented work we undertook an attempt to investigate in details with the MUM (MUons+Medium) muon propagation code (Ref. ) the influence of $`v_{cut}`$ and model of ionization energy loss upon result of MC simulation. Resulting errors are compared with uncertainties which come from parametrization accuracy both for sea-level muon spectrum and for muon cross sections. Our analysis allows to rank different uncertainties according to their importance, expose ones to which special attention should be put and choose the most adequate setting of MC parameters at simulation for this or that purpose.

We briefly describe the main features of the MUM code which was a basic instrument for reported investigation and give a short review of principal variables have been studied in Sec. II. Sec. III presents results of simulation for survival probabilities and muon flux intensities with different settings of simulation paremeters. In Sec. IV we analyse results of simulations presented in Sec. III, compare them with other published results and give some basic inferences. Sec. V gives general conclusions.

## II Method

Description of the MUM code which was used to obtain presented results has been published in Ref. . Here we only dwell briefly on those features of the code which seem to be important in the frame of discussed problem:

* the most recent results on parametrizations for muon cross sections are used in the code;
* we tried to decrease the “methodical” part of systematic error due to interpolation, numerical integration, etc., down to as low level as possible, especial attention was put on simulation procedures for free path and fraction of energy lost;
* the most important parameters, like value of $`v_{cut}`$, model for ionization loss, kind of medium, parametrizations for muon cross sections are changeable and represent input parameters for initiation routine;

The accuracy of simulation algorithm used in MUM was shown in Ref.  to be high enough to perform systematically significant analysis which is presented below.

In order to study how different factors influence upon result of simulation we performed several sets of simulations both for propagation of monoenergetic muon beam and muons sampled by real sea-level spectrum (in the later case we limited ourselves by simulation only vertical muons) through pure water down to depths from $`D=`$ 1 km to $`D=`$ 40 km. Of course, depths of more than several kilometers of water for vertical muons do not concern any real detector but simulations for larger depths allow us to study general appropriatenesses for large muon ranges which correspond, for instance, to nearly horizontal muons. Several runs were done for standard rock ($`A=`$ 22, $`Z=`$ 11, $`\rho =`$ 2.65 g cm-3), as well. Muons whose energy decreased down to $`E`$ = 0.16 GeV (the Cherenkov threshold for muon in water) were considered as stopped ones. We tested different settings of parameters which were as follows.

* $`v_{cut}`$, which changed within a range of $`v_{cut}=`$ 10-4 to $`v_{cut}=`$ 0.2. Actually, “inner” accuracy of the MUM code becomes somewhat worse at $`v_{cut}>`$ 5$`\times `$10-2, especially if fluctuations in ionization are not simulated (Ref. ). So, results for $`v_{cut}`$ = 0.1 for $`v_{cut}`$ = 0.2 are presented here only to illustrate some general qualitative appropriatenesses.
* Model for ionization loss which was treated both as completely “continuous” and “stochastic” with simulation of energy lost if knock-on electron energy $`\mathrm{\Delta }Ev_{cut}E`$.
* Parametrization for vertical sea-level conventional atmospheric muon spectrum. Two spectra were tested, namely one proposed in Ref.  (basic):

$$\frac{dN}{dE}=\frac{0.175E^{2.72}}{cm^2ssrGeV}\left(\frac{1}{1+{\displaystyle \frac{E}{103GeV}}}+\frac{0.037}{1+{\displaystyle \frac{E}{810GeV}}}\right),$$
(3)
and widely used Gaisser spectrum (Ref. ):

$$\frac{dN}{dE}=\frac{0.14E^{2.7}}{cm^2ssrGeV}\left(\frac{1}{1+{\displaystyle \frac{E}{104.6GeV}}}+\frac{0.054}{1+{\displaystyle \frac{E}{772.7GeV}}}\right).$$
(4)
* Parametrization for total cross section for absorbtion of a real photon of energy $`\nu =s/2m_N=vE`$ by a nucleon at photonuclear interaction which was treated both according to Bezrukov-Bugaev parametrization proposed in Ref  (basic):

$$\sigma _{\gamma N}=[114.3+1.647\mathrm{ln}^2(0.0213\nu )]\mu b$$
(5)
and ZEUS parametrization (Ref. ):

$$\sigma _{\gamma N}=(\mathrm{\hspace{0.25em}63.5}s^{0.097}+145s^{0.5})\mu b$$
(6)
($`\nu `$ and $`s`$ in GeV and GeV2, correspondingly).
* A factor $`k_\sigma `$ which all muon cross sections were multiplied by to test influence of uncertainties in cross sections parametrization (and, concequently, in energy loss) upon result. We used $`k_\sigma `$ = 1.0 as a basic value but in some cases set also $`k_\sigma `$ = 0.99 and $`k_\sigma `$ = 1.01, which corresponds to decrease and increase of total energy loss by 1%, respectively. Note that it is an “optimistic” evaluation, the real accuracy of existing parametrization for muon cross sections is worse (see Refs. ).

For each run we fixed the muon spectra at final and several interim depths. The differences between obtained spectra were a point of investigation.

## III Results

### A Propagation of monoenergetic muon beams

At the first set of simulations we propagated monoenergetic muon beams of 4 fixed initial energies $`E_s`$ = 1 TeV, 10 TeV, 100 TeV and 10 PeV down to depths $`D`$ = 3.2 km, 12 km, 23 km and 40 km, respectively, through pure water. In each case propagation of 106 muons was simulated. Fig. 1 shows resulting survival probabilities $`p=N_D/N_s`$ (where $`N_s`$ = 106 is initial number of muons and $`N_D`$ is number of muons which have survived after propagation down to depth $`D`$) vs. $`v_{cut}`$ for final and five interim depths. Two curves are given on each plot for two models of ionization. Also results for $`k_\sigma `$ = 1.00 $`\pm `$ 0.01 and for parametrization (6) are presented as simulated with the most accurate value $`v_{cut}`$ = 10-4.

The following appropriatenesses are visible.

* In most cases with the exception of some plots of the lower row and the left column in Fig. 1 (which corresponds to low survival probabilities and low muon initial energies, respectively) uncertainty in our knowlege of muon cross sections gives the principal effect which essentially exceeds ones from other tested parameters.
* The difference between survival probabilities for two models of ionization is the less appreciable the larger muon energy is. It is quite understandable because at muon energies $`E<`$ 1 TeV ionization represents the great bulk of total energy loss, and vice versa, it becomes minor at $`E>`$ 1 TeV. Thus, contribution which is given by ionization at higher energies is small and, the more, its fluctuations do not play an important role. For muons with initial energies $`E`$ 1 TeV fluctuations in ionization become important only at very last part of muon path and “are not in time” to produce some noticeable effect.
* Generally, parametrizations (5) and (6) do not show a noticeable difference in survival probabilities, in most cases it is within statistical error or exceeds it only slightly.
* Increase of $`v_{cut}`$ gives effect of both signs in survival probabilities: function $`p(v_{cut})`$ grows at the beginning of muon path and falls at the last part. The same “both-sign” dependencies are observed for ionization model.
* For $`v_{cut}`$ 0.02 – 0.05 there is almost no dependence of survival probability on $`v_{cut}`$ with the exception of very last part of muon path where survival probability becomes small. Generally, dependence $`p(v_{cut})`$ is the less strong the larger initial muon energy is.

The last item is illustrated complementary by Fig. 2 and Fig. 3 which show that for all initial energies $`E_s`$ simulated survival probability does not depend, in fact, on $`v_{cut}`$ until 90% (for $`E_s`$ = 1 TeV) to 99.5% (for $`E_s`$ = 10 PeV) muons have been stopped.

It was shown above what is result of simulations with different models of ionization and values of $`v_{cut}`$. It was a special point of interest for us to track how does it influence upon behaviour of survival probability.

Fig. 4 shows how muon spectrum resulting from monoenergetic muon beam with initial energy $`E_s`$ = 1 TeV transforms when propagating through pure water down to the depth of 3.2 km. Simulated results for four settings of parameters are presented by four columns of plots. The first three columns represent spectra obtained with ionization included in SEL for $`v_{cut}`$ = 10-4, 10-2 and 0.2. 4th column contains spectra simulated with entirely “continuous” ionization and $`v_{cut}`$ = 10-4. The spectra grouped into the first column represent some “ethalon” because they were simulated with the most accurate tuning both for $`v_{cut}`$ and ionization model. The first three columns demonstrate that “compactness” of spectra at the same depth is the higher the more value of $`v_{cut}`$ is. Put your attention to the right edge of spectra which shifts toward low energies when $`v_{cut}`$ increases (it is the most noticeably for $`v_{cut}`$ = 0.2). The reason is that at any depth energy of the most energetic muons in simulated beam is determined by CEL. These muons due to statistical fluctuations did not undergo interactions with $`vv_{cut}`$ and, consequently, lost energy only by CEL which increases when $`v_{cut}`$ increases. That is why the maximum energy in simulated muon beam is lower for large values of $`v_{cut}`$. Fraction of muons which did not undergo an “catastrophic” act with $`vv_{cut}`$ till given depth grows with increase of $`v_{cut}`$ because free path between two sequential interactions with $`vv_{cut}`$ grows approximately as $`\overline{L}v_{cut}`$. It leads, in particular, to distinctly visible picks in spectra for $`v_{cut}`$ = 0.2 consisted just of muons which lost energy only by CEL. Also, some deficit of low energy muons appears if one sets $`v_{cut}`$ to a large value. In this case left edge of spectrum is provided only with muons which interacted with large fraction of energy lost while for smaller $`v_{cut}`$ an additional fraction of muons comes here. As a result simulated spectrum of initial monoenergetic muons at given depth is more narrow if $`v_{cut}`$ is large and, on the contrary, more wide if $`v_{cut}`$ is small.

Now it is easy to understand how value of $`v_{cut}`$ influences on simulated survival probabilities. When simulated muon beam goes through medium loosing energy both in CEL and SEL processes, its spectrum is constantly shifting to the left (energy decreases). For $`v_{cut}`$ = 10-4 the left part of spectrum reaches $`E`$ = 0 at a smaller depth comparing with larger $`v_{cut}`$ (because in this case spectrum is wider) and survival probability starts to decrease. At the same depth survival probability for $`v_{cut}`$ = 10-2 and $`v_{cut}`$ = 0.2 is still equal to 1. Thus, for the first part of path the survival probability is always larger for large $`v_{cut}`$. At some depth (which is equal to approximately 2.8 km for considered case) compactness of spectra simulated with large $`v_{cut}`$ starts to play an opposite role. Due to this compactness and higher CEL muons stop faster comparing with accurate simulation. So, at the final part of the beam path simulated survival probability for large $`v_{cut}`$ decreases faster comparing with accurate simulation and, for instance, for $`v_{cut}`$ = 0.02 the rest of muon beam which reaches the depth of $`D`$ = 2.72 km (37% of initial number of muons) completely vanishes within the next 30 m of path, while some fraction of muons simulated with $`v_{cut}`$ = 10-4 (0.07%) escapes down to the depth of $`D`$ = 3.2 km.

Qualitatively the same effect leads to the same results if one treats ionization as completely “continuous” energy loss. Again, spectra becomes more narrow since fluctuations in ionization do not work and, as a consequence, survival probability becomes significantly higher comparing with simulation with accurate treatment of ionization at the beginning of muon beam path and falls down essentially faster at the final part of path.

Results presented in Sec. III A show the significant influence which both model of ionization and value of $`v_{cut}`$ have over survival probability for monoenergetic muon beam. But for practical purposes the more important is how this factors do work for real atmospheric muons with a power spectrum?

### B Propagation of muons sampled according to a power sea-level spectrum

In Fig. 5 we present intensity of vertical atmospheric muon flux $`I`$ at different depths of pure water $`D`$ from 1 km to 20 km vs. $`v_{cut}`$ as simulated with muons sampled according to sea-level spectrum (3). Simulation continued until 104 muons reached given depth. Curves for two models of ionization are shown for each depth along with results for $`k_\sigma `$ = 1.00 $`\pm `$ 0.01 at $`v_{cut}`$ = 10-4, parametrization (6) at $`v_{cut}`$ = 10-4, sea-level muon spectrum (4) at $`v_{cut}`$ = 10-4 and all energy loss treated entirely as CEL (for depths $`D`$ 5 km only).

General appropriatenesses for real muon spectrum are qualitatively the same as observed for monoenergetic muon beams.

* For all depths at which neutrino telescopes are located it was found to be better to take into account fluctuations of energy loss simulated by any model than to treat energy loss as completely “continuous”: muon flux intensity computed with non-stochastical model of energy loss is always less comparing with stochastical model, the difference reaches 10 % at 3 km w.e. and 40% at 5 km w.e. (at the depth of 20 km of pure water vertical muon flux computed with ignorance of fluctuations is only 10 % of simulated flux).
* Like in a case for monoenergetic beams 1%-uncertainty in muon cross sections plays the principal role for resulting error in simulated muon depth intensity. This error has a tendency to grow with depth from $`\pm `$2.5% at depth of 1 km w.e. to $`\pm `$15% at 20 km w.e.. But a particular case of this uncertainty, namely difference between parametrizations for total cross section for absorbtion of a real photon by a nucleon at photonuclear interaction from Refs. , does not lead to a significant difference in resulting intensity.
* Difference between muon spectra (3) and (4) leads to uncertainty from -4% ($`D`$ = 1 km) to 16% ($`D`$ = 20 km).
* Error which appears due to simplified, entirely “continuous” ionization lies, commonly, at the level of 2–3 %.
* Dependence of simulated muon flux intensity upon $`v_{cut}`$ is the most weak one comparing with other studied error sources. Function $`I(v_{cut})`$ is almost a constant if $`v_{cut}`$ 0.02–0.05 and changes in a range $`\pm `$1–2% which is very close to statistical error.

We also tried to reveal how value of $`v_{cut}`$ and model for ionization influence upon differential muon depth spectra. No differences were detected which would exceed statistical error. It is illustrated by two figures. Fig. 6 displays simulated mean energies for vertical muon flux at different depths $`D`$ of pure water vs. $`v_{cut}`$ as simulated with muons sampled according to sea-level spectrum (3). Two dependencies are presented at each depth for two model of ionization loss. No appropriatenesses are visible on the plot. In Fig. 7 we present simulated differential muon spectra at four depths $`D`$ = 1 km, 5 km, 10 km, 20 km. Two spectra are displayed on each plot obtained i) for the most accurate value of $`v_{cut}`$ = 10-4 and ii) for the most rough case $`v_{cut}`$ = 0.2. Spectra are normalized to 104 muons and at each depth are divided into three parts along muon energy to keep linear scale at Y-axis. Again, no statistically significant differences are marked.

## IV Discussion

Results obtained in this work are evidence of accuracy in parametrizations for muon cross sections and spectra to be the principal source of uncertainties when simulating muon flux at depths where neutrino telescopes are located. Both for monoenergetic muons and muons sampled according to a power spectrum it gives uncertainty at least of 2–4% to 10-15% in resulting intensity of muon flux, depending upon depth. Unfortunately, this level has at present to be considered as a limit for accuracy of muon propagation codes. Influence of model for ionization exceeds this limit only for monoenergetic muons with sea-level muon energies $`E`$ 10 TeV and only if level of observation is at very last stage of muon path where major fraction of initial muon energy has been lost. Actually, due to steep shape of atmospheric muon power spectrum, an essential part of muons reaches detector location being just on the last part of path. Therefore effect remains noticeable also for real atmospheric muons but in this case uncertainty was found in this work to be much less: 2–3%, which is in an excellent agreement with Refs. , while Ref.  predicts much more significant difference (up to 17%). We suppose this disagreement may result from the fact that “small transfer grouping” technique used for simulation in Ref.  treats muon cross sections to be constant between two interactions in contrast with algorithm used in the MUM code. As was shown in Ref. , in this case switching off the fluctuations in ionization leads to an additinal amplification in effective energy loss for muon energies $`E`$ 1 TeV and, consequently, simulated muon flux intensity must decrease at relatively small depths where muon spectrum is formed by muons with sea-level energies just in a range $`E_s`$ 102–103 GeV. In Ref.  the same simplification was used but reported reults were obtained by simulation with $`v_{cut}`$ = 10-3. At this value of $`v_{cut}`$ role of correct treatment for free path is not significant. Choice of value for $`v_{cut}`$ is of even less importance and again, it is more critical if one investigates monoenergetic muon beam but for power spectrum alteration in $`v_{cut}`$ within $`v_{cut}`$ 0.02–0.05 leads only to 1–2% differences in simulated muon flux intensities which, again, is in a good agreement with level of errors reported in Ref. . Differences between muon flux intensities simulated for different models of ionization and values of $`v_{cut}`$, as obtained in given work and in Refs. , are presented in Fig. 8 and Fig. 9.

Since computation time which is neccessary for simulation depends strongly upon used model of ionization loss and upon $`v_{cut}`$, it seems to be reasonable for most purposes to set $`v_{cut}`$ to $`v_{cut}`$ 0.01–0.05. Anyway, in this case resulting error will be less comparing with error caused by accuracy of parametrization for muon cross sections and sea-level muon spectrum, if one treats free path $`L`$ by an accurate method (Ref. ). Note that “inner” accuracy of the MUM code was found to be worse than “outer” one (which results from uncertainties of muon cross sections) if ionization is set to “continuous” model and $`v_{cut}>`$ 10-2. For fluctuated ionization inner accuracy remains high enough at least till $`v_{cut}`$ = 0.05 (Ref. ). So, using MUM one should set $`v_{cut}`$ to $`v_{cut}`$ 0.01 if ionization is entirely continuous and to $`v_{cut}`$ 0.05 for ionization included in SEL. For test runs it is possible to set $`v_{cut}`$ even to larger values. On the other hand for some methodical purposes it may be useful to simulate fluctuations in knock-on electron production and choose more fine $`v_{cut}`$, e.g., if one wants to exclude an additional error when comparing results of simulations for different models of muon sea-level spectrum with each other.

It is impossible to consider all particular cases and give some strict conformity between setting of parameters at muon MC propagation code and problem to be solved. We had for an object to create a code which would allow to change easy the most important parameters and to perform an analysis which would allow anyone to choose these parameters according to one’s purposes and taste. We hope this object has been achieved with this article.

## V Conclusions

We have presented the detailed investigation for dependence of computed depth muon flux upon the most important parameters of muon propagation MC algorithm, which was done with the MUM (MUons+Medium) code (Ref. ). Contributions of different simplifications and uncertainties to the resulting error have been analysed and ranked. We hope that presented work can be useful for further development of muon transport algorithms which are neccessary for adequate analysis of muon data obtained at underground/water/ice neutrino detectors.

###### Acknowledgements.

We would like to express our gratitude to I. Belolaptikov, A. Butkevich, R. Kokoulin, V. Kudryavzev and V. Naumov for useful discussion and essential remarks which were taken into account in the final version of the text. One of us (I.S.) is also grateful to L. Bezrukov and Ch. Spiering for the attention and support.


text:
# Can dusty Lyman break galaxies produce the submillimeter counts and background? Lessons from lensed Lyman break galaxies

## 1 Submillimeter emission from Lyman break galaxies

Measurements of the cosmic star formation rate density (SFRD) based on surveys for Lyman break galaxies (LBGs) are affected by extinction, and attempts to correct for this effect lead to a substantial upwards revision of the SFRD Since the light absorbed in the ultraviolet (UV) is reradiated in the far-infrared (FIR), LBGs affected by extinction must emit FIR radiation, which will contribute to the submillimeter background.

The expected submillimeter emission from the population of LBGs has recently been estimated based on the observed (if not entirely understood) correlation between the spectral index $`\beta `$ in the UV (defined by the relation $`f_\lambda \lambda ^\beta `$) and the ratio of FIR to UV flux in a sample of local galaxies.  Applying this relation to the LBG population, Adelberger & Steidel  found that the submillimeter counts and integrated background can be accounted for. As noted by these authors, these estimates are uncertain, since the validity of the $`\beta `$-FIR/UV correlation at high redshift has not been established, and the distribution of $`\beta `$ values and its dependence on magnitude, and the luminosity function at faint magnitudes are poorly constrained. It is therefore necessary to verify this analysis by direct submillimeter observations of LBGs.

## 2 Lessons from strongly lensed Lyman break galaxies

The predicted $`850\mu \text{m}`$ fluxes for most LBGs based on their UV properties  are $`1\text{mJy}`$ or less, which is too faint for current instrumentation. Submillimeter observations of individual LBGs are significantly easier if strongly lensed LBGs are targeted. We have used SCUBA on the JCMT to observe two strongly lensed LBGs: the object cB58 at $`z=2.72`$ lensed by the cluster $`\mathrm{MS}\mathrm{\hspace{0.17em}1512}+36`$  and the object G1 at $`z=4.92`$ lensed by the cluster $`\mathrm{MS}\mathrm{\hspace{0.17em}1358}+62`$ Rest-frame UV colors indicate significant reddening in both of these objects.  For both galaxies, the method of predicting the FIR emission based on the UV properties $`\beta `$  implies strong $`850\mu \text{m}`$ emission (Table 1). The results of our SCUBA measurements are given in Table 1. The object cB58 is detected at the $`4.7\sigma `$ level; the object G1 was not detected. In both cases the procedure of predicting the submillimeter flux from the observed color and magnitude in the rest-frame UV  overpredicts the submillimeter emission. The magnitude of the discrepancies is illustrated in Fig. 1. For G1 the discrepancy is not significant, given the scatter in the $`\beta `$-FIR/UV relation, but the factor 14 discrepancy for cB58 is highly significant.

Since this discrepancy is so large, it cannot be attributed to observational uncertainties. However, an element of uncertainty in the analysis is introduced by differential lensing, if the effective amplification factor of cB58 in the UV is a factor 10 higher than that in the FIR. A large discrepancy can only be introduced if most of the UV emission comes from the most strongly amplified portions of the source near the caustic, while most of the FIR emission comes from more weakly lensed regions. This situation would require a very different distribution of FIR and UV emission. However, an extinction correction based on UV properties of one region of a galaxy will not be able to predict the submillimeter emission in a completely different region of the system. Therefore, if differential lensing plays a major role, the physical basis of using the $`\beta `$-FIR/UV correlation disappears. For cB58, a lensing model based on new HST data shows that the amplification factor is at least a factor of 5 at every position, and that the intensity-weighted amplification factor in the UV is about a factor of 25 (on both sides of the fold arc, so that the total amplification is approximately a factor 50). Thus differential lensing can account for a discrepancy of at most a factor of 5, but probably much less, since UV and FIR emission should have a similar morphology for the $`\beta `$-FIR/UV correlation to work.

## 3 Discussion and conclusions

These results demonstrate that attempts to produce the submillimeter counts and background based on an extinction correction applied to the LBG population  are fraught with considerable uncertainty, and fail in the case of the two lensed LBGs discussed here. While our sample is small, the results argue against the validity of the low redshift $`\beta `$-FIR/UV correlation in the case of LBGs. Even if the local correlation were valid for high redshift galaxies as well, it still would not be able to produce the brighter submillimeter galaxies. These objects have luminosities putting them in the class of the ultraluminous infrared galaxies (ULIGs), which do not follow the $`\beta `$-FIR/UV correlation, as shown by recent HST-STIS data of a sample of nearby ULIGs. Since these galaxies already account for $`50\%`$ of the submillimeter background, it is not possible for the LBGs to produce a dominant fraction (let alone all) of the submillimeter background. A more likely situation is that the LBGs account for 25 to 50% of the submillimeter background as indicated by the faint structure in the HDF at $`850\mu \text{m}`$ but that the dominant part of the submillimeter background is made by a small number of ULIGs, which are either not detected in LBG surveys, or if they are detected, cannot be reliably corrected for extinction, in the same way that local ULIGs cannot be extinction corrected based on UV data.

In summary therefore, these results support the view that the brighter $`850\mu \text{m}`$ galaxies making at least 50% of the submillimeter background, form a population which cannot be reliably reproduced by extinction corrections applied to the LBG population.