text:
# Nanolensing of gamma-ray bursts

## 1 Introduction

At present we do not know what constitutes the bulk of the material Universe, i.e. the dark matter. Many dark matter candidates have been proposed, ranging from elementary particles to macroscopic objects such as black holes – see, e.g., Trimble (1987), Ashman (1992), Carr (1994). To date it has been possible to eliminate some suggestions (e.g. brown dwarfs: Tinney 1999), by virtue of a clear conflict between models and data, but a positive identification has not been achieved. One proposal whose implications have not yet been thoroughly explored is that of Hawkins (1993, 1996), who argued, on the basis of photometric monitoring of quasars, that the Universe contains a near-critical density of planetary-mass objects (cf. Press and Gunn 1973). Gravitational lensing by such a population is capable of explaining much of the observed variability (see also Schneider and Weiß 1987; Schneider 1993) — though this should not be taken to imply that intrinsic variations are absent. This suggestion is a radical departure from the now-standard picture of a universe dominated by elementary particles (i.e. Cold Dark Matter – see, e.g., Peebles 1993; Blumenthal et al 1984; Davis et al 1985), and has received little theoretical attention. Building on earlier work (Schneider and Wagoner 1987; Rauch 1991; Seljak and Holz 1999; Metcalf and Silk 1999), Minty, Heavens and Hawkins (2001) described a statistical test of the model based on the light-curves of distant supernovae, but existing data do not allow these ideas to be usefully implemented. A test based on surface brightness variability in low-redshift galaxies has been described by Lewis and Ibata (2001), but this idea also requires data which are not yet available, as does the test based on monitoring quasars seen through low-redshift clusters of galaxies (Walker and Ireland 1995; Tadros, Warren and Hewett 1998). Data from the MACHO and EROS experiments exclude planetary-mass compact objects as a significant contributor to the Galactic dark matter (Alcock et al 1998). However objects which are sufficiently compact that they qualify as strong gravitational lenses at cosmological distances are not necessarily strong gravitational lenses when they are located in the Galactic halo (Walker 1999; see also Draine 1998, Rafikov and Draine 2001). The Galactic microlensing experiments therefore do not directly test the quasar nanolensing hypothesis.

It is obvious that one could attempt to devise further, more sophisticated tests of the nanolensing interpretation based on existing quasar data, but it is a priori unlikely that any such test will yield a definitive result. The main reason for this is source size: the low amplitude of variations seen in quasar optical light-curves means that these sources are necessarily larger than the scale-size of the hypothesised caustic structure, thus smoothing the magnification pattern, and in the process destroying the high spatial-frequency information. (This would be less of a barrier to progress if we had a reliable model for the physical structure of quasars, so that detailed predictions for their lensed appearance could be given with some confidence.) For example the approach of Dalcanton et al (1994; see also Canizares 1982), based on the statistics of equivalent widths of quasar emission lines, is quite insensitive in this regime where the continuum source is resolved by the lenses. To make progress we therefore require a numerous population of small, highly-luminous sources at large distances; gamma-ray bursters constitute such a population. Gravitational lensing of gamma-ray bursts has previously been considered by a number of authors (e.g. Paczyński 1986b; McBreen and Metcalfe 1988; Mao 1992; Gould 1992; Blaes and Webster 1992; Nemiroff and Gould 1995), but their considerations are relevant to lenses which are either much more massive, or much less massive than the planetary-mass range which we consider here. The present work also differs substantively from previous investigations in that the evolution (expansion) of the source structure plays a critical role in our analysis.

Following the launch of the Compton Gamma-Ray Observatory, BATSE (Burst And Transient Source Experiment) soon discovered that gamma-ray bursts are isotropically distributed on the sky, and that they depart from Euclidean source counts at low flux levels (Fishman et al 1994). These discoveries immediately shifted attention away from interpretations based on Galactic neutron stars, which had previously been popular, towards a variety of models in which energies amounting to a significant fraction of $`\mathrm{M}_{}c^2`$ are rapidly released by sources at cosmological distances (see Paczyński 1995). Discovery of X-ray afterglows from these transient events (Costa et al 1997) by BeppoSAX allowed the sources to be accurately positioned, and thus led to the discovery of optical and radio afterglows (van Paradijs et al 1997; Djorgovski et al 1997; Frail et al 1997). Spectroscopy of these optical transients has in some cases revealed absorption lines from gas at redshifts $`z_s1`$, thus firmly establishing the distance scale of the bursters to be cosmological (Metzger et al 1997; van Paradijs, Kouveliotou and Wijers 2000).

Although we still do not know what process injects the burst energy, it is broadly agreed that the radiation observed from GRBs arises from a relativistically expanding source. Relativistic expansion is required in order that we should see gamma-rays at all – at least at energies $``$ MeV – otherwise the inferred source-size is so small, and the photon energy density so high, that two-photon pair-production converts essentially all the gammas into material particles (Cavallo and Rees 1978; Fenimore, Epstein and Ho 1993; Baring and Harding 1997). These considerations require expansion with Lorentz factors of at least $`\mathrm{\Gamma }10^2`$.

At an early stage in the modelling of gamma-ray bursts in a cosmological context, it was pointed out that gamma-ray emission should arise as the the ambient medium undergoes shock compression by the expanding material (Rees & Mészáros 1992). This picture – the “blast-wave”, or “external shock” model – now provides the accepted context for modelling the broad-band (X-ray, optical and radio) afterglows of bursts (Paczyński and Rhoads 1993; Mészáros and Rees 1997; Granot, Piran and Sari 1999), but has fallen into disfavour as an interpretation of the prompt gamma-ray burst. The reason for the demise of this model is simply that it cannot, in itself, accommodate the rapid, large-amplitude variability which is manifest in the gamma-ray data (Fenimore, Madras and Nayashkin 1996; Sari and Piran 1997; see also Dermer and Mitman 1999; Fenimore, Ramirez-Ruiz and Wu 1999). This point has promoted the idea that the gamma-rays arise from internal shocks within the relativistic outflow, so that the temporal variations of the burst reflect the input power variations of the source (e.g. Rees and Mészáros 1994; Kobayashi, Piran and Sari 1997).

Because we are investigating an hypothesis in which apparent variability arises external to the source, as a result of gravitational lensing, arguments which assume that the observed variations are intrinsic immediately lose most, if not all of their force. Given this there are at least three motivations to return to the blast-wave interpretation of the prompt gamma-ray emission: first, this picture allows the initial deposition of energy to be impulsive, and is therefore suitable for a broad range of source models (independent of the specific physics of the energy input) with no fine-tuning required; secondly, the low radiative efficiency of internal shocks (Spada, Panaitescu and Mészáros 2000) appears to be inconsistent with the event energetics in at least some cases (Paczyński 2001); and thirdly, gamma-ray emission from an external shock ought to be present at some level. This emission can be sensibly described by a self-similar, decelerating blast-wave model, similar to those developed for the afterglow emission (e.g. Granot, Piran and Sari 1999), but radiative rather than adiabatic. Such a description leads us to expect a circular source with very strong limb brightening and consequently we anticipate significant flux variations if the limb crosses a caustic. This thin, bright and rapidly expanding ring is a near-perfect instrument for revealing any caustic structure that might be present along the line-of-sight to the source — see also Loeb and Perna (1998), and Mao and Loeb (2001), who considered microlensing of optical afterglows by stellar-mass objects in the low optical-depth limit.

In this paper we use the source model just described to explore the hypothesis that the Universe contains a high density of planetary-mass objects, with an optical depth to gravitational lensing of order unity. Often the characteristic angular scale of the lenses (in arcseconds) is indicated in the nomenclature given to associated phenomena (e.g. microlensing), and following this convention dictates the name “nanolensing” for the effects discussed herein. Although gravitational lensing is the main focus of our attention here, in studying the variations which might be introduced by this process we are also implicitly addressing the physics of the sources themselves. We start by presenting our source model in §2, and in §3 we turn to aspects of the gravitational lensing; light-curves which result from the marriage of these elements are presented and compared with data in §4. The role of parallax is explored in §5, followed by a discussion of related issues concerning both the lenses and the sources (§6), and our conclusions are given in §7.

## 2 Source model

Following Rees and Mészáros (1992), we adopt a spherically-symmetric, ultra-relativistic blast-wave model of the gamma-ray burst phenomenon, in which a total energy of $`10^{52}E_{52}\mathrm{erg}`$ resides in ejecta which expand with an initial Lorentz factor $`\gamma =10^3\gamma _310^3`$ into a homogeneous medium of density $`n\mathrm{cm}^3`$. At first the ejecta coast ($`\gamma `$ constant), but this lasts only for a brief period $`0.1\gamma _3^{8/3}(E_{52}/n)^{1/3}\mathrm{s}`$, and subsequently the ejecta decelerate, with some fraction of the thermalised power appearing as radiation. We will treat these as distinct, self-similar phases of evolution, characterised by the index $`m`$, where

$$\mathrm{\Gamma }^2=\mathrm{\Gamma }_0^2\left(\frac{R_0}{R}\right)^m$$
(1)
at radius $`R`$, and $`\mathrm{\Gamma }=\sqrt{2}\gamma `$ is the Lorentz factor of the shock-wave which precedes the ejecta (Blandford and McKee 1976). In the coasting phase we evidently have $`m=0`$; while $`m=3`$ for adiabatic evolution, and $`m=12`$ for a fully radiative blast-wave – i.e. one in which all of the thermalised energy is promptly radiated away (Cohen, Piran and Sari 1998). The adiabatic approximation is appropriate for the afterglow, when the radiation time-scale is much longer than the expansion time-scale; for the gamma-ray burst itself we adopt the fully radiative solution, $`m=12`$. In this circumstance the gamma-rays arise from a thin shell immediately behind the shock front.

It is convenient to take $`R_0`$ as the radius at which the transition between coasting and semi-radiative solutions takes place, and to model the evolution of the blast wave as if there were an instantaneous transition between these solutions, as the shock crosses this radius. We can estimate the value of $`R_0`$ by noting that the shock will start to decelerate when a significant fraction of the initial energy of the blast has been thermalised (Rees and Mészáros 1992):

$$R_01.2\times 10^{16}\left(\frac{E_{52}}{n\mathrm{\Gamma }_3^2}\right)^{1/3}\mathrm{cm}.$$
(2)
Where numerical estimates are required we adopt the values $`E_{52}=\mathrm{\Gamma }_3=n=1`$ throughout. For some estimates we need to know, in addition, the Hubble constant, which we take to be $`H_0=70\mathrm{km}\mathrm{s}^1\mathrm{Mpc}^1`$.

In this paper we consider only bolometric radiation properties, so that it is not necessary to specify the radiation mechanism.

### 2.1 Kinematics

Because the expansion is spherically-symmetric, the observed source structure is axisymmetric, and at any given time it can be expressed as a function of a single variable, such as the apparent radius, $`r`$, as measured by a distant observer, relative to the line-of-sight to $`R=0`$. Because the expansion is relativistic, it is essential to incorporate light travel-time in any computation of the appearance of the source. In terms of the time, $`t`$, measured by a distant observer at the same redshift (with $`t=0`$ corresponding to the start of the expansion), there is a maximum apparent radius, $`r_{max}`$, from which photons can be received by the observer:

$$r_{max}^{2m+2}=\mathrm{\Gamma }_0^2R_0^m\left[\left(\frac{2m+2}{m+2}\right)ct\right]^{m+2},$$
(3)
and this defines the limb of the source at any time. In general we expect the source to be at non-zero redshift, $`z_s`$, and to allow for this one simply makes the replacement $`tt/(1+z_s)`$, both here and in subsequent formulae. In cases where a numerical estimate is required, we adopt the value $`z_s=5`$ throughout this paper; this corresponds to the median redshift in the GRB source population model of Bromm and Loeb (2002). The radius at $`r=r_{max}`$ is given by $`R_{max}`$:

$$R_{max}^{m+1}=2ct\mathrm{\Gamma }_0^2R_0^m\left(\frac{m+1}{m+2}\right),$$
(4)
and if we define $`qR/R_{max}`$, then for any point on the emitting surface, at any given time, we have

$$\stackrel{~}{r}^2=\frac{q}{m+1}\left[(m+2)q^{m+1}\right],$$
(5)
where $`\stackrel{~}{r}r/r_{max}`$. Notice that there are two values of $`q`$ which correspond to each value of $`\stackrel{~}{r}`$ (other than $`\stackrel{~}{r}=1`$); one of these values lies in the range $`0q<1`$, while the other lies in the range $`1<qQ`$, with $`Q^{m+1}=m+2`$.

Knowing the apparent size of the source as a function of time, we can easily find the apparent expansion speed, $`\beta _{ap}`$ (in units of $`c`$), just by differentiating equation 3:

$$\beta _{ap}^{2m+2}=\mathrm{\Gamma }_0^2\left[\frac{(m+2)R_0}{(2m+2)ct}\right]^m.$$
(6)
### 2.2 Intensity profile

By virtue of being almost coincident with the shock front, the emitting surface exhibits the geometry of a shell moving with Lorentz factor $`\mathrm{\Gamma }`$, as described in §2.1. However, the emitting particles are part of the post-shock flow so their bulk Lorentz factor is $`\gamma =\mathrm{\Gamma }/\sqrt{2}`$, and in consequence the limb of the source ($`\stackrel{~}{r}=1`$) does not correspond to the peak of the observed intensity. This point does not appear to have been recognised previously (cf. Granot and Loeb 2001; Gaudi, Granot and Loeb 2001). For a thin emitting shell, the peak surface brightness corresponds to polar angle $`\psi =\pi /2`$ (as measured in the co-moving frame) relative to the surface normal. For an emitting shell of fractional radial thickness $`\mathrm{\Delta }1`$, the bolometric surface brightness can be described by

$$I(r,t)=\frac{𝒟^4}{\sqrt{\mathrm{cos}^2\psi +\mathrm{\Delta }/4}},$$
(7)
with $``$ being the bolometric intensity emitted normal to the surface, as measured in the rest frame of the emitting shell, and $`𝒟`$ the Doppler factor. (Recall that $`q(\stackrel{~}{r})`$ is double-valued, so that there are implicitly two contributions to the right-hand-side of equation 7.) It is easy to show that rays reaching the observer at a given instant satisfy

$$\mathrm{cos}\psi =\frac{(2m+3)q^{m+1}(m+2)}{(2m+1)q^{m+1}+(m+2)},$$
(8)
and that the ratio of observed-to-emitted photon energies is

$$𝒟=\frac{4(m+1)\gamma _{max}q^{(m+2)/2}}{(2m+1)q^{m+1}+(m+2)},$$
(9)
where $`\gamma _{max}\sqrt{2}:=\mathrm{\Gamma }_0(R_0/R_{max})^{m/2}`$. Now the total energy density in the post-shock gas is just $`4\gamma ^2nm_pc^2`$ (Blandford and McKee 1976), and in the frame of this material the shock itself moves at speed $`c/3`$, so if all of this power is promptly radiated, the emergent normal intensity (i.e. at $`\psi =0`$) is

$$\frac{1}{3\pi }\gamma ^2nm_pc^3.$$
(10)
The above set of equations provides a complete prescription for calculation of the bolometric intensity profile; the fully radiative case, corresponding to $`m=12`$, is shown in figure 1 for a fractional shell thickness of $`\mathrm{\Delta }=4\times 10^3`$. This choice for $`\mathrm{\Delta }`$ is appropriate to, for example, synchrotron emission at $`200`$ keV if the magnetic field energy-density is 10% of the equipartition value. The essential feature of this profile is that it is very strongly limb-brightened, with a sharp peak at $`r0.99295r_{max}`$ (corresponding to $`q^{m+1}:=(m+2)/(2m+3)`$, i.e. $`\psi =\pi /2`$). Source intensity profiles for any index $`3m12`$ exhibit this limb-brightening effect, but as $`m`$ increases the limb-to-centre intensity ratio grows according to $`(32/27)(m+2)^{3m/(m+1)}`$, and the peak-to-centre intensity ratio equals $`(2m+3)^{3m/(m+1)}/8\sqrt{\mathrm{\Delta }}`$; the peak also shifts to larger radii with increasing $`m`$.

The received flux is computed in the usual way, i.e.

$$F=d\mathrm{\Omega }I2\pi r_{max}^2_0^1d\stackrel{~}{r}\stackrel{~}{r}I.$$
(11)
Using equation 5 it is straightforward to rewrite the integral over $`\stackrel{~}{r}`$ as an integral over $`q`$ such that $`0qQ`$. In this way we find that the observed flux should vary as $`Ft^{(2m)/(m+1)}`$, so that $`F`$ grows as $`t^2`$ during the early, coasting phase ($`m=0`$), and then declines as $`t^{10/13}`$ during the self-similar radiative ($`m=12`$) evolution. A simple model light-curve can thus be constructed by applying these two solutions in their appropriate regimes (small/large $`t`$, respectively), and switching between the two at the point where they cross. That is the procedure we shall follow here — our calculations are intended to be illustrative and high precision is not needed. One important point to note is that the luminosity in the self-similar deceleration phase declines too slowly to yield a finite radiated energy at late times, and so the model clearly has a limited range of applicability.

## 3 Lens model

If the Universe contains a mean mass density in lenses $`\mathrm{\Omega }_{lens}`$ (in units of the critical density), then the optical depth to gravitational lensing is $`\tau 0.11\mathrm{\Omega }_{lens}`$ for a source at redshift $`z_s=1`$, and $`\tau 0.62\mathrm{\Omega }_{lens}`$ for a source at redshift $`z_s=5`$ (Turner, Ostriker and Gott 1984). These results are for a universe with a mean total density of $`\mathrm{\Omega }=1`$ (the Cosmological Constant, $`\mathrm{\Lambda }=0`$), and for a given redshift the ratio $`\tau /\mathrm{\Omega }_{lens}`$ increases as $`\mathrm{\Omega }`$ decreases. The study by Schneider (1993) demonstrated that $`\mathrm{\Omega }_{lens}0.5`$ is required (see his figure 6b) if the variability reported by Hawkins (1993) is attributed to gravitational lensing. The nanolensing interpretation of quasar variability is also consistent with the currently favoured combination of cosmological model parameters ($`\mathrm{\Lambda }0.7`$, $`\mathrm{\Omega }0.3`$) provided that essentially all of the matter is in the form of lenses, i.e $`\mathrm{\Omega }_{lens}0.3`$ (Minty 2001). The details of the cosmological model are unimportant here; for simplicity of calculation we have therefore adopted a model universe with $`\mathrm{\Omega }_{lens}=\mathrm{\Omega }=1`$.

The magnification pattern which is introduced by the lenses takes the form of a caustic network, in which the appearance of any source is influenced by a large number of lenses (Paczyński 1986a); an example of such a network is shown in figure 2. In this figure we see a large number of fold caustics – a fold is the lowest order catastrophe, at which two images merge – joined by cusps (catastrophes where three images merge). In the present paper we illustrate the effects of GRB nanolensing through analytic and numerical calculations of a source crossing a single fold caustic (this section). We then (§4) simulate the global lightcurves expected for a source expanding across a caustic network such as that shown in figure 2.

In the immediate vicinity of a fold caustic the mapping between rectangular source coordinates $`(x,y)`$, and image coordinates $`(X,Y)`$, is well approximated by a quadratic form (see Schneider, Ehlers & Falco 1992):

$$x=\frac{1}{2}aX^2,y=\frac{1}{2}Y,$$
(12)
with the $`x`$ axis oriented perpendicular to the caustic, and the curvature, $`a`$, being a positive constant. Here the factor $`1/2`$ in the $`xX`$ mapping is arbitrary (because $`a`$ remains unspecified), but the factor of $`1/2`$ in the $`yY`$ mapping is required in order that the fold be convergence free – a pure shear lens mapping – as all empty-beam lens mappings must be in General Relativity. There are two solutions (images) to these equations for any source located at $`x>0`$, and it is straightforward to determine their magnifications:

$$\mu _\pm =\pm \frac{1}{2\sqrt{2ax}}.$$
(13)
In addition to this image pair, there will typically be a large number of other images present, but the observed flux variations around the time of caustic crossing are completely dominated by the two very bright images which are described by equations 12 and 13. If a burst occurs at any location $`x<0`$, then neither of these two images will be present until the source expands so as to touch the caustic, at which point the observed flux will rise abruptly. On the other hand, if a burst occurs at $`x>0`$ then both images will be present at all times.

There are two separate cases which we must consider, depending on which side of the caustic the burst occurs. For a point source at $`x<0`$ there are $`2n+1`$ images, whereas for $`x>0`$ there are $`2n+3`$; we therefore denote these two regions as the one-image and three-image regions, respectively, with the understanding that there are many more images present, but only the pair created/annihilated at $`X=0`$ is of interest here. If the source takes the form of a thin ring, as is approximately the case for the model we have adopted (§2; see figure 1), then the lensed images appear as shown in figure 3. In our subsequent development we denote the $`x`$ coordinate of the origin of the burst by $`s`$, so that the sign of $`s`$ determines whether one or three images are initially present.

### 3.1 Burst in single-image region ($`s<0`$)

Because the source intensity distribution is axisymmetric, it is useful to consider the properties of an infinitesimally-thin ring, radius $`r`$, of uniform surface-brightness, under the mapping given by equation 12. If the burst occurs at $`s<0`$, then images of the source first appear when the source reaches a radius $`r=|s|`$, and for $`r>|s|`$ the total magnification of the ring (summed over both images) is given by

$$\mu _1=\frac{1}{\pi }\sqrt{\frac{u}{ar}}(\varphi |u)$$
(14)
where $``$ denotes the Incomplete Elliptic Integral of the first kind; $`\mathrm{cos}2\varphi :=s/r`$, and $`u2r/(r+s)`$. Here “magnification” means the total area occupied by the two images of the annulus, divided by the area of the annulus itself. Simple approximations to this exact result are available: $`\mu _1(13ϵ/8)/2\sqrt{a|s|}`$ for $`0<ϵ1`$, where $`ϵr/|s|1`$; and $`\mu _1𝒦(1/2)/\pi \sqrt{ar}`$ for $`r\mathrm{}`$, where $`𝒦`$ is the Complete Elliptic Integral of the first kind, and $`𝒦(1/2)1.8541`$. The exact form of the magnification, given by equation 14, is graphed in figure 4 for the case $`s=1/a`$.

### 3.2 Burst in three-image region ($`s>0`$)

For a burst which occurs in the three-image region ($`s>0`$), the total magnification of an annulus is given by

$$\mu _3=\frac{1}{\pi }\sqrt{\frac{u}{ar}}𝒦(u)$$
(15)
for $`r<s`$. On the other hand for $`r>s`$ we have (as with eq. 14)

$$\mu _3=\frac{1}{\pi }\sqrt{\frac{u}{ar}}(\varphi |u).$$
(16)
Simple approximations can be found for these expressions in the following regimes: $`\mu _31/\sqrt{2as}`$, for $`rs`$; $`\mu _3\mathrm{ln}(32/|ϵ|)/2\pi \sqrt{as}`$, where $`|ϵ|1`$ (and $`ϵ=r/s1`$); and $`\mu _3𝒦(1/2)/\pi \sqrt{ar}`$ for $`r\mathrm{}`$. The exact results are shown in figure 4, for the case where $`s=1/a`$, alongside the corresponding magnification curve for a burst occurring in the single image region.

### 3.3 The median lens

For some purposes (e.g. estimating a characteristic length scale in the source plane) it is important to know whereabouts on the line-of-sight the typical lens is located. This need is most sensibly addressed by determining the redshift at which the optical depth to gravitational lensing is one half of the total optical depth for the source under consideration. Using equation 2.13b of Turner, Ostriker and Gott (1984) this can be quite straightforwardly achieved, and the result is that for a source at redshift $`z_s`$ the median lens is located at redshift $`z`$:

$$1+z=\sqrt{1+z_s}.$$
(17)
Thus for a source at redshift $`z_s=3`$ the median lens redshift is unity, and at low redshift ($`z_s1`$) the median lens is located half-way to the source.

One application of the foregoing is the relationship between transverse dimensions measured in the source plane and in the observer’s plane. Appropriate angular-diameter distances for our circumstance are evidently “empty-beam”, distances (for $`\mathrm{\Omega }=1`$), and these are given in table 1 of Turner, Ostriker and Gott (1984). Denoting observer-lens and lens-source angular-diameter distances by $`D_d`$ and $`D_{ds}`$, respectively, we introduce the “lever arm”, $`\varpi D_{ds}/D_d`$. It is this ratio which, for any given lens, determines the relationship of a transverse distance in the source plane to the corresponding distance in the observer’s plane. It is straightforward to show that the median lever arm (i.e. the lever arm of the median lens) is

$$\varpi =(1+z_s)^{3/4},$$
(18)
and we will make use of this result in §5.

### 3.4 Variability time-scale

In combination with our kinematic model of the source structure (§2.1), the statistical properties of the caustic network determine the predicted variability time-scale of GRBs as a function of source redshift and lens mass. At each occasion when the limb of the source crosses a caustic, a peak is introduced into the light-curve, so the characteristic time-scale between peaks is just the time taken by the source to expand in area by $`1/\mathrm{\Sigma }`$, where $`\mathrm{\Sigma }`$ is the number of caustics per unit area in the source plane. In the low optical depth regime there are four folds in each projected Einstein Ring. At high optical depth, where a caustic network develops, the caustic density is higher by a factor equal to the mean magnification, $`(1\tau )^2`$, and we therefore estimate the caustic density by $`\mathrm{\Sigma }4\tau /(1\tau )^2`$ per unit Einstein Ring area in the source plane. The circumstance $`\tau =1`$ corresponds to $`z_s10.7`$ in our model universe, and the caustic density becomes very large, and variability time-scale correspondingly short, for sources around this redshift.

Using our self-similar ($`m=12`$) source model, we can now estimate the expected variability time-scale, $`t_{var}`$ via

$$t_{var}^1\mathrm{\Sigma }\frac{\mathrm{d}}{\mathrm{d}t}\pi r_{max}^2.$$
(19)
It is worth noting that $`r_{max}^2t^{(m+2)/(m+1)}`$ (see eq. 3), so that the apparent area of the source increases nearly linearly with time, and the variability time-scale is thus not expected to evolve significantly during a burst. The results are shown in figure 5 as a function of lens mass, for two different values of the source redshift ($`z_s=1,5`$), and an initial Lorentz factor $`\mathrm{\Gamma }_0=10^3`$. Considering that bursts may last up to several hundred seconds, and that the available temporal resolution is below 100 ms in existing data, almost the entire temporal range in figure 5 is open to study. We thus recognise that a cosmological population of planetary-mass lenses should introduce variability to gamma-ray burst profiles, suggesting a powerful test of the nanolensing interpretation of quasar variability.

The model does not make a precise prediction for the time-interval between caustic crossings, because the appropriate lens mass is not tightly constrained (nor indeed are the parameters of the source, such as $`\mathrm{\Gamma }_0`$). Hawkins (1993) initially gave an estimate of $`10^3\mathrm{M}_{}`$, but Schneider’s (1993) analysis clearly favours $`10^4\mathrm{M}_{}`$ or even smaller values. Using three-dimensional ray shooting simulations, Minty (2001) finds that these data suggest lens masses of $`10^510^4\mathrm{M}_{}`$. Referring to figure 5, if we adopt $`\mathrm{\Gamma }_0=10^3`$ and lens masses of order $`10^5\mathrm{M}_{}`$, then the variability time-scale is expected to be $`60`$ seconds for a source at redshift $`z_s1`$, and $`4`$ seconds for a source at $`z_s5`$.

## 4 Light curves

We have taken two distinct approaches to the study of lensed light-curves: detailed calculation of the behaviour around the time of a fold-caustic-crossing event, exhibiting the type of profiles which are expected for individual “pulses” (cf. Norris et al 1996); and simulation of light-curves for GRBs seen through a caustic network. For any given burst the latter reflects the random structure of the caustic network in the vicinity of the source. We note here that the propagation times for the various images in our model differ only by amounts $`1`$ ns, and this is negligible for our purposes.

### 4.1 Individual “pulses”

Making use of the source intensity profile computed in §2, and the magnification curves derived in the previous section, it is now a straightforward exercise to compute the flux history, $`F(t)`$, from

$$F(t)dr\mathrm{\hspace{0.17em}2}\pi rI(r,t)\mu ,$$
(20)
where the total magnification, $`\mu `$, can be written as the sum of the “unlensed” images (actually, a large number of lensed images whose total flux is roughly constant over the time-scale of the caustic crossing), plus the bright pair of images associated with the fold caustic under consideration: $`\mu =1+\mu _{1,3}`$. (The choice of unity for the constant term in this relation is somewhat arbitrary.) As with equation 11, this formula can easily be rewritten as an integral over the variable $`q`$. After substituting for $`I(r,t)`$ (equations 7–10) and $`\mu _{1,3}`$ (equations 14–16) it is straightforward to evaluate equation 20 numerically; to do this we used the Mathematica software package. The results are shown in figure 6, for bursts occurring in both one- and three-image regions. For these calculations we have assumed $`|s|=1/a`$, and we have taken the fractional radial thickness of the emitting shell to be $`\mathrm{\Delta }=4\times 10^3`$. For each light-curve, the limb of the source first touches the caustic at $`t=1.5`$ seconds (an arbitrary choice).

Some aspects of figure 6 require immediate comment. The difference in overall normalisation of the two curves at early times simply reflects the ratio of total magnifications at the onset of the burst – see figure 4 – roughly $`(1+1/\sqrt{2})/(1+0)1.7`$. This property is a direct result of the simple lens model we have adopted (i.e. the quadratic relationship between source and image coordinates), and it should be remembered that this approximation is accurate only in the immediate vicinity of the caustic. Consequently it is the behaviour seen in the light-curves around $`t=1.5`$ s which is of prime interest, rather than the global properties. The essential features to note are, therefore, that: (i) the presence of a fold caustic introduces a significant peak into the light-curve, close to the time at which the limb of the source first touches the caustic, and (ii) this peak can be either cuspy, and approximately time-symmetric, or rounded and asymmetric. These features could, in fact, have been anticipated simply by examining the magnification profiles shown in figure 4, bearing in mind that the source structure is basically just a thin ring (figure 1).

### 4.2 Simulated light-curves

Although the circumstance we are considering is one in which a population of gravitational lenses is distributed along the entire line-of-sight, we have made use of simulations which are two-dimensional, with all of the lenses located in a single plane. This is partly because three-dimensional ray-tracing simulations are much slower than their two-dimensional counterparts, but mainly because a well-tested two-dimensional code was available to us (Wambsganss, Paczyński & Katz 1990). Two extensions of this code were required for our purposes: convolution of the magnification map with our model source intensity profile (§2); and a progressive “zoom” which recomputes the magnification map on larger scales as the source expands. The latter feature was implemented with an incremental “zoom” factor of 1.2, and our computational grid was 2048 pixels on a side, so that our source radius was in the range 800 to 1000 pixels. Despite this our resolution is only just adequate for the task, as the reader can easily verify by examining the radial width of the source intensity distribution (figure 1). An estimate of the Full-Width at Half-Maximum (FWHM) of the source intensity peak can be obtained from

$$\mathrm{FWHM}\left(\frac{m+2}{2m+3}\right)^{(m+2)/(m+1)}\frac{2\sqrt{3\mathrm{\Delta }}}{2m+3},$$
(21)
and for $`m=12`$, $`\mathrm{\Delta }=4\times 10^3`$ this yields (coincidentally) $`\mathrm{FWHM}4\times 10^3`$. With a source radius of 800 pixels the intensity profile is thus 1.5 times over-sampled.

To illustrate the type of lightcurves which might be expected, we have undertaken simulations for sources at redshift $`z_s=1`$ and $`z_s=5`$, corresponding to optical depths to gravitational nanolensing of $`\tau =0.11`$ and $`\tau =0.62`$, respectively. The results of these simulations are shown in figures 7 and 8, for lenses of mass $`10^8,10^5\mathrm{M}_{}`$, respectively. Because the source evolves in a self-similar (scale-free) manner, it is not necessary to recompute magnification curves when the adopted lens mass is changed: different lens masses simply require a rescaling of the time axis in panel (a) of these figures, thus changing the characteristic time-scale of any variability, but not its amplitude. The choice of lens mass was made, for each of the two redshifts, so as to maximise the variability within each of the synthetic bursts.

All of our simulations follow the same basic pattern: starting from an initial magnification which is below unity, the mean magnification (i.e. average magnification over the whole source) exhibits rapid variability, the amplitude of which decreases as the source expands. These features are readily understood. Initially the source is likely to be demagnified, because most of the area of the source plane is occupied by regions with magnification below unity. Subsequently, as the source expands across the caustics, variability is seen. However, accompanying the increase in source radius, is an increase in the thickness of the high intensity ring and the modulation of the mean magnification of the source is consequently reduced. Simultaneously the time-averaged magnification trends towards its large-scale average value as the source covers an increasing area.

There are two important points which are evident on examining figures 7 and 8. First, the time-scale of the variations in the simulated light-curves is in broad agreement with the rough estimates given in §3.4 (see figure 5). Secondly, the depth of modulation seen in our simulated magnification curves is quite small, being typically of order 10%, and in the light-curves themselves, where the secular evolution is strong, it is not always apparent that the flux is varying on short time-scales.

### 4.3 Comparison with observations

Many high-quality gamma-ray burst light-curves were acquired by the BATSE experiment on the Compton Gamma-Ray Observatory (e.g. Fishman et al 1994), and the substructure within these light-curves has been analysed by Norris et al (1996: N96), in the case of long, bright bursts. The principal conclusion of this analysis was that these bursts can be regarded as being made up of a number of “pulses” having a typical separation of order one second (but which are not periodic). N96 found that the typical pulse shape is cuspy, similar to the peak near $`t1.5`$ s in the top curve in figure 6.

Comparing these findings to the results of §§3,4, we can identify three key points. First, the observed time-scale of order one second between “pulses” is comparable to the predicted time-scale (§§3.4,4), but only if the lenses have low mass ($`10^8\mathrm{M}_{}`$), or the sources are at high redshift. Secondly, the typical “pulse” shape observed by N96 is similar to the $`31`$, i.e. image-annihilating, caustic crossing discussed in §3.2, and shown in figures 4 and 6. However, the theory presented earlier in this paper gave no reason to suppose that $`31`$ crossings should be preferred over $`13`$ crossings, and in fact one does expect both to be present. Both are present in the simulated light-curves shown in figures 7 and 8. One could appeal to an observational bias which favours $`31`$ caustic crossings over their $`13`$ counterparts, because the high and cuspy peaks of the former render them more visible in the light-curves (see, particularly, figure 4). However, this does not seem to us to be a strong argument. Finally, although the single fold-crossing events modelled in §3 are obvious in the light-curves, it is plain that the simulations in §4 do not yield sufficient modulation to be able to explain the profound variations seen in the data.

On the last of the above points we note that there is no discrepancy between the calculations undertaken in §3 and the simulations of §4. The difference between the computed modulation depth in the two cases simply reflects the fact that the caustics present in the simulations exhibit more curvature in the lens mapping than was assumed in §3.

## 5 Parallax

A key aspect of the lensing model we are investigating here is that of parallax. The importance of parallax phenomena depends on the spatial separation of the detectors which are employed, in comparison with the spatial scale on which the magnification pattern changes. In the low optical depth regime, the latter scale is given by

$$\left[4\frac{GM}{c^2}\frac{D_sD_d}{D_{ds}}\right]^{1/2}60\sqrt{M_4}\mathrm{AU},$$
(22)
where $`D_s`$ is the angular-diameter distance of the source from the observer, and to arrive at the numerical result we have assumed $`z_s=3`$ and considered the median lens. Hence for detectors separated by a few Astronomical Units, parallax might be expected to show up only for lenses of mass $`M10^7\mathrm{M}_{}`$ (Nemiroff and Gould 1995). However, in the model we are considering the expanding source crosses caustics, and in the vicinity of caustics the magnification pattern varies significantly on spatial scales much smaller than given by equation 22. Indeed at caustic crossings the transverse flux gradient is limited only by the source structure, and parallax may be observable even over very modest baselines. Specifically, if the width of the high intensity ring of the source is taken to be $`4\times 10^3`$ of the apparent radius (see §4.2), and the latter is approximately 6 AU at $`t=10`$ seconds ($`E_{52}=n=\mathrm{\Gamma }_3=1`$), then parallax may be evident over transverse length scales as small as $`4\times 10^{11}`$ cm. In this circumstance the main instrumental requirement is for good temporal resolution in the detectors, and high signal-to-noise ratio, so that the times at which caustic crossings occur can be precisely determined (Grieger, Kayser and Refsdal 1986; Hardy and Walker 1996). Gamma-ray bursts are typically recorded with temporal resolution of tens of milliseconds, or better, often at high signal-to-noise, so that the available timing precision is generally rather good. Moreover, we wish to emphasise that, provided there are significant increases in flux associated with a caustic crossing, the ability to precisely time such an event is essentially independent of lens mass and this parameter does not influence the detectability of the two parallax phenomena we shall describe.

### 5.1 Displaced burst locations

Suppose we observe a burst with a pair of identical detectors whose spatial separation is small in comparison with the scale on which the lens magnification pattern changes. The caustic network will appear very similar as seen from the two different locations, but with a slight shift on the sky in a direction parallel to the sky-projection of the vector separation between the detectors, $`\stackrel{}{b}_{12}=\stackrel{}{x}_1\stackrel{}{x}_2`$, and by an amount proportional to $`b_{12}`$. Consequently the two light-curves are almost identical in structure, but the caustic-crossings are seen to occur at slightly different times at the two detectors. Thus if burst time-of-arrival is used to determine which direction the wavefront has arrived from, then the present model predicts that the inferred location of the source is shifted from the true location (cf. McBreen and Metcalfe 1988).

In the absence of any lensing phenomena, under the assumption of a source at infinity (plane wavefront), a measured difference in arrival time, $`\mathrm{\Delta }t_{12}=t_2t_1`$, between detectors tells us that the burst came from a direction defined by the unit vector $`\widehat{n}`$, such that $`\widehat{n}\stackrel{}{b}_{12}=c\mathrm{\Delta }t_{12}`$. Thus the source is constrained to lie on a circle of angular radius $`\alpha :=\mathrm{cos}^1(c\mathrm{\Delta }t_{12}/b_{12})`$, centred on the direction $`\widehat{b}_{12}`$. If a timing offset, $`\delta t_{12}`$, is introduced then the corresponding angular offset, $`\delta \alpha `$, is given by $`\mathrm{sin}\alpha \delta \alpha =c\delta t_{12}/b_{12}`$. The timing difference which is expected as a result of parallax, in our model, is roughly $`\delta t_{12}\varpi \mathrm{sin}\alpha (b_{12}/\beta _{ap}c)`$, for any given caustic crossing. (The precise value of $`\delta t_{12}`$ depends, of course, on the orientation of $`\stackrel{}{b}_{12}`$ relative to the orientation of the caustic.) Thus, for a burst which exhibits only a single caustic crossing event, the expected angular offset is

$$|\delta \alpha |\frac{\varpi }{\beta _{ap}}.$$
(23)
Some caution should be exercised in applying this result, because in our case, where a network of caustics is present, the folds cannot be readily associated with individual lenses, and because the lenses are distributed along the whole line of sight it is not clear how the effective value of $`\varpi `$ should be estimated. We shall nevertheless employ equation 23 as if the values of $`\varpi `$ were in one-to-one correspondence with lenses, and use the median value of $`\varpi `$ as given in eq. 18. For a source at redshift $`z_s=3`$, at an observed time $`t10`$ s after the start of the burst, equation 6 leads us to expect an apparent expansion speed of $`\beta _{ap}320`$, and projecting this into the observer’s plane yields the result $`|\delta \alpha |4`$ arcmin. This estimate is valid for a burst which exhibits only a single caustic crossing; if the burst has temporal substructure with $`N_p`$ peaks in the light-curve, then we expect $`|\delta \alpha |`$ to be smaller by a factor of order $`\sqrt{N_p}`$. We now turn to the question of whether such systematic errors are either evident in existing data, or are significantly constrained by them.

Comparing burst arrival times between three spacecraft localises any source to one of two possible intersections between three separate loci (one locus derived from each spacecraft pair). Because the observable quantity is the burst arrival time, and each locus reflects the difference between a pair of arrival times, only two of these loci provide independent information on the burst location; in other words the location is not determined with any redundancy if only three spacecraft are employed. If four or more interplanetary spacecraft are available for burst triangulation, then redundant positioning is possible, and systematic errors can therefore be revealed using only the triangulation data. Many bursts were redundantly positioned using the first generation of interplanetary gamma-ray burst detectors (Atteia et al 1987), and in two of these cases there were highly significant discrepancies (ten standard deviations) between the localisations. However, in at least one of these two cases (GRB 790329) the cause of the discrepancy appears to be understood (Laros et al 1985), and we adopt the conservative assumption that the other discrepant burst (GRB 790116) is also attributable to effects other than the parallax phenomenon under discussion here. The typical level of agreement between redundantly determined burst locations excludes errors at the level $`|\delta \alpha |1^{}`$ (Laros et al 1985), but this result does not strongly constrain the model we have presented because offsets as large as this are expected only in exceptional circumstances. Specifically, offsets larger than a degree can only occur if the apparent expansion speed is $`\beta _{ap}<60`$ and the source is at low redshift ($`z_s<1`$) and a caustic crossing event occurs.

There are no four-spacecraft triangulations available from the modern (BATSE and post-BATSE) era of gamma-ray burst studies, and we have therefore searched the literature for bursts which have been accurately located both by triangulation and by an independent method. The best sample of this type that we found is comprised of nine gamma-ray bursts which were observed by BATSE, Ulysses and either the Pioneer Venus Orbiter (Laros et al 1998) or the Mars Observer (Laros et al 1997) – for which triangulated positions are therefore available – all of which were independently positioned by the rotation modulation collimator of the WATCH experiment (Sazonov et al 1998). The smallest error associated with the triangulated positions ($`1\sigma `$ statistical plus systematic) for these bursts is 20 arcseconds, while the median error is close to 3 arcminutes. Moreover the WATCH error circles all have radii exceeding 14 arcminutes ($`1\sigma `$ statistical plus systematic). Considering that we expect deviations of order an arcminute the WATCH error circle, in particular, is so large that at first sight it seems hopeless to attempt to constrain the model via astrometry. The situation is not quite as bad as it first appears because pairs of Inter-Planetary Network (IPN) loci sometimes intersect at a very acute angle, in which case the location of the intersection point is very sensitive to any errors in $`\alpha `$, thus allowing us to gauge the magnitude of the IPN errors. There are four examples of this type of configuration amongst our sample of bursts. However, none of the implied errors is highly significant, when compared to the estimated measurement error, nor do these results provide any powerful constraints on the theoretical model. We have therefore relegated the details of our astrometric analysis of this sample to an appendix, and in the present section we confine ourselves to a brief discussion of the outcome.

For the four bursts for which we could gauge the actual error in $`\alpha `$, the mean error value is only $`1.2`$ times the estimated error ($`1\sigma `$ statistical-plus-systematic), and is thus consistent with expectations. It is therefore appropriate to quote the results of our analysis in the form of upper limits. To do this we take an upper limit of $`3\sigma `$ in each of the four cases; in order to compare with our theoretical prediction (equation 23) we then convert to an upper limit on the typical parallax error on each caustic crossing, by multiplying the $`3\sigma _\alpha `$ limit by $`\sqrt{N_p}`$. The values of $`N_p`$ for BATSE451 and BATSE1698 were taken from Norris et al (1996): $`N_p=4`$ and $`N_p=7`$, respectively. For BATSE2387 and BATSE907, examination of the archival light-curves (http://cossc.gsfc.nasa.gov) reveals that $`N_p=1`$, and $`N_p5`$, respectively. The resulting limits on parallax for a single caustic crossing are thus deduced to be $`|\delta \alpha |<2.7,\mathrm{\hspace{0.17em}37},\mathrm{\hspace{0.17em}2.6},\mathrm{\hspace{0.17em}37}`$ arcmin, from BATSE451, 907, 1698 and 2387, respectively. Even the lowest of these limits is only comparable to the predicted value ($`4`$ arcmin for a caustic crossing occurring at $`t=10`$ s, and a source at redshift $`z_s=3`$), and so we conclude that existing gamma-ray burst astrometry does not significantly constrain the model we have presented.

### 5.2 Burst profiles

For bursts lasting a few seconds or more it is usually the case that the light-curves contain several peaks. In a nanolensing model these would be identified with caustic crossings, with each peak marking the time at which the limb of the expanding source first touches the caustic. Now the relative timing offsets introduced by parallax are, of course, different for different caustics, so we expect that parallax should alter the time-interval between the peaks of a given burst, as measured by physically separated (but otherwise identical) detectors. The effect is evident in figure 9, which shows two simulations of an expanding source, as described in §4.2, having slightly different initial burst locations. (In respect of the actual calculations undertaken this is equivalent to displacing the observer, and thus simulates the effect of parallax.) One of the caustic crossings visible in figure 9 is seen at almost the same time ($`t0.013`$) in both light-curves, whereas the other is seen at different times in the two simulations. The expected timing difference is simply the quantity $`\delta t_{12}`$, estimated in §5.1, and this evaluates to $`|\delta t_{12}|1`$ s, for $`z_s3`$, $`b_{12}\mathrm{sin}\alpha 3`$ AU and $`t10`$ s. This value is much larger than the limiting time-resolution of the Ulysses and BATSE instruments (Hurley et al 1992; Fishman et al 1994) and should therefore be detectable if the signal-to-noise ratio of the data is sufficiently high. Such measurements offer a very powerful test of the model we are proposing, both because the anticipated size of the effect should render it measurable, and because no other model mimics this behaviour. In particular, the test is more powerful than one based on apparent burst locations (§5.1), because there is no requirement for absolute timing information; only relative timing is important and essentially all spacecraft yield relative timing information with high accuracy.

To emphasise the clarity of this test let us review what should be seen in the absence of any parallax effects. Two detectors observing in the the same energy band should record source flux variations which, once the light-curves have been shifted so as to align the times of burst onset, are in direct proportion to each other. Detector efficiency and spacecraft geometry ordinarily do not change significantly on time-scales $`100`$ s, so a single constant of proportionality should hold throughout each burst. In the case of spacecraft in low Earth orbit, such as the Compton Gamma-Ray Observatory, some fraction of the sky is occulted by the Earth, and it is possible for ingress/egress to occur during a burst. But this is a rare phenomenon; moreover it is easily recognised, and therefore not a concern in the present context. Particle backgrounds, and their trends with time, inevitably differ between any two detectors, but in normal circumstances they also change relatively little on the time-scale of a gamma-ray burst and a linear trend in the background over the duration of the burst is a suitable model for most purposes. Electron precipitation events constitute a variable particle background in low Earth-orbit which in some respects can mimic a gamma-ray burst (Horack et al 1992); however, these events are readily distinguished from gamma-ray bursts by comparing count rates amongst the BATSE detectors. Finally, detector dead-time introduces non-linearity (different for different detectors) which can significantly affect the recorded light-curve, causing the count-rate to saturate if the burst is bright. Saturation will not affect the peak location, because the recorded count-rate is a monotonic function of the flux, but will displace the burst centroid; however, saturation can be recognised trivially, from the count rate, and affected bursts can be excluded from consideration. In summary, none of the effects described above should be confused with the parallax signal we have discussed.

In practice, testing for parallax effects will not be done with a pair of matched detectors, but with whatever instrumentation is available in the IPN, so we will be comparing data from detectors with different energy responses. This presents us with a potential problem in that many of the “pulses” which are seen in GRB light-curves exhibit spectral evolution (N96), and therefore have different shapes in different energy bands. Consequently, with imperfect knowledge of the (evolving) burst spectrum or the energy response of either of the detectors, the pulse will appear to be slightly shifted in time, and this effect can mimic a parallactic timing offset. Fortunately, not all pulses are of this type (N96); a good fraction of pulses are narrow, time-symmetric events, whose peak and centroid locations are insensitive to photon energy (N96). This subset of pulses allows parallax searches to be attempted even with detectors which are not closely matched in their energy response.

A large number of GRB light-curves are freely available at http://cossc.gsfc.nasa.gov/batse (the BATSE archive) with 64 ms resolution in each of four energy bands. For comparison with these data we have searched the literature for published GRB light-curves from the Ulysses spacecraft (Hurley et al 1992). This comparison offers very large baselines (up to 6.3 AU), for which the timing offsets should be correspondingly large. Although Ulysses has detected a large number of GRBs, including hundreds which were also detected by BATSE (Hurley et al 1999a,b), we were able to find only 3 published Ulysses light-curves for GRBs in the BATSE catalogue. Two of these – BATSE2151 (GRB930131, the “Superbowl Burst”, Hurley et al 1994), and BATSE143 (GRB910503; Hurley 1992) – have very modest projected separations between the spacecraft (1.4 AU and 1.6 AU, respectively). We have therefore focussed our attention on BATSE1663 (GRB920622), for which the Ulysses-BATSE projected separation (i.e. transverse to the line-of-sight to the burst) was 3.8 AU. This burst has been analysed in detail by Greiner et al (1995).

We digitised the Ulysses light-curve directly from figure 1 of Greiner et al (1995) and this light-curve is reproduced in figure 10, together with the BATSE light-curve for the same interval. The BATSE data we plot here are the sum of the two lowest energy channels only, corresponding to the photon energy range 25–100 keV; this range is similar to the sensitive range of the Ulysses detector (15–150 keV; Hurley et al 1992). The BATSE data have had a linear background subtracted from them; the background model was determined by taking the median values of 1 minute intervals starting 3 minutes before and 4 minutes after burst onset. (This procedure could not be used for the Ulysses data because the published light-curve does not span a large enough time interval.) The two light-curves were then aligned in time such that their cross-correlation is maximised, and an estimate of the Ulysses detector background was made by maximising the cross-correlation as a function of assumed background, at constant temporal offset. (We note that varying the assumed background yielded insignificant changes in the temporal alignment.) The estimated Ulysses background was in this way found to be 107 counts/bin; the peak value of the normalised cross-correlation between light-curves is 0.974.

Although the temporal alignment between the two spacecraft was determined to a precision of 10 ms, the accuracy of the procedure is poorer than this. We estimated the uncertainties associated with our method by cross-correlating the Ulysses data with the two highest energy BATSE channels. These data are for the same burst, but represent photons of energies higher than those to which Ulysses is sensitive. This procedure yielded an alignment which differed by 70 ms, and a background rate which differed by 2 counts/bin, when referred to the results of the cross-correlation with the lowest energy BATSE data. Our method can be reasonably expected to be in error by less than the shifts exhibited here.

The flux scale in figure 10 is chosen to be that of the Ulysses data, and the BATSE count rates have been rescaled accordingly (and we have then added a constant background of 107 counts/bin to the scaled BATSE data). This choice is appropriate because the BATSE data are of much higher signal-to-noise ratio than the Ulysses data, and the Ulysses flux scale is thus germane to the statistical significance of any comparison between the two. We have also rebinned the BATSE data from the original value of 64 ms to the 250 ms of the Ulysses data, and for clarity we have offset the BATSE data by $`100`$ counts/bin. It can be seen from figure 10 that the BATSE and Ulysses light-curves are very similar, as expected, although there some evident differences in the $`t=1619`$ s region, and the Ulysses light-curve seems to manifest more fluctuations in the region $`t>20`$ s than would be expected on the basis of the BATSE light-curve. Quantitatively we can say that the BATSE data form an acceptable model of the Ulysses data, with a $`\chi ^2`$ statistic of 76 for the difference between the two, over the first 19 seconds of the burst, with $`(19/0.25)3=73`$ degrees of freedom.

However, the key question for us here is not the overall similarity of the two light-curves, but whether or not the temporal substructure occurs at different times at the two spacecraft. In figure 11 we concentrate our attention on the central region of the burst, where the BATSE light-curve manifests obvious substructure. In particular the regions $`t=7.512`$ s and $`t=14.519`$ s show a number of peaks which one can use to test the model. We can compute the combined $`\chi ^2`$ of these intervals: it is $`\chi ^2=48.5`$, with 36 degrees of freedom. Clearly the BATSE-derived model does not perform as well in these intervals with temporal substructure as it does elsewhere; however, the model is still acceptable, with this value of $`\chi ^2`$ being realised by chance in 8% of trials. Moreover, a simple $`\chi ^2`$ test does not reveal whether the model is performing poorly because the temporal substructure is shifted in time, or for some other reason such as the different energy response of the Ulysses detector compared to the combination of BATSE channels 1 and 2.

Examination of figure 11 suggests that both of these effects may play a role. In particular we note that the feature in the low-energy BATSE light-curve just after $`t=18`$ s is clearly chromatic, with a centroid which occurs earlier at higher energies, and we must therefore allow that in the Ulysses data this feature may exhibit a profile different from the BATSE model. By contrast, at $`t=15`$ s, the BATSE model shows a well-defined event, for which we can tentatively identify a counterpart peak in the Ulysses light-curve, but this counterpart occurs 0.4 s later. Although there is some difference in the profiles of this event between high- and low-energy BATSE light-curves, it does not seem reasonable to attribute the BATSE-Ulysses shift to chromatic effects because the differences are modest whereas the shift is comparable to the width of the feature. Another case where the data suggest a temporal offset occurs around $`t=11`$ s. Here we find two narrow, approximately time-symmetric peaks in the BATSE model. In the Ulysses data there is no question about which peaks to identify as their counterparts – the two unresolved peaks either side of $`t=11`$ s – but the second peak appears earlier in the Ulysses data than in the BATSE model. It does not seem reasonable to attribute this to chromatic effects because the profiles of the low- and high-energy BATSE data are quite similar. However, we still cannot be confident of having found an example of the effect we are looking for, because (i) the peak is so narrow that it is undersampled in the Ulysses light-curve, and (ii) the modest signal-to-noise ratio of these peaks in the Ulysses data mean that any comparison cannot yield highly statistically significant results.

In summary, we have compared published Ulysses data with archival BATSE data for a single burst having a large projected separation between spacecraft. This comparison yielded some hints of parallax, but the evidence is of low statistical significance and therefore not compelling. By reprocessing the raw Ulysses data – which have a temporal resolution eight times finer than the published light-curve – it should be possible to make a more meaningful comparison between the two datasets. In particular the higher temporal resolution would be valuable in studying the peaks close to $`t=11`$ s, which are unresolved in the published Ulysses light-curve. A further improvement on our analysis would be to derive a model light-curve using the point-by-point spectral data (“colours”) derived from BATSE, together with the spectral response matrices of the two spacecraft.

## 6 Discussion

Although the burst (BATSE1663) studied in §5.2 was chosen for its relatively large projected separation between the spacecraft, it is not extreme in this respect, and there are $`80`$ instances (Hurley et al 1999a,b) of larger projected separations amongst the bursts detected by both BATSE and Ulysses. Indeed the principal criterion for using the BATSE1663 data was that the data are published. This burst merited detailed study, and hence publication of the light-curve, because it happened to occur within the COMPTEL field-of-view (Greiner et al 1995). Amongst the many other bursts in which one might consider searching for parallax, how is one to choose? The expected temporal offset is large (see §5.1) if (i) the projected baseline $`b_{12}\mathrm{sin}\alpha `$ between spacecraft is large, and (ii) the source redshift is low (hence $`\varpi 1`$), and (iii) the apparent expansion speed of the source ($`\beta _{ap}`$) is low. Unfortunately only the first of these criteria can be securely determined directly from the data at hand. However, if all other parameters are held fixed, then we expect that high expansion speed and high source redshift would both lead to large numbers of caustic crossing events, so it makes sense to select bursts which exhibit a small number of peaks in their light-curves. Ideally the overlap amongst these peaks should be small, so that the properties of each peak can be characterised more-or-less independently of the others. These criteria also help to avoid the potential problem of peak confusion: if parallax displaces a burst peak by an amount comparable to the separation between peaks in one of the light-curves, then it becomes difficult to identify counterparts between the light-curves from the different spacecraft. Finally, and most obviously, a significant detection absolutely requires high signal-to-noise ratio, so bright bursts are strongly preferred.

If we could increase the baseline between the two detectors beyond the few AU which characterises interplanetary spacecraft, our model predicts that the amplitude and location of the substructure peaks would differ more and more as the baseline increased, and it would eventually become difficult to identify counterpart peaks between the two profiles. Ultimately, if the separation were to reach values $`60\sqrt{M_4}`$ AU, the entire caustic pattern would differ as seen from the two locations, and there would consequently be little resemblance between the recorded burst profiles. Although there seems to be no immediate prospect of separating a pair of detectors by such a large distance, it is nevertheless possible to reach this regime, and thus to observe completely different profiles, if a gamma-ray burst is “macroscopically” gravitationally lensed. More specifically, if a gravitational lens forms multiple images which are split by an angle very much greater than $`15\sqrt{M_4}`$ nano-arcseconds, then the caustic crossings seen in each image will bear little resemblance to one another. This, of course, would be true even for micro(arcsecond)lensing due to stars at cosmological distances. Thus in the case of lensing by galaxies, and aggregates of larger mass, any lensed “echo” should exhibit a different temporal profile to that of the counterpart signal which arrives first (cf. Williams and Wijers 1997; Paczyński 1986b; Blaes and Webster 1992). It is expected that gravitational lensing by galaxies should yield echoes of roughly 0.1% of bursts (Paczyński 1986b; Mao 1992), but this effect has never been clearly detected. The null results to date (Nemiroff et al 1994; Marani et al 1999) are not highly statistically significant, nor is it trivial to recognise echoes even if they are simply scaled copies of the “original” – see Wambsganss (1993), and Nowak and Grossman (1994). This means that the lack of any identified echoes, while consistent with the present model, is not a decisive argument in its favour. On the other hand, if examples were found with complex temporal substructure faithfully reproduced in an echo, then it follows that the observed substructure is not due to gravitational nanolensing.

It has previously been noted that if the Universe is populated with low mass lenses then interference phenomena may manifest themselves in gamma-ray burst data (Gould 1992; Stanek, Paczyński and Goodman 1993; see also Deguchi and Watson 1986). The lens masses for which this effect is relevant are usually thought to be many orders of magnitude smaller than the planetary-mass bodies considered here. In the case of lensing by a fold, where the angular splitting of the images vanishes for a source located precisely on the caustic, interference fringes can be realised for larger masses than in the case of isolated gravitational lenses. However, strong interference fringes also require that the source be at most comparable in size to the Fresnel scale – $`\sqrt{D_s\lambda }`$ for wavelength $`\lambda `$ – otherwise the fringe patterns from different regions of the source manifest their maxima at quite different locations, and the fringe visibility (contrast) is much reduced. We have already estimated (§5) the thickness of the high intensity limb of the source to be $`4\times 10^{11}`$ cm, so with $`D_s4`$ Gpc we can expect interference phenomena to be important only for $`\lambda 10^5`$ cm, i.e. longward of the far ultraviolet region, for the source model we have employed.

It is worth noting that the presence of a caustic network on nano-arcsecond scales has no immediate implications for the observable properties of the GRB afterglow, because at this late stage in its evolution the observable emission from the blast-wave extends over angular scales which are so large that the net magnification is very close to the mean magnification. Microlensing of the afterglow is, however, possible (Loeb & Perna 1998; Mao & Loeb 2001), because the characteristic angular scale of the magnification pattern is much larger in this case.

The two main discrepancies, which we noted in §4, between the modelled nanolensing variations and the observed temporal structure in GRB light-curves are the small amplitude and long time-scale of the nanolensing fluctuations. Given that the model parameters appropriate to real sources are typically not known with any accuracy, one could try to evade the second of these discrepancies by considering sources at high redshift. Specifically, as the source redshift approaches $`z_s=10.6`$, corresponding to nanolensing optical depth $`\tau =1`$, the caustic density becomes very large indeed, and the variability time-scale correspondingly decreases (§3.4). Thus the observed characteristic time-scale of order one second between “pulses” (N96) matches the nanolensing variability time-scale only for very low mass lenses ($`10^8\mathrm{M}_{}`$) at $`z_s=1`$, but for sources at $`z_s10`$ the implied lens mass would correspond to that ($`10^4\mathrm{M}_{}`$) deduced from the data on quasars (Schneider 1993; Minty 2001). This idea is, however, quickly disposed of because as the angular scale of the caustic network shrinks so does the width of the high magnification regions close to the caustics. In turn this means that a source of given size will exhibit smaller variations as it crosses the caustics, exacerbating the other notable discrepancy between the nanolensing simulations and the observed variability. This aspect of gravitational lensing at high optical depth, where the depth of modulation vanishes for non-point-like sources as $`\tau 1`$, is well known (Deguchi and Watson 1987).

We wish to draw attention to the fact that all of the results reported in the present paper are for shear-free environments, and if a large external beam shear is present then the properties of any nanolensing can be quite different. High shear is expected if, for example, nanolensing occurs close to a microlensing caustic. Calculations appropriate to this circumstance will be reported elsewhere (Walker and Lewis 2003).

## 7 Conclusions

Motivated by an existing interpretation of quasar variability in terms of gravitational nanolensing, we have examined the implications of this model for the observed properties of GRBs. Using a self-similar blast-wave model to represent the source (no intrinsic variability) we find that the light-curves of some of the caustic crossings resemble those of the “pulses” commonly seen in GRBs, and for high redshift sources the time-scale of the predicted nanolensing variations is consistent with the GRB data. However, the predicted depth of nanolensing modulation is far too small to explain the deep variations observed in GRBs, and this problem is exacerbated if the GRBs are at high redshift. These results mean that the GRB data do not exclude the nanolensing interpretation of quasar variability; conversely, the simplest (shear free) nanolensing model cannot explain the observed GRB variability. Despite this failure, there are weak indications in the published InterPlanetary Network data that nanolensing may actually be responsible for some of the observed variations of GRBs: the light-curves for one GRB show hints of parallax. This effect is uniquely associated with lens-induced variations, and thus motivates a careful examination of existing IPN data.

MAW much appreciates the hospitality of the Raman Research Institute, where this work was begun. We have benefitted greatly from many helpful discussions with Dipankar Bhattacharya, Sunita Nair and P. Sreekumar, all of whom came within a hair’s breadth of being co-authors on the present paper. GFL would like to thank Joachim Wambsganss for providing his ray-tracing code and David Bowie for his Low album. We thank the referee, Robert Nemiroff, for helpful comments.

## Appendix A Gamma-ray burst astrometry

In this Appendix we give details of the astrometric analysis of the sample of bursts referred to in §5.1. The sample was drawn from the nine bursts observed by BATSE, Ulysses and the Mars Observer spacecraft (Laros et al 1997), plus the 37 bursts observed by BATSE, Ulysses and Pioneer Venus Orbiter (Laros et al 1998). BATSE was capable of localising sources in its own right, by comparing the count rates recorded in different detectors, each pointing in a different direction. This was not, however, a simple facility to implement in practice (Pendleton et al 1999) and the resulting localisations are only accurate to a few degrees — this point has been directly verified by comparing the positions of solar flares, as located by BATSE, with the position of the Sun (Brock et al 1992). Consequently although BATSE localisations are useful in lifting the degeneracy between the two intersection points of the IPN loci, they are not accurate enough to reveal systematic errors in the triangulation procedure at the arcminute level we are interested in.<sup>1</sup><sup>1</sup>1At this point we should draw attention to the anomalous case of BATSE 2475, reported by Laros et al (1997), which apparently displayed an IPN intersection that was significantly discrepant even with respect to the rather crude BATSE localisation. Having no way of understanding this result, Laros et al (1997) were obliged to conclude that Mars Observer had not detected the GRB detected by the other spacecraft, and had in addition responded to another (unspecified) event, which in turn was not detected by the other spacecraft in the network, despite having the characteristics of a GRB. This coincidence seems rather contrived. An alternative explanation for the timing/position anomaly of this burst can be found in terms of parallax (§5 of the present paper), provided that BATSE 2475 exhibited a low expansion speed, with $`\beta _{ap}40`$. (This conclusion assumes a low redshift burst, which nevertheless manifests a caustic crossing event, and makes use of the BATSE localisation given by Meegan et al .) However, we note that the reported discrepancy between BATSE and Mars Observer count rates for BATSE 2475 cannot be so readily explained. To constrain any such errors we therefore employed the localisations which were obtained by the WATCH experiment, quite independently of the triangulated positions. This experiment employed a rotation modulation collimator to locate sources, typically with accuracies better than half a degree (at the $`1\sigma `$ level; Sazonov et al 1998). The additional requirement that each burst have been observed by WATCH necessarily decreases the sample size and we are left with only nine bursts in our final sample; these bursts and the corresponding astrometric information are listed in table A1, while the localisations themselves are presented graphically in figure 12.

It is immediately apparent from table A1 that the WATCH error circles are too large to determine positional errors at the level of a few arcminutes (which is the upper end of the range expected for $`\delta \alpha `$ in the blast-wave model — see §5.1). However, as can be seen from figure 12, there are a number of instances in this sample where the loci derived from burst timing intersect at a very acute angle, and in these cases the intersection point of the loci is displaced on the sky by an amount $`\delta \alpha `$, so that errors of the anticipated size might possibly be revealed by comparison with the WATCH localisations. To test this possibility we have employed a penalty, $`\mathrm{\Pi }(\delta \alpha )`$, which is a function of the offset, $`\delta \alpha `$, in the radius, $`\alpha `$, of one of the loci:

$$\mathrm{\Pi }(\delta \alpha )=\left(\frac{\delta \alpha }{\sigma _\alpha }\right)^2+\left(\frac{\xi }{\sigma _w}\right)^2,$$
(A1)
where $`\xi =\xi (\delta \alpha )`$ is the angular separation between the intersection of the IPN loci and the WATCH location; $`\sigma _w`$ is the estimated error in the WATCH location (one third of the $`3\sigma `$ statistical error, plus the systematic error, as quoted by Sazonov et al 1998), and $`\sigma _\alpha `$ is the estimated error in the determination of $`\alpha `$ (one third of the $`3\sigma `$ statistical-plus-systematic error quoted by Laros et al 1997, 1998). The location ($`\delta \alpha _{min}`$) of the minimum value, $`\mathrm{\Pi }_{min}`$, of this penalty function can, in some cases, give us an estimate of the value of $`\delta \alpha `$. One point to note is that the WATCH error distribution is ellipsoidal, but we are approximating it by a circular distribution.

Generally speaking, of the three loci derived from burst timing, one (corresponding to the BATSE-Ulysses baseline) has by far the smallest errors and our adopted procedure is therefore to treat this locus as an absolute constraint on the burst location. We then discard the annulus with the largest estimated error, and use the remaining locus to determine the intersection point. The radius of this second locus is then allowed to vary by an amount $`\delta \alpha `$, over a range $`\pm 5\sigma _\alpha `$, and for each value of $`\delta \alpha `$ we compute the separation, $`\xi `$, and thus the value of the penalty, $`\mathrm{\Pi }(\delta \alpha )`$. In general the value of $`\mathrm{\Pi }_0\mathrm{\Pi }(0)`$ will not coincide with the minimum value of the function, $`\mathrm{\Pi }_{min}=\mathrm{\Pi }(\delta \alpha _{min})`$, and if $`\mathrm{\Pi }_0\mathrm{\Pi }_{min}1`$ then it is clear from the definition of $`\mathrm{\Pi }`$ that this means $`\delta \alpha _{min}`$ is distinguishable from zero. The ratio $`\delta \alpha _{min}/\sigma _\alpha `$ tells us whether the actual error which we have detected is significantly larger than the expected error in $`\alpha `$. The results of this procedure are given in table A1, where it can be seen that there are four bursts (BATSE 451, 907, 1698 and 2387) for which $`\mathrm{\Pi }_0`$ is significantly greater than $`\mathrm{\Pi }_{min}`$. With the possible exception of BATSE 2387, the errors which we measure in these cases (i.e. $`\delta \alpha _{min}`$) are only as large as one would expect, given the estimated error, $`\sigma _\alpha `$, in each case. Of course it remains possible that the systematic error we are looking for is in fact a major contribution to the estimated value of $`\sigma _\alpha `$, and we are unable to exclude this possibility with the information at hand.


format: text

text:
# Monte Carlo Studies of Three-Dimensional Bond-Diluted Ferromagnets

## 1 Introduction

The influence of quenched, random disorder on phase transitions is of great importance in a large variety of fields cardybook , ranging from experiments with absorbed monolayers Schwenger94  in condensed matter physics to conceptual questions in non-perturbative quantum gravity des .

For pure systems exhibiting a continuous phase transition, Harris Harris74  derived the criterion that random disorder is a relevant perturbation when the critical exponent of the specific heat of the pure system is positive, $`\alpha _{\mathrm{pure}}>0`$. In this case one expects that the system falls into a new universality class with critical exponents governed by a “disordered” fixed point. For $`\alpha _{\mathrm{pure}}<0`$ the behaviour of the pure system should persist, and $`\alpha _{\mathrm{pure}}=0`$ is a special, marginal case. In two dimensions (2D) this scenario has been confirmed by various methods for many different systems selke . In three dimensions (3D) extensive computer simulation studies have concentrated mainly on the site-diluted Ising model selke ; Folk2001 .

If a pure system with a first-order phase transition is subject to quenched disorder, the transition is softened and may even turn into a continuous one ImryWortis79 . This is always the case in 2D Aizenman89 ; for numerical verifications see Refs. \[9-13\]. In higher dimensions, a tricritical point may appear at a finite concentration of impurities cardy , separating “non-softened” first-order and “softened” second-order regimes cardy\_statphys . Numerically such a scenario has recently been observed for the 3D site-diluted 3-state Potts model Ballesteros00 . Since the first-order transition of the pure version of this model is very weak villanova , however, the characterization of the tricritical point remained inconclusive.

In this report we give an overview on recent results obtained from extensive Monte Carlo (MC) computer simulations of the 3D 2-state (Ising) ccp\_ising  and 4-state Potts \[19-21\] models with bond-dilution. The motivation to study the 4-state Potts model derives from the fact that, in the pure case, this model is known to exhibit a fairly strong first-order transition, such that a disorder-induced softening to a second-order transition would give clear support for the theoretical picture sketched above. Modeling the disorder by bond-dilution enables in the Ising case a test of the expected universality with respect to the type of disorder. Furthermore, for both models this choice facilitates a quantitative comparison with recent high-temperature series expansions meik1 ; meik2  for general random-bond $`q`$-state Potts models.

## 2 Model and Simulation Setup

The 3D bond-diluted $`q`$-state Potts model is defined by the Hamiltonian

$$\beta H=\underset{ij}{}K_{ij}\delta _{\sigma _i,\sigma _j};\sigma _i=1,\mathrm{},q,$$
(1)
where the sum extends over all pairs of neighbouring sites on a cubic lattice of size $`L^3`$ with periodic boundary conditions, and the couplings $`K_{ij}`$ are distributed according to the distribution

$$\mathrm{}(K_{ij})=p\delta (K_{ij}K)+(1p)\delta (K_{ij}),$$
(2)
where $`KJ/k_BT`$ is the inverse temperature in natural units. The dilution parameter $`p`$ is thus the concentration of magnetic bonds in the system, i.e., $`p=1`$ corresponds to the pure case. Below the percolation threshold $`p_c=\mathrm{0.248\hspace{0.17em}812\hspace{0.17em}6}(5)`$ Lorenz1998  one does not expect any finite-temperature phase transition since without any percolating cluster in the system long-range order is impossible.

The model (1), (2) was studied by means of large-scale MC simulations using the Swendsen-Wang (SW) cluster algorithm SW87  in the regime of second-order transitions, and multicanonical simulations muca1 ; muca2  in the regime where the first-order transition of the pure 4-state Potts model persists, i.e., at weak dilution close to $`p=1`$. Thermodynamic quantities were averaged over a large number of quenched disorder realisations, ranging between 2 000 and 5 000. The stability of the disorder averages has been checked by monitoring running averages as a function of the number of random samples. In fact, some care is necessary because too small a number of disorder realisations would lead to typical values rather than average ones Derrida84etc , and these two values are different if the probability distribution over the disorder realisations exhibits a long tail.

## 3 Results

### 3.1 3D Bond-Diluted Ising Model

The phase diagram of the 3D bond-diluted Ising model has been obtained numerically from the locations of the maxima of a diverging quantity such as the magnetic susceptibility $`\overline{\chi }_L`$ depicted in Fig. 1. We focused on $`\overline{\chi }_L`$ because the stability of the disordered fixed point implies a negative specific-heat exponent in a random system chayes . The error in this quantity is hence typically larger than that in the susceptibility.

For an accurate determination of the maxima of the susceptibility we used the histogram reweighting technique with $`N_{\mathrm{MCS}}=\mathrm{2\hspace{0.17em}500}`$ MC sweeps (MCS) and between $`N_s=\mathrm{2\hspace{0.17em}500}`$ and $`\mathrm{5\hspace{0.17em}000}`$ disorder realisations. The choice of $`N_{\mathrm{MCS}}`$ is justified by the increasing behaviour of the energy autocorrelation time $`\tau _e`$ as a function of $`p`$ and $`L`$. At the critical point of a second-order phase transition one expects a finite-size scaling (FSS) behaviour $`\tau _eL^z`$, where $`z`$ is the dynamical critical exponent. From fits to the data shown in Fig. 2 we obtained values of $`z0.27,0.38,0.41`$, and $`0.59`$ for $`p=0.4,0.55,0.7`$, and 1.0 (the pure case), respectively. The critical slowing down thus weakens for the disordered model and becomes less and less pronounced when the concentration of magnetic bonds $`p`$ decreases. The largest autocorrelation time observed for the disordered model was around $`\tau _e9`$ for $`p=0.7`$ and $`L=96`$. For each dilution and each size we thus collected at least $`250`$ effectively uncorrelated measurements of the physical quantities, i.e., $`N_{\mathrm{MCS}}>250\tau _e`$. On the other hand, for small $`pp_c`$ it is necessary to increase the number of disorder realisations because of the vicinity of the percolation threshold.

In order to check the quality of the averaging techniques, we first studied the stability of the susceptibility versus the number of MC sweeps involved in the thermal average. For the largest size considered our results are given in Table 1 for different samples as well as for the disorder average. With $`\mathrm{2\hspace{0.17em}500}`$ MCS, the accuracy of the results for a given sample is not perfect, of course, but the precision of the average over disorder is quite good on the other hand. The disorder average procedure has been investigated by computing the susceptibility $`\chi _j`$ for different samples, $`1jN_s`$. As can be seen in Fig. 3, the dispersion of the values of $`\chi `$ is not very large because the fluctuations in the (running) average value disappear already after a few hundreds of realisations.

The phase diagram as obtained from the susceptibility maxima for the largest lattice size is shown in Fig. 4. We do not include the results from high-temperature series expansions meik1  in this figure since they would just fall on top of the MC data. For comparison we have drawn, however, a simple mean-field (MF) estimate of the transition point,

$$K_c^{\mathrm{MF}}(p)=K_c(1)/p,$$
(3)
with $`K_c(1)=\mathrm{0.443\hspace{0.17em}308\hspace{0.17em}8}(6)`$ talapov . This rather crude approximation only holds in the low-dilution regime, $`p>0.8`$. On the other hand, the single-bond effective-medium (EM) approximation of Turban turban ,

$$K_c^{\mathrm{EM}}(p)=\mathrm{ln}\left[\frac{(1p_c)e^{K_c(1)}(1p)}{pp_c}\right],$$
(4)
gives very good agreement with the simulated transition line over the full dilution range. Since $`K_c(1)`$ and $`p_c`$ are input parameters this relation is trivially exact in the vicinity of both the pure system and the percolation threshold.

Due to the competition between different fixed points, one of the main problems encountered in previous studies of the disordered Ising model was the question whether one measures effective or asymptotic exponents. Although the change of universality class should happen theoretically for an arbitrarily low disorder, it can be very difficult to measure the new critical exponents because the asymptotic behaviour cannot always be reached practically. Another difficulty comes from the vicinity of the ratios $`\gamma /\nu `$ and $`\beta /\nu `$ in the pure and disordered universality classes. Indeed, for the 3D Ising model these values are:

$`\gamma /\nu `$ $`=`$ $`1.966(6),\beta /\nu =0.517(3),\nu =0.6304(13)\text{pure case }\text{Guida1998},`$ (5)
$`\gamma /\nu `$ $`=`$ $`1.963(5),\beta /\nu =0.519(3),\nu =0.6837(53)\text{disordered case }\text{Ballesteros1998}.`$ (6)
Recent field theoretical determinations of critical exponents in the disordered case are presented in Refs. Varnashev2000 ; Pelissetto2000 , and for an excellent review of various experimental, theoretical and numerical estimates in the last two decades, see Ref. Folk2001 . Thus, from standard FSS techniques, the critical exponent $`\nu `$ only will allow us to discriminate between the two fixed points. This exponent can be evaluated from the FSS behaviour of the derivative of the magnetisation w.r.t. temperature which is expected to behave as $`d\mathrm{ln}\overline{m}/dKL^{1/\nu }`$. From this power-law behaviour, we have extracted the effective size-dependent exponent $`(1/\nu )_{\mathrm{eff}}`$ which is plotted in Fig. 5 against $`1/L_{\mathrm{min}}`$ for different bond concentrations $`p`$, with $`L_{\mathrm{min}}`$ denoting the smallest lattice size used in the fits.

We clearly see that in the regime of low dilution ($`p`$ close to $`1`$), the system is influenced by the pure fixed point. On the other hand, when the bond concentration is small, the vicinity of the percolation fixed point induces a decrease of $`1/\nu `$ below its expected disordered value. This is plausible since the percolation fixed point is characterized by $`1/\nu 1.12`$ Lorenz1998 .

### 3.2 3D Bond-Diluted 4-State Potts Model

Let us now turn to the 4-state Potts model which exhibits in the pure case a rather strong first-order phase transition. In order to map out the phase diagram of the diluted model we considered all concentrations $`p`$ in the interval $`[0.28,1]`$ in steps of $`0.04`$ and determined again the locations of the maxima of the susceptibility for a given lattice size $`L`$. The resulting phase diagram is depicted in Fig. 6, where we show for comparison again the simple mean-field prediction (3) and the effective-medium approximation (4), using $`K_c(1)=0.62863(2)`$ Chatelain01 . On the scale of Fig. 6, estimates from high-temperature series expansions up to order 18 are hardly distinguishable, for a comparison see Ref. meik2 .

In a second step, the order of the phase transitions was investigated. Here, in order to satisfy our criterion $`N_{\mathrm{MCS}}>250\tau _e`$, the number of MC sweeps had to be increased up to $`\mathrm{15\hspace{0.17em}000}\mathrm{30\hspace{0.17em}000}`$, which is rather large compared to the values used in the Ising case. A first indication is given by the FSS behaviour of the autocorrelation time $`\tau _e`$ at the transition point. A glance on the log-log plot of Fig. 7 shows a crossover around $`p=0.80`$ from a power-law behaviour for strong disorder (small $`p`$) to a clear exponential behaviour for weak disorder ($`p1`$), as is typical for a first-order phase transition. In fact, in the latter case one expects $`\tau _e\mathrm{exp}(2\sigma _{od}L^2)`$, where the (reduced) interface tension $`\sigma _{od}`$ parameterizes the free-energy barrier separating the coexisting ordered and disordered phases.

In the first-order regime we performed multicanonical simulations and estimated the interface tension from

$$\sigma _{od}=\frac{1}{2L^2}\mathrm{log}\frac{P_{\mathrm{max}}}{P_{\mathrm{min}}},$$
(7)
where $`P_{\mathrm{max}}`$ is the maximum of the probability density reweighted to the temperature where the two peaks are of equal height, and $`P_{\mathrm{min}}`$ is the minimum in between, see Fig. 8. The linear extrapolations of $`\sigma _{od}`$ in $`1/L`$ in the lower part of Fig. 8 imply non-vanishing interface tensions only for $`p=0.84`$ and above. For $`p0.76`$, $`\sigma _{od}`$ seems to vanish in the infinite-volume limit, being indicative of the expected softening to a second-order phase transition. The tricritical point would thus be located around $`p=0.760.84`$, in good agreement with the estimate of $`p=0.80`$ derived from the analysis of autocorrelation times.

To confirm the softening for $`p0.76`$ we have performed a detailed FSS study at $`p=0.56`$ with lattice sizes ranging up to $`L=96`$ and the number of realisations varying between 2 000 and 5 000 Chatelain01 . The choice of $`p=0.56`$ is motivated by our observation that in this range of dilutions the corrections to asymptotic FSS of the effective transition points are minimal, cf. Fig. 9.

The log-log plot for $`\overline{\chi }_{\mathrm{max}}`$ in Fig. 10 shows that for this quantity the corrections to asymptotic FSS seem to become quite small above $`L=30`$, and fits of the form $`a_\chi L^{\gamma /\nu }`$ starting at $`L_{\mathrm{min}}>30`$ yield $`\gamma /\nu =1.50(2)`$. Using the data for $`L<30`$ only, on the other hand, we obtained perfect fits assuming percolation exponents, $`\gamma /\nu 2.05`$ Lorenz1998 , cf. Fig. 10. Similarly, the FSS of the quantity $`(d\mathrm{ln}\overline{m}/dK)_{K_{\mathrm{max}}}L^{1/\nu }`$ gives for $`L_{\mathrm{min}}>30`$ an estimate of the exponent $`1/\nu =1.33(3)`$, consistent with the stability condition of the random fixed point ($`1/\nu D/2=1.5`$) chayes . The same procedure was applied to the magnetization $`\overline{m}L^{\beta /\nu }`$, but here the associated critical exponent turned out to be not yet stable. We therefore also considered the FSS behaviour of higher (thermal) moments of the magnetization, $`\overline{\mu ^n}`$, which should scale with an exponent $`n\beta /\nu `$. The results for the first moments exhibit, however, again much stronger corrections to scaling than we observed for $`\overline{\chi }`$ or $`d\mathrm{ln}\overline{m}/dK`$, leading to quite a conservative final estimate of $`\beta /\nu =0.65(5)`$. We nevertheless note that our results do not fit satisfactorily the scaling law $`2\beta /\nu =d\gamma /\nu `$. The reason could be the strong corrections to scaling at the random fixed point which are hard to cope with for medium-sized systems Chatelain01 .

## 4 Conclusions

By performing large-scale Monte Carlo simulations we have investigated the influence of bond dilution on the critical properties of the 3D Ising and 4-state Potts models. In the 3D Ising case the universality class of the disordered model is modified by disorder but its precise characterization turned out be difficult because of the competition between the different fixed points which induce crossover effects, even for relatively large lattice sizes.

Applying similar techniques to the 3D 4-state Potts model we obtained clear evidence for softening to a continuous transition at strong disorder, with estimates for the critical exponents of $`\nu =0.752(14)`$, $`\gamma =1.13(4)`$, and $`\beta =0.49(5)`$ at $`p=0.56`$. The analysis of both the autocorrelation time and the interface tension leads to the conclusion of a tricritical point around $`p=0.80`$.

## 5 Acknowledgements

We would like to thank Meik Hellmund and Loic Turban for helpful discussions. The mutual visits within this collaboration were financially supported by the joint PROCOPE exchange programme of the DAAD and EGIDE. C.C. thanks the EU network “EUROGRID: Discrete Random Geometries: From solid state physics to quantum gravity” for a post-doctoral position in Leipzig, and W.J. thanks the German-Israel-Foundation (GIF) for support. The numerical work would have been impossible without the computer-time grants 2000007 of the Centre de Ressources Informatiques de Haute-Normandie (CRIHAN), hlz061 of NIC, Jülich, and h0611 of LRZ, München. We are grateful to all institutions for their generous support.

## Index