format: text

text:
# Poincaré cycle of a multibox Ehrenfest urn model with directed transport

## I introduction

Physical laws governing the microscopic processes are mostly reversible in time. In macroscopic world, however, people often experience time-irreversible phenomena in their daily life. To understand why the reversible microscopic processes lead to irreversible macroscopic manifestations one refers to the Poincaré Theorem, which states that a system having a finite energy and confined to a finite volume will, after a sufficient long time – the so called Poincaré cycle, return to an arbitrarily small neighborhood of almost any given initial state . The key point is to note that the typical value of a Poincaré cycle for even a moderate-sized system is far beyond the meaningful time scale one can measure or experience, thus the irreversiblity is realized.

Usually to describe a macroscopic system one has to know only a few parameters, such as volume, pressure, and temperature. However, to describe the same system in terms of its microscopic constituents, one has to deal with a large number of parameters, such as the momenta and positions of a huge amount of particles, which are impossible to calculate in practice. Based on this reason together with the fact that the macroscopic laws are insensitive to the microscopic details (of system history), it is natural for people to adopt the probability (ensemble) description in statistical mechanics, which deals with the equilibrium state (a macroscopic state that has stationary value of state parameters) of a macroscopic system. In this kind of description the macroscopic quantities are defined as the ensemble average of their microscopic correspondences. This definition connects the microscopic and macroscopic worlds.

To study how a system approaches its equilibrium state one also uses probability description, where the evolution of the system is treated as a stochastic process. One famous model for simulating such a process was proposed by Ehrenfest one century ago , which is an $`N`$-ball, $`2`$-urn problem. In the beginning $`N`$ numbered balls are distributed arbitrarily in either urn $`A`$ or urn $`B`$. At each time step one ball is picked out at random and then put into the other urn. This simple model can be exactly solved to give an explicit Poincaré cycle. This model was then generalized by several authors to mimic more complicated situations encountered in real physical phenomena . An attractive feature of these urn model problems is that they are easy to formulate but not always easy to solve. The solutions obtained have, therefore, sometimes led to new mathematical techniques and insights . Recently, some new urn models were proposed and solved analytically or numerically. Their results provide very good descriptions on granular and glass systems .

In this paper, we obtain the exact solution of a generalized urn model. Hereafter we call it “periodic urn model”. In this model, one considers $`N`$ distinguishable balls which are distributed in $`M`$ urns. These $`M`$ urns are arranged along a circle and numbered one by one to form a cycle, that is, we define the $`(M+1)`$th urn as the 1st urn (See Fig. 1). To begin with, the initial distribution of the $`N`$ balls in the $`M`$ urns is given by $`|m_{1,0},m_{2,0},\mathrm{},m_{M,0}|𝐦_0`$, where $`m_{i,0}`$ is the number of balls in the $`i`$th urn at the start. At each time step one ball is picked out of the $`N`$ balls such that every ball has an equal probability of being picked up. The ball is then placed into the next numbered urn. The state that the $`i`$th urn contains $`m_i`$ balls is represented by $`|m_1,m_2,\mathrm{},m_M|𝐦`$, which we name it state vector. Hereafter we call a distribution string $`𝐦`$ (without knowing the numbering of the balls) a configuration of the system. Otherwise, if we also know the location of each numbered ball, we call such a distribution a microstate of the system.

After $`s`$ steps, the transition probability from state $`|𝐦_0`$ to state $`|𝐦`$ can be written as $`𝐦|P^s|𝐦_0`$, where $`P`$ represents the operation in one step, and the orthonomality conditions of these state vectors have been assumed.

According to the above description, the transition probabilities corresponding to the $`s`$th step and the $`(s1)`$th step satisfy the recursion relation

$`m_1,m_2,\mathrm{},m_M|P^s|𝐦_0`$ (1)
$`={\displaystyle \underset{i=1}{\overset{M}{}}}{\displaystyle \frac{m_i+1}{N}}\mathrm{},m_i+1,m_{i+1}1,\mathrm{}|P^{s1}|𝐦_0,`$ (2)
where $`m_{M+1}=m_1`$ as has been mentioned before. In addition, any state that does not satisfy the constraint

$$m_1+m_2+\mathrm{}+m_M=N.$$
(3)
is an unphysical state and has null contribution in the sum. Hereafter we will use Eq. (2) as the basis to derive all the results we want to know.

This paper is organized as follows. In the next section we calculate the average number of balls in an urn at any time. In section III we introduce a generating function that has $`N`$ variables and solve the problem completely. In section IV the solution of the model will be applied to the calculation of Poincaré cycle. Finally in section V we give the summary of this paper.

## II Average number of balls in an urn

The first thing we want to know is how many balls on average will appear in the $`i`$th urn after first $`s`$ steps. We define the average (the expectation value) of a quantity $`A`$ (which depends on the state vector $`|𝐦`$ at each step, written as $`A(𝐦)`$) after $`s`$ steps as

$$A_s=\underset{\{𝐦\}}{}A(𝐦)𝐦|P^s|𝐦_0,$$
(4)
where $`\{𝐦\}`$ include all the configurations satisfying the constraint (3).

Let $`A(𝐦)=m_i`$, then from Eq. (2) and (4) we have

$`m_i_s`$ $`=`$ $`{\displaystyle \underset{\{𝐦\}}{}}m_i𝐦|P^s|𝐦_0,`$ (5)
$`=`$ $`{\displaystyle \underset{\{𝐦\}}{}}{\displaystyle \underset{j=1}{\overset{M}{}}}{\displaystyle \frac{m_i(m_j+1)}{N}}`$ (7)
$`\times \mathrm{},m_j+1,m_{j+1}1,\mathrm{}|P^{s1}|𝐦_0,`$
$`=`$ $`{\displaystyle \underset{j=1}{\overset{M}{}}}\left({\displaystyle \frac{m_im_j_{s1}}{N}}{\displaystyle \frac{m_j_{s1}}{N}}\delta _{j,i}+{\displaystyle \frac{m_j_{s1}}{N}}\delta _{j+1,i}\right),`$ (8)
$`=`$ $`{\displaystyle \underset{j=1}{\overset{M}{}}}{\displaystyle \frac{m_im_j_{s1}}{N}}{\displaystyle \frac{m_i_{s1}}{N}}+{\displaystyle \frac{m_{i1}_{s1}}{N}},`$ (9)
$`=`$ $`\left(1{\displaystyle \frac{1}{N}}\right)m_i_{s1}+{\displaystyle \frac{m_{i1}_{s1}}{N}},`$ (10)
here we have used the constraint (3).

Now we are ready to solve $`m_i_s`$. The recurrence relation (10) can be written as

$$_s=P_{ave}_{s1},$$
(11)
where $`_s`$ is a $`M\times 1`$ column vector defined by

$$_s=\left[\begin{array}{cccc}m_1_s& & & \\ m_2_s& & & \\ \mathrm{}& & & \\ m_M_s& & & \end{array}\right],$$
(12)
and $`P_{ave}`$ is a $`M\times M`$ matrix written as

$$P_{ave}=\left[\begin{array}{cccc}1\frac{1}{N}& 0& \mathrm{}& \frac{1}{N}\\ \frac{1}{N}& 1\frac{1}{N}& \mathrm{}& 0\\ \mathrm{}& \mathrm{}& \mathrm{}& \mathrm{}\\ 0& 0& \mathrm{}& 1\frac{1}{N}\end{array}\right].$$
(13)
By means of the recurrence relation (11), $`_s`$ can be deduced:

$$_s=P_{ave}^s_0,$$
(14)
where

$$_0=\left[\begin{array}{cccc}m_{1,0}& & & \\ m_{2,0}& & & \\ \mathrm{}& & & \\ m_{M,0}& & & \end{array}\right]$$
(15)
represents the initial state. $`P_{ave}^s`$ can be calculated if one knows the eigenvalues $`\lambda _m`$ and eigenvectors $`Q_m`$ of $`P_{ave}`$. They are given by (See Fig. 2)

$$\lambda _m=1\frac{1}{N}+\frac{1}{N}q_m^{},Q_m=\frac{1}{\sqrt{M}}\left[\begin{array}{cccc}q_m& & & \\ q_m^2& & & \\ \mathrm{}& & & \\ q_m^M& & & \end{array}\right],$$
(16)
where

$$q_m=\mathrm{exp}\left(\frac{2m\pi i}{M}\right),m=1,2,\mathrm{},M.$$
(17)
Denote $`R`$ as the $`M\times M`$ matrix of the eigenvectors $`Q_m`$

$`R`$ $`=`$ $`[Q_1,Q_2,\mathrm{},Q_M]`$ (18)
$`=`$ $`{\displaystyle \frac{1}{\sqrt{M}}}\left[\begin{array}{cccc}q_1& q_2& \mathrm{}& q_M\\ q_1^2& q_2^2& \mathrm{}& q_M^2\\ \mathrm{}& \mathrm{}& \mathrm{}& \mathrm{}\\ q_1^M& q_2^M& \mathrm{}& q_M^M\end{array}\right]`$ (24)
and $`\mathrm{\Lambda }`$ as the diagonal matrix of $`P_{ave}`$’s eigenvalues $`\lambda _m`$

$$\mathrm{\Lambda }=\left[\begin{array}{cccc}\lambda _1& 0& \mathrm{}& 0\\ 0& \lambda _2& \mathrm{}& 0\\ \mathrm{}& \mathrm{}& \mathrm{}& \mathrm{}\\ 0& 0& \mathrm{}& \lambda _M\end{array}\right]$$
(25)
then we obtain

$$P_{ave}^s=R\mathrm{\Lambda }^sR^1=R\mathrm{\Lambda }^sR^{}=R\mathrm{\Lambda }^sR^{},$$
(26)
where we have used the following properties of $`R`$

$$R=R^t,(R^1)=R^{}=(R)^{},R_{mn}=q_n^m=q_1^{mn}.$$
(27)
Now the average number of balls in the $`i`$th urn after $`s`$ steps can be determined:

$`m_i_s`$ $`=`$ $`{\displaystyle \frac{1}{M}}{\displaystyle \underset{j=1}{\overset{M}{}}}{\displaystyle \underset{k=1}{\overset{M}{}}}{\displaystyle \underset{l=1}{\overset{M}{}}}q_j^i\lambda _j^s\delta _{jk}q_l^km_l_0,`$ (28)
$`=`$ $`{\displaystyle \frac{1}{M}}{\displaystyle \underset{j=1}{\overset{M}{}}}{\displaystyle \underset{l=1}{\overset{M}{}}}q_1^{j(il)}\lambda _j^sm_{l,0},`$ (29)
where $`m_l_0=m_{l,0}`$ is the initial number of balls in the $`l`$th urn.

Let us now consider a simple example. Suppose initially all the $`N`$ balls are in the first urn, that is,

$$m_{1,0}=N,m_{2,0}=m_{3,0}=\mathrm{}=m_{M,0}=0,$$
(30)
then according to Eq. (29), we have

$$m_1_s=\frac{N}{M}\underset{j=1}{\overset{M}{}}\lambda _j^s$$
(31)
Fig. 3 shows the results for $`N=50`$ at $`M=2,5,10,25`$. The $`M=2`$ case is the original Ehrenfest model, in which the average number of balls in the first urn decays to $`N/2`$ in a period of steps of order $`N`$. For any $`M>2`$ case, however, we observe that before the system arrives its true equilibrium (here we mean the value of $`m_i_s`$ for each $`i`$ does not change anymore), $`m_1_s`$ undergoes several oscillations, which seems to be a unique feature of this model and have never been found in other kinds of urn models — to our knowledge. Furthermore, in Fig. 3 the $`M=25`$ case shows that before the appearance of the first peak of $`m_1_s`$ there is a period during which $`m_1_s`$ is almost zero. This phenomenon together with the oscillations mentioned before seem to be typical results when both $`M`$ and $`N`$ are large.

To understand the behaviors of $`m_1_s`$ when $`N`$ is large, note first that in this limit Eq. (31) can be approximated by

$$m_1_s=\frac{N}{M}\underset{j=1}{\overset{M}{}}\mathrm{exp}\left[\frac{s(q_j^11)}{N}\right],$$
(32)
hence $`m_1_s/N`$ becomes a universal function of $`s/N`$ for a fixed $`M`$. In Fig. 4 we plot the $`m_1_s/N`$ curves for $`M=30`$ (which is large enough in practice) and $`N=5,15,30,60`$. The $`N=5`$ case is shown here too to compare with those cases with large $`N`$. As one can see, except for the $`N=5`$ case, the $`m_1_s/N`$ curves corresponding to different $`N`$’s merge and become one universal function of $`s/N`$. In the beginning this function decays from 1 exponentially like that in the $`M=2`$ case. However, this curve does not decay to its stationary value $`1/M`$ directly, instead, before it growing up again during a period it becomes too tiny to be visible. In practice this tiny value can be ignored thus hereafter we treat this period as the “empty-urn” period. The values of $`m_1_s/N`$ in this period for our chosen cases are of the order $`10^4`$. The precise duration of the empty-urn period depends on how accurate we treat the urn as empty. After this period the curve grows up again and oscillates several times with evanescent amplitudes and finally approaches the stationary mean value $`1/M`$.

Further exploration shows that there can have more than one empty-urn period. And the total number of these periods, their durations, and the local maxima and minima of each $`m_1_s/N`$ curve are all determined by $`M`$. Fig. 5 shows the results for $`N=60`$ at $`M=30,60,90`$. These three $`M`$’s are chosen intentionally to have a simple integer proportion $`1:2:3`$ in order to help us inspect relevant results without much effort. The first maxima of $`m_1_s/N`$ for these three $`M`$’s are located at $`s/N=30,60,90`$, which are exactly these $`M`$’s. In addition, the locations of the first maxima for $`M=60`$ and $`90`$ are just the same locations of the second and third maxima for $`M=30`$, respectively. Moreover, for the $`M=60`$ and $`M=90`$ cases the total number of empty-urn periods are more than one.

Starting from Eq. (32), we now derive two useful approximations of $`m_1_s/N`$ to help us understand these observations. First, define

$$\tau \frac{s}{N}$$
(33)
and expand the exponential functions in Eq. (32) as power series of $`\tau `$, we have

$`{\displaystyle \frac{m_1_s}{N}}`$ $`=`$ $`{\displaystyle \frac{e^\tau }{M}}{\displaystyle \underset{j=1}{\overset{M}{}}}\left(1+\tau q_j^1+{\displaystyle \frac{1}{2}}\tau ^2q_j^2+\mathrm{}\right)`$ (34)
$`=`$ $`e^\tau \left[1+{\displaystyle \frac{\tau ^M}{M!}}+{\displaystyle \frac{\tau ^{2M}}{(2M)!}}+\mathrm{}\right].`$ (35)
In deriving Eq. (35), we have used the fact that $`_{j=1}^Mq_j^k=_{j=1}^Mq_1^{jk}=0`$ except for $`k=rM`$ with $`r`$ an integer. Note that the form of the $`j`$th term ($`j>0`$) appearing in the last line of Eq. (35) is the same as the probability of $`n`$ successful trials in a Poisson process :

$$P(n)=\frac{e^\tau \tau ^n}{n!}.$$
(36)
Here $`n=jM`$, and $`\tau =n`$ is the expectation value of $`n`$. This is not an accident and can be easily understood. The quantity $`m_1_s/N`$ represents not only the average number of balls in the first urn divided by $`N`$ but also the probability of finding a certain ball, say, ball 1, in the first urn. At each time step ball 1 has the probability of $`p=1/N`$ being picked out and moved to the next urn. Now since $`N`$ is large, $`p`$ is small. In this limit if we do the same operation $`s`$ times (here we assume $`s`$ is also large and define $`\tau =ps=s/N`$), then the probability for ball 1 to move $`n`$ steps forward ($`n`$ successful trials) from the first urn is given by Eq. (36). Furthermore, since our system has a circulating property, the probability of finding ball 1 in the first urn after $`s`$ steps consists of the following possibilities: (1) ball 1 has never been selected (the corresponding probability is $`(1p)^se^{ps}=e^\tau `$), and (2) it has been picked up $`M`$ times, and (3) it has been chosen $`2M`$ times, and so on. The summation of all these contributions gives us the expression of the last line of Eq. (35).

Applying the saddle-point method to each $`e^\tau \tau ^{jM}`$ term and Stirling formula  to each $`(jM)!`$ term, we obtain

$$e^\tau \tau ^{jM}\left(\frac{jM}{e}\right)^{jM}e^{(\tau jM)^2/2jM},$$
(37)
$$(jM)!=\sqrt{2\pi jM}\left(\frac{jM}{e}\right)^{jM}.$$
(38)
We thus have

$$\frac{m_1_s}{N}e^\tau +\frac{e^{(\tau M)^2/2M}}{\sqrt{2\pi M}}+\frac{e^{(\tau 2M)^2/4M}}{\sqrt{4\pi M}}+\mathrm{}.$$
(39)
Eq. (39) tells us that: (1) The $`m_1_s/N`$ curve consists of an exponentially decaying term $`e^\tau `$ and a series of Gaussian terms of different heights. (2) The center of the $`j`$th Gaussian term is located at $`jM`$, which is also a good approximation to the location of the $`j`$th local maximum of the $`m_1_s/N`$ curve. The separation between two successive maxima is thus $`\mathrm{\Delta }\tau =M`$. (3) The height and the standard deviation of the $`j`$th Gaussian term are $`1/\sqrt{2\pi jM}`$ and $`\sqrt{jM}`$, respectively. (4) If $`j<<M`$ we have $`\sqrt{jM}<<\mathrm{\Delta }\tau `$, and the overlap between two successive Gaussian terms can be neglected. In this situation the $`m_1_s/N`$ curve in the neighborhood of $`jM`$ can be well approximated by the $`j`$th Gaussion term. Similarly, the $`j`$th local minimum between the $`j`$th and $`(j+1)`$th local maxima are approximated by

$$\left(\frac{1}{\sqrt{j}}+\frac{1}{\sqrt{(j+1)}}\right)\frac{e^{M/4(2j+1)}}{\sqrt{2\pi M}}.$$
(40)
When $`j`$ is too large, Eq. (39) becomes inaccurate, and we must find another expression for better descriptions. Note that Eq. (32) can be rewritten as

$`{\displaystyle \frac{m_1_s}{N}}`$ $`=`$ $`{\displaystyle \frac{1}{M}}[1+2e^{\tau (1\mathrm{cos}\theta _1)}\mathrm{cos}(\tau \mathrm{sin}\theta _1)`$ (42)
$`+2e^{\tau (1\mathrm{cos}\theta _2)}\mathrm{cos}(\tau \mathrm{sin}\theta _2)+\mathrm{}],`$
where we have defined

$$\theta _j=\frac{2\pi j}{M}=j\theta ,\theta =\frac{2\pi }{M}.$$
(43)
Remember that $`M`$ is a large number, so

$$\mathrm{sin}\theta _1\theta ,\mathrm{\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}1}\mathrm{cos}\theta _1\frac{\theta ^2}{2}.$$
(44)
We thus have the approximation

$$\frac{m_1_s}{N}\frac{1+2e^{\tau \theta ^2/2}\mathrm{cos}\left(\tau \theta \right)+2e^{2\tau \theta ^2}\mathrm{cos}\left(2\tau \theta \right)}{M}$$
(45)
for describing the behaviors of $`m_1_s/N`$ curve in the region $`\tau >>M^2/2\pi ^2`$, that is, $`\tau (1\mathrm{cos}\theta _1)\tau \theta ^2/2>>1`$. From this expression we also get the result

$$\underset{s\mathrm{}}{lim}m_1_s=\frac{N}{M},$$
(46)
which must be true after the system has reached its stationary state.

Up to now we have been focusing our attention on only one single urn. In fact, we can also understand the behaviors of the system in a global manner. First, let’s define the “phase angle” of the $`k`$th urn (see Fig. 1) as

$$\varphi _k=(k1)\theta .$$
(47)
Now we define the “center-of-mass” (COM) of the $`N`$-ball, $`M`$-urn system as

$$\text{COM}\underset{k=1}{\overset{M}{}}e^{i\varphi _k}m_k_s.$$
(48)
According to Eq. (48) the COM is in general a complex number, say, $`Z=Re^{i\varphi }`$. Here $`R=|Z|`$ satisfies $`0RN`$, and $`\varphi `$ is the phase angle of $`Z`$. The variation of the the phase angle $`\varphi `$ with respect to $`s`$ or $`\tau =s/N`$ represents how fast on average these $`N`$ balls circulate around the $`M`$ urns, and the norm $`R`$ here gives us the information of the extent of the distribution of the $`N`$ balls. If we consider the normalized center-of-mass $`\mathrm{COM}/Nre^{i\varphi }`$ instead of COM, we conclude that the trace of $`re^{i\varphi }`$ must be contained inside the unit circle. Substitute (31) and (47) into (48) and use (16) and (17), we get

$`{\displaystyle \frac{\text{COM}}{N}}`$ $`=`$ $`{\displaystyle \frac{1}{M}}{\displaystyle \underset{j=1}{\overset{M}{}}}{\displaystyle \underset{k=1}{\overset{M}{}}}q_1^{(j1)(k1)}\lambda _j^s`$ (49)
$`=`$ $`\lambda _1^s`$ (50)
$``$ $`e^{\tau \theta ^2/2}e^{i\tau \theta },`$ (51)
hence

$$r=e^{\tau \theta ^2/2},\varphi =\tau \theta .$$
(52)
Here we see that the $`\mathrm{COM}/N`$ curve is approximately described by a spiral circulating inside a unit circle. The angular frequency of the circulating motion of the COM with respect to $`s`$ is $`\theta =2\pi /M`$ (clockwise), consistent with both the definition of our game and the oscillation behaviors of the $`\mathrm{COM}/N`$ curve discussed before. Furthermore, when $`\tau >M^2/2\pi ^2`$ we have $`r<e^1`$, which indicates that the system has a wide extent of the ball-distribution, also consistent with the previous discussions. Some examples are illustrated in Fig 6 and Fig. 7.

## III State matrix and generating function

Now we calculate $`𝐦|P^s|𝐦_0`$ – the transition probability from $`|𝐦_0`$ to $`|𝐦`$ after $`s`$ steps. Once one knows the exact solution of $`𝐦|P^s|𝐦_0`$, any quantity can be calculated explicitly.

Define

$$S_{\mathrm{𝐦𝐦}^{}}=𝐦|P|𝐦^{},$$
(53)
then we have:

$$𝐦|P^s|𝐦_0=(S^s)_{\mathrm{𝐦𝐦}_0}.$$
(54)
Here $`S`$ is a $`H_N^M\times H_N^M`$ matrix, we name it state matrix, and $`|𝐦`$ is a $`H_N^M`$ column vector, here

$$H_N^M=C_{M1}^{N+M1}=\frac{(N+M1)!}{N!(M1)!}.$$
(55)
Like before, $`S^s`$ can be calculated by means of its eigenvalues and eigenstates. According to Eq. (2), the matrix $`S`$ has components

$`S_{\mathrm{𝐦𝐦}^{}}`$ $`=`$ $`{\displaystyle \underset{i=1}{\overset{M}{}}}{\displaystyle \frac{m_i+1}{N}}\delta _{m_1,m_1^{}}\delta _{m_2,m_2^{}}\mathrm{}`$ (57)
$`\times \delta _{m_i+1,m_i^{}}\delta _{m_{i+1}1,m_{i+1}^{}}\mathrm{}\delta _{m_M,m_M^{}},`$
where $`m_{M+1}=m_1`$ have been assumed. The eigenvalue equation can be written as

$$\underset{\{𝐦^{}\}}{}S_{\mathrm{𝐦𝐦}^{}}\varphi _𝐦^{}=\gamma \varphi _𝐦,$$
(58)
or more explicitly

$`{\displaystyle \underset{i=1}{\overset{M}{}}}`$ $`{\displaystyle \frac{m_i+1}{N}}\varphi _{m_1,m_2,\mathrm{},m_i+1,m_{i+1}1,\mathrm{},m_M}`$ (59)
$`=`$ $`\gamma \varphi _{m_1,m_2,\mathrm{},m_i,m_{i+1},\mathrm{},m_M}.`$ (60)
In the above expression we have set $`\varphi _{m_1,m_2,\mathrm{},m_M}=0`$ for those $`|𝐦=|m_1,m_2,\mathrm{},m_M`$ that do not satisfy the constraint (3).

It is not an easy task to diagonalize $`S`$ directly. Thus we adopt another strategy. We first construct a generating function for $`\varphi _{m_1,m_2,\mathrm{},m_M}`$ and then transform the matrix eigenvalue equation (58) to its differential equation form. We find that the differential equation can be solved analytically.

By introducing variables $`x_1,x_2,\mathrm{},x_M`$, the generating function can be defined as

$$f(x_1,x_2,\mathrm{},x_M)\underset{\{𝐦\}}{}\varphi _{m_1,m_2,\mathrm{},m_M}x_1^{m_1}x_2^{m_2}\mathrm{}x_M^{m_M}.$$
(61)
Hereafter we also use the following expression

$$f(𝐗)=\underset{\{𝐦\}}{}\varphi _𝐦X^𝐦,$$
(62)
where $`𝐗`$ and $`X^𝐦`$ are defined by

$$𝐗\left[\begin{array}{ccccc}x_1& & & & \\ x_2& & & & \\ \mathrm{}& & & & \\ x_M& & & & \end{array}\right],X^𝐦x_1^{m_1}x_2^{m_2}\mathrm{}x_M^{m_M}.$$
(63)
To proceed further, note that $`f(𝐗)`$ satisfies the following two relations:

$`_{x_i}f(𝐗)`$ $`=`$ $`{\displaystyle \underset{\{𝐦\}}{}}(m_i+1)\varphi _{m_1,m_2,\mathrm{},m_i+1,\mathrm{},m_M}X^𝐦,`$ (64)
$`x_if(𝐗)`$ $`=`$ $`{\displaystyle \underset{\{𝐦\}}{}}\varphi _{m_1,m_2,\mathrm{},m_i1,\mathrm{},m_M}X^𝐦,`$ (65)
as can be easily checked. Multiplying $`X^𝐦`$ on both sides of Eq.(60), summing over all $`\{𝐦\}`$, and using the results of (65) , we get

$$\underset{i=1}{\overset{M}{}}\frac{x_{i+1}}{N}_{x_i}f(𝐗)=\gamma f(𝐗),$$
(66)
or equivalently

$$\underset{i=1}{\overset{M}{}}x_{i+1}_{x_i}\mathrm{ln}[f(𝐗)]=N\gamma ,$$
(67)
which is the desired differential equation form of the eigenvalue equation (60). Define

$$x_{q_j}=x_1q_j+x_2q_j^2+\mathrm{}+x_Mq_j^M,$$
(68)
we find

$$\underset{i=1}{\overset{M}{}}x_{i+1}_{x_i}\mathrm{ln}(x_{q_j})=q_j^1=q_j^{}.$$
(69)
This implies that the complete solution of $`\mathrm{ln}[f(𝐗)]`$ can be written as

$$\mathrm{ln}[f_𝐧(𝐗)]=\underset{j=1}{\overset{M}{}}n_j\mathrm{ln}(x_{q_j}),$$
(70)
which gives us

$$f_𝐧(𝐗)=\underset{j=1}{\overset{M}{}}x_{q_j}^{n_j}X_𝐪^𝐧.$$
(71)
Here $`f_𝐧(𝐗)`$(a homogeneous $`N`$th power function) and the eigenvalue $`\gamma _𝐧`$ are characterized by $`𝐧=[n_1,n_2,\mathrm{},n_M]`$ and $`𝐪=[q_1,q_2,\mathrm{},q_M]`$, satisfying

$$N=\underset{j=1}{\overset{M}{}}n_j,\gamma _𝐧=\frac{1}{N}\underset{j=1}{\overset{M}{}}n_jq_j^{}=\frac{𝐧𝐪^{}}{N}.$$
(72)
Denoting the $`𝐧`$-th eigenvector of $`S`$ as $`\varphi (𝐧)`$, Eq. (58) and (62) now become

$$\underset{\{𝐦^{}\}}{}S_{\mathrm{𝐦𝐦}^{}}\varphi _𝐦^{}(𝐧)=\gamma _𝐧\varphi _𝐦(𝐧)$$
(73)
and

$$f_𝐧(𝐗)=\underset{\{𝐦\}}{}\varphi _𝐦(𝐧)X^𝐦=X_𝐪^𝐧.$$
(74)
To diagonalize $`S`$ we first define an orthogonal transformation matrix $`U`$:

$$U_{\mathrm{𝐦𝐧}}=\varphi _𝐦(𝐧),$$
(75)
where $`\varphi _𝐦(𝐧)`$ according to Eq. (74) is the coefficient of $`X^𝐦`$ that appears in the expansion of $`f_𝐧(𝐗)=X_𝐪^𝐧`$.

We now are ready to solve the matrix $`U^1`$. Multiplying $`X^𝐦=x_1^{m_1}x_2^{m_2}\mathrm{}x_M^{m_M}`$ on both sides of

$$\underset{\{𝐧\}}{}\varphi _𝐦(𝐧)U_{\mathrm{𝐧𝐥}}^1=\underset{\{𝐧\}}{}U_{\mathrm{𝐦𝐧}}U_{\mathrm{𝐧𝐥}}^1=\delta _{\mathrm{𝐦𝐥}}$$
(76)
and summing over all possible $`\{𝐦\}`$, we get

$$\underset{\{𝐧\}}{}f_𝐧(𝐗)U_{\mathrm{𝐧𝐦}}^1=\underset{\{𝐧\}}{}X_𝐪^𝐧U_{\mathrm{𝐧𝐦}}^1=X^𝐦.$$
(77)
Furthermore, define two vectors $`𝐘`$ and $`𝐗_𝐪`$ as

$`𝐘`$ $`=`$ $`\left[\begin{array}{cccc}y_1& & & \\ y_2& & & \\ \mathrm{}& & & \\ y_M& & & \end{array}\right]\left[\begin{array}{cccc}q_1& q_2& \mathrm{}& q_m\\ q_1^2& q_2^2& \mathrm{}& q_m^2\\ \mathrm{}& \mathrm{}& \mathrm{}& \mathrm{}\\ q_1^M& q_2^M& \mathrm{}& q_m^M\end{array}\right]\left[\begin{array}{cccc}x_1& & & \\ x_2& & & \\ \mathrm{}& & & \\ x_M& & & \end{array}\right],`$ (90)
and

$$𝐗_𝐪\left[\begin{array}{cccc}x_{q_1}& & & \\ x_{q_2}& & & \\ \mathrm{}& & & \\ x_{q_M}& & & \end{array}\right],$$
(91)
we have

$$𝐘=\sqrt{M}R𝐗.$$
(92)
Using the same notation and remember that $`R^1=R^{}`$, we find

$$𝐗=\frac{1}{\sqrt{M}}R^{}𝐘=\frac{1}{M}𝐘_𝐪^{}.$$
(93)
These results further lead to

$`X^𝐦`$ $`=`$ $`{\displaystyle \frac{1}{M^N}}(Y_𝐪^{})^𝐦={\displaystyle \frac{1}{M^N}}Y_𝐪^{\stackrel{~}{𝐦}}={\displaystyle \frac{1}{M^N}}f_{\stackrel{~}{𝐦}}(𝐘),`$ (94)
$`=`$ $`{\displaystyle \frac{1}{M^N}}{\displaystyle \underset{\{𝐧\}}{}}\varphi _𝐧\left(\stackrel{~}{𝐦}\right)Y^𝐧={\displaystyle \frac{1}{M^N}}{\displaystyle \underset{\{𝐧\}}{}}\varphi _𝐧\left(\stackrel{~}{𝐦}\right)X_𝐪^𝐧`$ (95)
$`=`$ $`{\displaystyle \frac{1}{M^N}}{\displaystyle \underset{\{𝐧\}}{}}f_𝐧(𝐗)\varphi _𝐧\left(\stackrel{~}{𝐦}\right).`$ (96)
Comparing Eq. (96) with (77), we get:

$$U_{\mathrm{𝐧𝐦}}^1=\frac{1}{M^N}\varphi _𝐧\left(\stackrel{~}{𝐦}\right),$$
(97)
where we have used the relations:

$`Y_𝐪^{}^𝐦`$ $`=`$ $`y_{q_1^{}}^{m_1}y_{q_2^{}}^{m_2}y_{q_3^{}}^{m_3}\mathrm{}y_{q_{M1}^{}}^{m_{M1}}y_{q_M^{}}^{m_M},`$ (98)
$`=`$ $`y_{q_{M1}}^{m_1}y_{q_{M2}}^{m_2}y_{q_{M3}}^{m_3}\mathrm{}y_{q_1}^{m_{M1}}y_{q_M}^{m_M},`$ (99)
$`=`$ $`y_{q_1}^{m_{M1}}y_{q_2}^{m_{M2}}y_{q_3}^{m_{M3}}\mathrm{}y_{q_{M1}}^{m_1}y_{q_M}^{m_M},`$ (100)
$`=`$ $`Y_𝐪^{\stackrel{~}{𝐦}}=f_{\stackrel{~}{𝐦}}(𝐘),`$ (101)
and defined $`\stackrel{~}{𝐦}`$ as

$$\stackrel{~}{𝐦}[m_{M1},m_{M2},m_{M3},\mathrm{},m_1,m_M].$$
(102)
Finally we obtain the desired solution of $`𝐦|P^s|𝐦_0`$:

$`𝐦|P^s|𝐦_0`$ $`=`$ $`(S^s)_{\mathrm{𝐦𝐦}_0},`$ (103)
$`=`$ $`(U\mathrm{\Gamma }^sU^1)_{\mathrm{𝐦𝐦}_0},`$ (104)
$`=`$ $`{\displaystyle \frac{1}{M^N}}{\displaystyle \underset{𝐦^{}}{}}\gamma _𝐦^{}^s\varphi _𝐦(𝐦^{})\varphi _𝐦^{}(\stackrel{~}{𝐦}_0),`$ (105)
where $`\mathrm{\Gamma }`$ is the eigenvalue matrix of $`S`$, which has components $`\mathrm{\Gamma }_{\mathrm{𝐦𝐦}^{}}=\gamma _𝐦\delta _{\mathrm{𝐦𝐦}^{}}`$.

## IV poincaré cycle

In this section, we study the Poincaré cycle of our periodic urn model. For simplicity we first consider the situation that initially all the $`N`$ balls are stayed in the last urn, i.e.,

$$𝐦_0=[0,0,0,\mathrm{},N].$$
(106)
From (106) and (102) we have

$$𝐦_0=\stackrel{~}{𝐦}_0.$$
(107)
Now we want to know how many time steps on average are required for all of the $`N`$ balls to return to the last urn (the initial state). We thus have to calculate

$$𝐦_0|P^s|𝐦_0=\frac{1}{M^N}\underset{𝐦}{}\gamma _𝐦^s\varphi _{𝐦_0}(𝐦)\varphi _𝐦(\stackrel{~}{𝐦}_0).$$
(108)
Recall that $`\varphi _𝐦(𝐧)`$ is nothing but the coefficient of $`X^𝐦=x_1^{m_1}\mathrm{}x_M^{m_M}`$ appearing in the expansion of $`f_𝐧(𝐗)=X_𝐪^𝐧=x_{q_1}^{n_1}\mathrm{}x_{q_M}^{n_M}`$. From Eq. (74) and (106), we have

$$\varphi _{𝐦_0}(𝐦)=1,$$
(109)
$$\varphi _𝐦(\stackrel{~}{𝐦}_0)=\frac{N!}{m_1!m_2!m_3!\mathrm{}m_M!}\left(\begin{array}{cc}N& \\ 𝐦& \end{array}\right),$$
(110)
and hence

$$𝐦_0|P^s|𝐦_0=\frac{1}{M^N}\underset{𝐦}{}\left(\begin{array}{cc}N& \\ 𝐦& \end{array}\right)\gamma _𝐦^s𝒫(s).$$
(111)
Here $`𝒫(s)`$ represents the transition probability for the system to return to the initial state after $`s`$ steps. It does not preclude the possibility that the initial state has already been re-arrived before.

Since Poincaré cycle is defined as the time interval required for the event of first return to happen, so we have to do more calculations to extract what we really want. We define a function $`𝒬(s)`$ as the probability for the event of first return to happen at the $`s`$th step. The Poincaré cycle can thus be defined as

$$\mathrm{P}.\mathrm{C}.=\underset{s=0}{\overset{\mathrm{}}{}}s𝒬(s).$$
(112)
By definition $`𝒬(s)`$ relates to $`𝒫(s)`$ via the relation

$$𝒫(s)=𝒬(s)+\underset{k=1}{\overset{s1}{}}𝒬(k)𝒫(sk),$$
(113)
and hence $`𝒬(s)`$ can be calculated from $`𝒫(s)`$. To ease the calculation we now use again the generating function method. We first define two generating functions:

$$h(z)\underset{s=1}{\overset{\mathrm{}}{}}𝒫(s)z^s,g(z)\underset{s=1}{\overset{\mathrm{}}{}}𝒬(s)z^s,$$
(114)
and then we find from Eq. (113) that

$$g(z)=\frac{h(z)}{h(z)+1}.$$
(115)
These two generating functions also lead to

$$\underset{s=0}{\overset{\mathrm{}}{}}s𝒬(s)=\left(\frac{dg}{dz}\right)_{z=1},g^{}=\frac{h^{}}{(1+h)^2},$$
(116)
which can determine the Poincaré cycle.

We now calculate $`h(z)`$. From Eq. (111) and (129) we obtain

$`h(z)`$ $`=`$ $`{\displaystyle \frac{1}{M^N}}{\displaystyle \underset{𝐦}{}}\left(\begin{array}{cc}N& \\ 𝐦& \end{array}\right){\displaystyle \underset{s=1}{\overset{\mathrm{}}{}}}\left(\gamma _𝐦z\right)^s,`$ (119)
$`=`$ $`{\displaystyle \frac{1}{M^N}}{\displaystyle \underset{𝐦}{}}\left(\begin{array}{cc}N& \\ 𝐦& \end{array}\right)\left({\displaystyle \frac{\gamma _𝐦z}{1\gamma _𝐦z}}\right).`$ (122)
Since we know from Eq. (72) that $`\gamma _{𝐦_0}=1`$, thus when $`z1^{}`$ , $`h(z)`$ becomes singular:

$$\underset{z1^{}}{lim}h(z)=\frac{1}{M^N}\frac{z}{1z}+\text{regular function}.$$
(123)
In this limit, we obtain

$$\underset{z1^{}}{lim}g^{}=\underset{z1^{}}{lim}\frac{h^{}}{(1+h)^2}=M^N,$$
(124)
which gives us the desired Poincaré cycle:

$$\mathrm{P}.\mathrm{C}.=\underset{s=0}{\overset{\mathrm{}}{}}s𝒬(s)=M^N.$$
(125)
To understand this result, we refer to the ergodic theorem , which says that if one waits a sufficiently long time, the locus of the representative point of a system will cover the entire accessible phase space. For our periodic urn model, the “representative point” corresponds to the microstate of the arrangement of balls, the “accessible phase space” is the set of total microstates, and the “locus” means the evolution history of the system (See Fig. 1). Since in our periodic model every microstate has the same probability to appear (a fundamental assumption of statistical mechanics), and the set of total microstates contains $`M^N`$ microstates, therefore on average one has to wait $`M^N`$ time steps to see a microstate reappearing again.

What will be the Poincaré cycle if the initial state is different from the $`𝐦_0`$ in Eq. (106)? Let’s denote the initial state by $`𝐝`$:

$$|𝐝=|d_1,d_2,\mathrm{},d_M.$$
(126)
Now we have

$$𝒫(s)=𝐝|P^s|𝐝=\frac{1}{M^N}\underset{\{𝐦\}}{}\gamma _𝐦^s\varphi _𝐝(𝐦)\varphi _𝐦(\stackrel{~}{𝐝}),$$
(127)
and

$$h(z)=\frac{1}{M^N}\underset{\{𝐦\}}{}\varphi _𝐝(𝐦)\varphi _𝐦(\stackrel{~}{𝐝})\left(\frac{\lambda _𝐦z}{1\lambda _𝐦z}\right).$$
(128)
In general it is difficult to calculate $`\varphi _𝐝(𝐦)`$ and $`\varphi _𝐦(\stackrel{~}{𝐝})`$. However, we do not need to calculate them all. Remember that to determine the Poincaré cycle the the knowledge of the asymptotic form of $`h(z)`$ near $`z=1`$ is enough. This is given by

$$\underset{z1^{}}{lim}h(z)=\frac{\varphi _𝐝(𝐦_0)\varphi _{𝐦_0}(\stackrel{~}{𝐝})}{M^N}\frac{z}{1z}+\mathrm{regular}\mathrm{function}.$$
(129)
Substituting

$$\varphi _𝐝(𝐦_0)=\left(\begin{array}{cc}N& \\ 𝐝& \end{array}\right),\varphi _{𝐦_0}(\stackrel{~}{𝐝})=1$$
(130)
into (129), we find

$$\mathrm{P}.\mathrm{C}.=\underset{z1^{}}{lim}\frac{h^{}}{(1+h)^2}=\frac{M^N}{\left(\begin{array}{cc}N& \\ 𝐝& \end{array}\right)}$$
(131)
This result conforms with our intuition. Consider the difference between the concepts of configuration and microstate , one can easily understand this result because the degeneracy of the configuration $`𝐝`$ is given by

$`\left(\begin{array}{cc}N& \\ 𝐝& \end{array}\right)={\displaystyle \frac{N!}{d_1!d_2!\mathrm{}d_M!}}.`$

## V summary

In this work we propose a generalized Ehrenfest urn model of many urns arranged periodically along a circle. We solve an $`N`$-ball, $`M`$-urn problem explicitly. The evolution of the system is studied, and the average number of balls in a certain urn at any time step is calculated. We find that this mean value oscillates several times before it arrives the stationary value. We also obtained the Poincaré cycle. The result is consistent with the consideration of the total number of possible microstates of the system.

## Acknowledgments

The authors would like to acknowledge helpful discussions with professor T. F. Jiang. This work received support from the National Science Council, Republic of China through Grant No. NSC 90-2811-M-009-019.


text:
# Angular Distribution of Gamma-ray Bursts and Weak Lensing

## 1 Introduction

Direct identification of X-ray, optical and radio counterparts of long duration, $`t>2`$ sec, Gamma-Ray Bursts (GRB), and hence their host galaxies has recently resolved the GRB distance scale controversy: the observed redshifts span a wide range, from $`z0.1`$ to $`z5`$ (Reichart & Lamb, 2001). This suggests a number of uses of GRBs as cosmologically distributed probes. Since they are believed to be associated with compact remnants of massive stars, it has been suggested they be used to trace star formation rate obscured by dust (Totani, 1999; Djorgovski et al., 2001), star formation rate at very high redshifts (Lamb & Reichart, 2000), as probes of the metal enrichment of the interstellar medium (Fiore, 2001), intergalactic medium (Fiore et al., 2000), and galactic and intergalactic dust at high redshifts (Perna & Aguirre, 2000).

Here, we use GRBs as sources for weak lensing by the large scale structure; our ultimate goal is to probe the mass distribution on $`>\mathrm{\hspace{0.17em}10}h^1`$ Mpc scales, at a typical redshift of 0.3.

Unlike weak shear lensing, which is detected through shape distortion of lensed resolved galaxies (Mellier, 1999; Bartelmann & Schneider, 2001), weak lensing of unresolved point sources, like GRBs, would manifest itself through the angular (anti-)correlation between sources and lenses, the latter being the intervening clumpy mass distribution, which is assumed to be traced by visible galaxies (Fugmann, 1990; Rodrigues-Williams & Hogan, 1994; Bartelmann, 1995). The sign of the correlations is determined by the slope of the source number counts in the appropriate redshift interval. For example, bright QSOs have steep number counts implying that positive correlations between these and intervening galaxies should be expected. In fact, optically selected LBQS QSOs, with $`z1.0`$, and radio selected 1Jy QSOs, with $`0.5z<\mathrm{\hspace{0.17em}2.5}`$, are independently correlated with faint APM galaxies on angular scales of $`1^{}`$ (Williams & Irwin, 1998; Norman & Williams, 2000). The linear extent of the structures is $`>\mathrm{\hspace{0.17em}15}h^1`$ Mpc at the redshift of typical lenses, $`z0.10.4`$. Radio selected 1Jy sources are also correlated with IRAS galaxies (Bartelmann & Schneider, 1994; Bartsch et al., 1997), and Zwicky clusters (Seitz & Schneider, 1995) on angular scales of $`12^{}`$.

In these studies, the redshifts of the sources and probable redshift distribution of the galaxies do not overlap, insuring that physical associations do not contaminate the lensing signal. Qualitative signature of these correlations are those of lensing, however, it is hard to explain the results quantitatively: the observed correlations persist to scales of $`>\mathrm{\hspace{0.17em}1}^{}`$ with amplitude of $`\omega (1^{})0.02`$, whereas $`\omega (1^{})0.002`$ is expected if we live in a Universe with $`\mathrm{\Omega }_m0.3`$, and matter power spectrum not too different from the observed galaxy power spectrum on large scales (Dolag & Bartelmann, 1997; Sanz et al., 1997; Bartsch et al., 1997). If lensing induced correlations were indeed as small as predicted, they could not have been detected on $`1^{}`$ scales given the available number of sources (Bartelmann & Schneider, 1993; Bartsch et al., 1997).

These results by themselves are puzzling. However, they become problematic in view of the recent observations of cosmic weak shear lensing, which agree well with the currently accepted cosmological model, matter fluctuation spectrum, and biasing parameter close to unity. There are now several independent determinations of weak shear lensing in fields of up to $`30^{}`$ diameter (Bacon et al., 2000; Hoekstra et al., 2002; Kaiser et al., 2000; Maoli et al., 2001; Van Waerbeke et al., 2000; Wittman et al., 2000), and the agreement between them is impressive. Furthermore, when cast in terms of the best estimate for $`\sigma _8\mathrm{\Omega }_{m}^{}{}_{}{}^{0.5}`$, weak lensing results agree remarkably well with cluster normalization constraints (Mellier et al., 2001).

GRBs provide us with a different set of cosmologically distributed sources, which should also be affected by weak lensing, just as QSOs are. However, as we explain below, unlike QSOs, which are predicted and shown to exhibit correlations with galaxies, GRBs, due to their distribution in fluence and angular positional error, are predicted to be anti-correlated with the foreground lensing mass distribution. Thus, GRBs should provide a new and interesting test of weak lensing on large scales.

There are two potential difficulties in using GRBs in a weak lensing study.

First, GRB positions are not well localized, with errors ranging from a fraction of a degree to as high as $`30^{}`$ for some events with low fluences. Obviously, correlation scales that can be reliably probed are limited by the GRB position errors on the sky. To circumvent this problem we use GRB subsets with upper limits on error.

Second, since individual GRB redshifts are largely unknown, some of the GRBs at low redshifts will be physically associated with the galaxies which we use to trace the lensing mass, and will ‘contaminate’ the lensing induced signal. Currently, only about 26 GRBs, i.e. a very small fraction of all events have confirmed redshifts111http://www.aip.de/~jcg/grbgen.html, site maintained by J. Greiner, Astrophysical Institute Potsdam.. Of these, about 3 are at $`z<0.4`$. Guided by the known redshifts, several workers have proposed empirical relations connecting observable properties of GRBs, such as functions constructed out of time information, to their luminosities (Stern et al., 1999; Norris et al., 2000; Reichart et al., 2001). Redshift distribution implied by these indicate that the true redshift distribution is very broad, extending to $`z510`$. The fraction of low redshift events is expected to be small; for example, using the model of Reichart & Lamb (2001) and their Fig. 3 histogram, the fraction of events below $`z_s=0.2(0.4)`$ is roughly estimated to be 6% (13%). As tracers of dark matter on large scales we use APM galaxies and Abell-Corwin-Olowin (ACO) clusters, which peak at $`z0.2`$, and extend to a maximum $`z=0.4`$. So the fraction of GRBs arising in these structures should be small, $`<10\%`$. Furthermore, and more importantly, we predict and find anti-correlations between GRBs and intervening galaxies, and so physical associations, if any, would have diminished the amplitude of lensing induced anti-correlations. Thus our observed signal (Section 4) is a lower limit. As we show in Section 4.3, the observed amplitude of lensing induced anti-correlations is hard to explain within standard cosmological models, so at this point there is no motivation to carefully subtract the effect of physical associations from the observed anti-correlation signal; we neglect the effect of physical associations in our analysis.

The plan of the paper is as follows. After describing the data in Section 2 we use the fluence and error information of Current BATSE Catalog GRBs to predict the amplitude and sign of GRB-galaxy correlations, i.e. compute magnification bias (Section 3). In Sections 4 and 5 we estimate the strength of observed correlations between GRBs and APM galaxies and ACO clusters, respectively. Anti-correlations are found in both cases, as predicted, but the amplitude of the effect is higher than expected. In Sections 6 we discuss the various possibilities for reconciling observations with theory, and propose a speculative scenario. We summarize and discuss our findings in Section 7.

## 2 Data selection

### 2.1 Gamma-Ray Bursts

We use the Current BATSE GRB Catalog, ending with trigger number 8121, which occurred on May 26, 2000. From these events we select those that were not overwritten by a later more intense trigger, and that have non-zero fluences in the 50-100 keV and 100-300 keV energy channels. GRB fluxes and fluences are recorded in 4 channels, which cover energy ranges 20-50 keV, 50-100 keV, 100-300 keV, and $`>300`$ keV, respectively. We use the fluences in the middle two channels because the corresponding energy range has the peak flux, and coincides with the energy range of the nominal BATSE on-board burst trigger (http://www.batse.msfc.nasa.gov/batse/grb/catalog/4b/4br\_flux.html). These cuts leave us with 2038 events.

BATSE events are not well localized in the sky; positional errors come in two flavors, statistical and systematic. Statistical errors, $`e`$, are recorded in the BATSE Catalog, and range from a fraction of a degree to $`30^{}`$ degrees. (All quoted errors are $`\sqrt{e_x^2+e_y^2}`$.) The peak of the statistical error distribution is at about $`3^{}`$. Systematic errors, $`e_{sys}`$, were analyzed for the revised 4B Catalog (Paciesas et al., 1999), and were found to have a modified Gaussian distribution such that 78% of the events have errors of $`1.85^{}`$. For our purposes an approximation used for the 3B Catalog will suffice: we assume the systematic error distribution to be a Gaussian such that 68% of the events have systematic errors of $`1.6^{}`$. Assuming that statistical and systematic errors are independent, the total error is $`e_{tot}=\sqrt{e_{sys}^2+e^2}`$. The statistical errors and fluences of the 2038 bursts are plotted in Fig. 1. Solid squares represent 10 BATSE GRBs with known fluences and redshifts.

BATSE Catalogs have a declination dependent completeness level on the sky, with a 74% difference between the most ($`\delta =75^{}`$) and least ($`\delta =10^{}`$) completely covered declinations (Paciesas et al., 1999). This incompleteness will not matter for our analysis in Section 4, because we compare galaxy density around each GRB with the expected average density on the same APM plate. This incompleteness will be important in the cross-correlation with ACO galaxy clusters, Section 5, and will be accounted for.

### 2.2 APM Galaxy Catalog

APM Galaxy Catalog (Irwin et al., 1994) is digitized POSSI plate material, with magnitude limits 20.0 and 21.5 on ‘red’ and ‘blue’ POSSI plates in the Northern Hemisphere, and 21.0 and 22.5 on ‘red’ and ‘blue’ UKST plates in the Southern Hemisphere, with internal accuracy of $`0.1`$ mag for all but the faintest objects, and external, i.e. plate-to-plate accuracy of $`0.3`$ mag. These are tolerable errors for our purposes: internal accuracy is sufficiently good to trace the light fluctuations on each plate, and external accuracy is less important to us because we confine analysis for each GRB to a single APM plate. Similar to previous work (Williams & Irwin, 1998) we select galaxies as objects that are classified as extended on the ‘red’ plates by the APM. The magnitude range is set to probe the most distant regions reached by the APM: on POSSI, $`18.5m_R20.0`$, and on UKST, $`19.5m_R21.0`$.

In the Northern Hemisphere APM plates do not cover a band between $`20^{}`$ and $`+20^{}`$ Galactic Latitude, and in the Southern Hemisphere the coverage is even less complete. Furthermore, even though each plate is $`5.8^{}\times 5.8^{}`$, because of plate vignetting we do not use GRBs located beyond $`2.5^{}`$ from the plate center. Thus, we only use a subset of 732 GRBs in the APM analysis in Section 4.

## 3 Predicting the effect of Weak Lensing on the GRBs

The intervening mass distribution between us and GRB events is an uneven lens that stretches some areas of the background sky and shrinks others. The redshift distribution of the observed APM galaxies is such that they trace dark matter from $`z_{min}0.1`$ to a maximum of $`z_{max}0.4`$. (see Fig. 1 of Williams & Irwin 1998). Let the total average lensing optical depth of this slab of matter be $`\kappa _0`$:

$$\kappa _0=\rho _{crit}\mathrm{\Omega }_0_{z_{min}}^{z_{max}}\frac{(cdt/dz)(1+z)^3}{\mathrm{\Sigma }_{crit}(z,z_s)}𝑑z.$$
(1)
Here, $`cdt`$ is the thickness of the lensing slice at redshift $`z`$, $`\rho _{crit}\mathrm{\Omega }_0(1+z)^3`$ is its mass density, and $`\mathrm{\Sigma }_{crit}(z,z_s)`$ is the critical lensing surface mass density at $`z`$ for a source at $`z_s`$. Within this redshift range, the fluctuations in the projected surface mass density are $`\delta \sigma /\sigma `$, and are assumed to be small on large scales. For a specified cosmology and source and lens redshifts, these fluctuations are translated into (de-)magnifications on the sky, with respect to the ‘smooth Universe’ case: $`M=(1\kappa )^21+2\kappa `$, where $`\kappa =\kappa _0(\delta \sigma /\sigma )`$.

For sources with exactly known positions, like QSOs, this magnification field $`M`$, combined with the QSO number counts for the relevant range of redshifts, can be translated into an observed distribution of number density on the sky of sources down to some specified flux limit. If $`\alpha =d\mathrm{log}n_{QSO}(<m)/dm`$ is the slope of the number counts near the survey flux limit, then the number density of QSOs will be a factor of $`q=M^{2.5\alpha 1}`$ different from the ‘smooth Universe’ ($`M=1`$) case; $`q`$ is called source over- or underdensity, depending on whether it is greater or less than 1. The fact that it is different from 1 is called the magnification bias. Bright QSOs, whose number counts are steep, would appear to be $`q`$ times more abundant in the directions of mass concentrations, i.e. they would appear to be correlated on the sky with the nearby structure. Faint QSOs, on the other hand, have a shallow number counts slope, and so will be anti-correlated with the foreground lenses.

With GRBs, the calculation of $`q`$ is somewhat different, because in addition to limiting these in fluence $`f`$, one also needs to place an upper limit on GRB positional error $`e`$, so that the GRBs have a reasonable chance of being within the specified area of correlations. Even though positional errors are not directly affected by lensing, the observed GRB errors are well correlated with their fluences, therefore, as a GRB’s fluence is changed through magnification, so is its error. Thus, in order to correctly predict $`q`$ one needs to consider GRB distribution in $`e`$ and $`f`$, Fig. 1.

Let $`p(e|f)`$ be the normalized probability distribution of GRB sources in statistical positional error $`e`$ at a given fluence $`f`$. As $`f`$ we use the sum of fluences in Channels 2 and 3. Probability distribution $`p(e|f)`$ is estimated from the BATSE data itself, by binning the distribution presented in Fig. 1 by fluence and error, and then normalizing distributions in $`e`$ for each fluence separately. If the sources are limited in error by $`e_1`$ and $`e_2`$ from below and above, and in fluence by $`f_1`$ from below, then the number of such sources seen behind a smooth patch of lens with magnification $`M`$ is given by,

$$n_{GRB}(e_1,e_2,f_1,M)=\frac{1}{M}_{f_1}^{\mathrm{}}\left[n_{GRB,0}(f^{}/M)_{e_1}^{e_2}p(e|f^{})𝑑e\right]𝑑f^{},$$
(2)
where $`n_{GRB,0}(f)`$ is the distribution of sources in fluence one would see in a Universe with completely smooth mass distribution. This distribution is not observable; however, given the small typical magnifications it is reasonable to assume that the observed sources give a fair representation of $`n_{GRB,0}(f)`$.

The ratio of the number of sources observed and the number that would be seen if the mass were smoothly distributed everywhere, is the over- or underdensity,

$$q(e_1,e_2,f_1,M)=\frac{n_{GRB}(e_1,e_2,f_1,M)}{n_{GRB}(e_1,e_2,f_1,M=1)}$$
(3)
Figure 2 shows the result, for a few combinations of $`e_1`$, $`e_2`$, and $`f_1`$. The upper two lines (cross and plus symbols) are for GRB subsets limited in fluence only, with no restrictions imposed on error. The lower three lines are for GRB subsets limited in error, with no limits placed on fluence. We used all 2038 GRBs in the Current BATSE Catalog to construct the $`(M,q)`$ relation in Fig. 2. Had we used, say the 732 GRBs that are found on APM plates the shape of the curves above would have been the same, within the noise. Similarly, the plot is not very sensitive to the particular choice of fluence channel. Note that all lines go through $`(M,q)=(1,1)`$, as they should. Some of the $`q(M)`$ lines vary non-monotonically with magnification, and most lines vary erratically. This is due to the noise associated with the finite number of GRB points in Fig. 1. Shot noise is especially pronounced when GRBs with small errors and/or large fluences are considered, and when $`M<1`$: in that case the $`n_{GRB,0}`$ term in eq. 2 refers to a very small number of GRBs at highest fluences, and so shot noise fluctuations can make the resulting $`q`$ vary a lot.

Most of the subsets that can be constructed from the total GRB catalog using different fluence/error cuts are predicted to be anti-correlated with lenses. In general, overdensities would be expected only if there is a large ‘reservoir’ of sources just below the detection limit. The fact that the predicted GRB $`q`$’s tend to be $`<1`$ can be seen directly from Fig. 1: for almost any $`e`$, $`f`$ cut there is no large ‘reservoir’ of sources just below the $`e`$, $`f`$ limits ready to be magnified into the observed subset. The only subsets of GRBs which are predicted to be correlated with lenses are those with large fluences, $`f_1>2.5\times 10^5`$ erg/cm2. However, these subsets contain very small numbers of GRBs, so the uncertainties in the $`(M,q)`$ relation for these are probably large.

Dashed lines in Fig. 2 are fits to the $`(e_1,e_2,f_1)`$=$`(0.0^{},1.0^{},0.0)`$ and $`(0.5^{},1.0^{},0.0)`$ GRB subsets, which we will be considering in some detail in the next Section. We wanted fits of the form $`q=M^\beta `$, to mimic the QSO relation, $`q=M^{2.5\alpha 1}`$. In the present case, $`\beta =0.62`$ and $`0.45`$ provide adequate fits for the two subsets respectively; the corresponding $`\alpha `$’s are $`0.152`$, and $`0.22`$. The fits are only rough, but are completely adequate for our purposes.

Having made predictions as to the amplitude and sign of the GRB correlations with the foreground matter, we can now proceed to do the corresponding ‘observations’.

## 4 GRB–APM galaxy correlations

In principle, we want to determine the number density of GRBs as a function of projected mass excess, $`\delta \sigma /\sigma `$. However, since GRBs are rare and galaxies are plentiful, we instead estimate galaxy density in circles around GRBs. Our analysis later, Section 4.3 will take this difference into account.

Let the number of galaxies within $`\theta `$ degrees of a GRB, and in the magnitude ranges specified in Section 2.2, be $`n_{gal,D}`$. The latter should be normalized by average expected galaxy count within similarly sized random $`\theta `$-patches on the sky. In our analysis, because the APM plate-to-plate magnitude calibration can be rather uncertain, about 0.3 mag, we select control areas for normalization from the same plate as the corresponding GRB. Furthermore, because the object density on plates varies as a function of distance from center, the control patches are restricted to lie at the same distance from plate center as the GRB, $`d_{cen}`$. Fig. 3 illustrates our selection of control patches. For each GRB, $`N`$ such random positions are created. We scale $`N`$ linearly with $`d_{cen}`$, such that at $`d_{cen}=1^{}`$, $`N=100`$. (We tried using same $`N`$ for all GRBs; the results below did not change substantially.) The number of galaxies in $`\theta `$ circles around these are also recorded, and the average is calculated for every GRB, $`n_{gal,R}`$. The galaxy excess, $`[n_{gal,D}/n_{gal,R}]1`$, is then equal to $`b(\delta \sigma /\sigma )`$, where $`b`$ is the biasing parameter of APM galaxies with respect to the underlying mass.

### 4.1 Two GRB subsets most likely to show weak lensing effects

Using Fig. 2 we select GRB subsets that are most likely to show weak lensing induced signature. We choose subsets with $`(e_1,e_2)=(0.5^{},1.0^{})`$ and $`(0^{},1.0^{})`$. How do we decide on the size of the $`\theta `$-patch around each GRB? Obviously, each GRB should have a good chance of actually being located within these circles. The limiting factor here is the systematic positional error: even if a GRB has zero statistical error there is only a 25% chance that it will be actually within a $`0.5^{}`$ radius drawn around it. On the other hand, the size of the APM plate limits the size of $`\theta `$ from above; we settle on $`\theta =1.5^{}`$.

Given that the GRB and the $`N`$ control $`\theta `$-patches are equidistant from plate center, to avoid significant overlap between these we restrict our GRB subset further by selecting only those at $`d_{cen}\theta `$. Thus GRBs are limited to lie at $`1.5^{}d_{cen}2.5^{}`$. The $`(e_1,e_2)=(0.5^{},1.0^{})`$ and $`(0^{},1.0^{})`$ subsets contain 46 GRBs and 74 GRBs respectively, and we will refer to them by quoting these numbers.

We now ask if the galaxy density around GRBs in the two subsets is different from what one would expect if random points were used in place of GRBs. Instead of generating a set of random points for this purpose, we use the whole set of GRBs found at $`1.5^{}d_{cen}2.5^{}`$, regardless of positional error. There are 448 of these. Since most of these have rather large errors, their real positions are quite far from BATSE recorded positions, so in effect we have a set of randomly selected points.

In Fig. 4 the heavy solid histogram is the distribution of $`n_{gal,D}/n_{gal,R}`$, or, equivalently, $`b(\delta \sigma /\sigma )+1`$ for the 448 ‘random’ points. The dashed and dotted histograms are the 46 and 74 GRB subsets respectively. The averages of the three distributions are 0.998, 0.883, and 0.915, so GRBs in the two subsets are found in the directions of foreground regions that are, on average, 12% and 8% underdense in galaxies. In other words, we detect anti-correlations between GRBs and intervening galaxies, which is what was predicted in Section 3 and Fig. 2. We leave quantitative comparison with predictions until Section 4.3.

Let us now evaluate the statistical significance of this result. For each GRB we calculate $`N_>/N`$, where $`N_>`$ is the number of random $`\theta `$-patches, out of total $`N`$, that have less galaxies in them than the $`\theta `$patch around the real GRB. In other words, $`N_>/N`$ is the rank of the real GRB patch among its ‘random peers’; $`N_>/N=0.5`$ if GRBs are randomly distributed with respect to the foreground galaxies, but if GRBs have an excess of galaxies in their foregrounds then $`N_>/N>0.5`$. If GRBs occupy random positions with respect to the foreground galaxies, then the distribution of $`N_>/N`$ values is known—the cumulative distribution should be linear, and in fact for the whole set of 448 GRBs it is, see solid line in Fig. 5. The dashed and dotted histograms are for the 46 and 74 GRB subsets, respectively; their $`N_>/N`$ values are 0.381 and 0.423, implying that GRBs have a deficit of galaxies in their foregrounds. Using the Kolmogorov-Smirnof (KS) test these two distributions differ from the whole set of 448 GRBs at the 99.68% and 97.62% significance level. Note that the 46 subset \[$`(e_1,e_2)=(0.5^{},1.0^{})`$\] was predicted to be more strongly anti-correlated with galaxies than the 74 subset \[$`(e_1,e_2)=(0.5^{},1.0^{})`$\], which is what is observed, in spite of the 46 subset having less GRBs than the 74 subset.

Since the 46 and 74 GRBs in the two subsets are a part of the 448 ‘random’ points, these results could be a conservative estimate of statistical significance. However, if we use truly random points, the statistical significance of the results is not much higher.

### 4.2 Many GRB subsets

A further test of the statistical significance of the results for the two subsets in the last Section is provided by carrying out the same analysis on many GRB subsets, selected using a range of error and fluence criteria. Subsets are defined by $`\theta `$, the size of the patch around GRBs, statistical error range, $`e_1e_2`$, and the range of GRB’s $`d_{cen}`$. We use four different $`\mathrm{\Delta }e`$ ranges: $`0.5^{}`$, $`1.0^{}`$, $`1.5^{}`$, and $`2.0^{}`$, and include GRBs with errors up to $`10^{}`$. For each of these we use five different $`d_{cen}`$ ranges: $`0.0^{}2.5^{}`$, $`0.5^{}2.5^{}`$, …., $`2.0^{}2.5^{}`$. We try three values of $`\theta `$: $`0.5^{}`$, $`1.0^{}`$, and $`1.5^{}`$. The total number of GRB subsets is 615; the number of GRBs in these subsets varies from 1 to 269.

Of the 615 total cases considered, Table 1 lists those that deviate from the parent GRB set at more than 97% confidence level, as judged by the KS test. The two cases considered in Section 4.1 are marked with a star. Most of the other cases in Table 1 are related to these two subsets. Even though the subsets are not independent, one should still expect some small fraction of the 615 subsets to be significant at $`>97\%`$, so it is not surprising that we have these. It is interesting, however, that the subsets that do show significant deviations form the parent GRB population are the ones that we argued would be most likely to be affected by lensing. Furthermore, all cases reported in Table 1 are anti-correlations, consistent with weak lensing predictions of Section 3.

### 4.3 Can we account for the GRB-galaxy anti-correlations?

Here we attempt to account for the amplitude of anti-correlations of two GRB subsets considered in Section 4.1. As before, we assume that the whole set of 448 GRBs, found at distances $`1.5^{}d_{cen}2.5^{}`$ from plate centers is a collection of random points on the sky, since most of them have very large position errors. Hence, the solid line histogram of Fig. 4 is a fair representation of the distribution of average projected galaxy densities in $`\theta =1.5^{}`$ circles. This counts-in-cells distribution does not include fluctuation power on scales larger than the size of an APM plate. Because we are dealing with weak lensing regime, we can separate the effects of density fluctuations on different projected scales; our ‘observations’ and analysis do not deal with scales larger than $`2.5^{}`$. Fluctuation power below the scale of correlations, $`<\theta 1^{}1.5^{}`$, for example due to galaxy cluster cores, does not affect the correlations significantly. This is illustrated by Dolag & Bartelmann (1997) and Menard & Bartelmann (2002), who have shown that including the fluctuation power on non-linear scales changes the amplitude of lensing-induced correlations by $`10\%`$, for scales $`>10^{}`$.

We start with the assumption that without lensing GRBs and foreground APM galaxies are randomly distributed with respect to one another, on the sky. We now construct a synthetic lensed GRB subset. Suppose a large number of GRBs with the BATSE-observed distribution of fluence and error properties go off in the direction of a patch with some $`[n_{gal,D}/n_{gal,R}]`$ value, picked at random from the solid line distribution in Fig. 4. Using the various relations in Section 3 and the fit in Fig. 2, we can predict the overdensity of GRBs in that direction:

$$q=\left(1\frac{\kappa _0}{b}\left[\frac{n_{gal,D}}{n_{gal,R}}1\right]\right)^{2\beta }\left(1+2\beta \frac{\kappa _0}{b}\left[\frac{n_{gal,D}}{n_{gal,R}}1\right]\right).$$
(4)
The APM galaxies trace matter extending, at most, to $`z=0.4`$. The corresponding optical depth in a flat $`\mathrm{\Omega }_m=0.3`$ Universe is 0.025 for $`z_s=1`$, and 0.031 for $`z_s=3`$; we take $`\kappa _0=0.028`$, as an average. We set $`b=1`$, as estimated on large scales (Peacock, 2002). We use $`\beta =0.62`$ and $`0.45`$, for the 46 and 74 synthetic GRBs subsets, as determined in Section 3). If the number of GRBs going off is 1, then $`q`$ is the probability that it will be observed. Using the latter definition, a GRB is accepted into the synthetic lensed subset with probability $`q`$; if $`q>1`$, the GRB is accepted into the subset, and an additional one is accepted with probability $`q1`$. We repeat the process until we build up two separate subsets of 46 and 74 GRBs each; we repeat the procedure a 1000 times.

The average $`[n_{gal,D}/n_{gal,R}]`$ in the 46 and 74 synthetic subsets is close to 0.996, barely below what one expects for a randomly selected set of points on the sky, and far from 0.883 and 0.915 which are found for the real subsets. In fact, the fraction of synthetic subsets whose average $`[n_{gal,D}/n_{gal,R}]`$ is less than 0.883 and 0.915, is 0.3% and 0.75% respectively. Hence, this model fails to reproduce observations.

How much do we have to change the parameters in eq. 4 to reach agreement with observations? If $`\kappa _0`$ is increased by a factor of 10 (20), then the respective percentages for the 46 and 74 subsets become about 2% (8%), and 3% (10%), so, roughly speaking, the observed subsets can occur with $`<\mathrm{\hspace{0.17em}10}\%`$ probability. We conclude that observations could be deemed to agree with expectations if $`\kappa _00.4`$. How realistic is this value? The optical depth of the entire column of matter between the observer and a source at $`z_s=1.0`$ (3.0) is $`\kappa _0=0.065`$ ($`0.34`$). It is very hard to imagine how APM galaxies, with faint galaxy magnitudes of around 20-21 in typical optical bands can be faithful tracers of matter fluctuations at redshifts 1-3. An alternative is to keep $`\kappa _0`$ at 0.028, but require $`b0.1`$. This is inconsistent with observations that estimate $`b`$ on large scales to be within $`30\%`$ of unity (Peacock & Dodds, 1994; Gastanaga & Juszkiewicz, 2001). We will return to interpretation of our results in Section 6.

## 5 GRB–ACO cluster correlations

If nearby mass distribution is weakly lensing GRBs, then all tracers of the nearby mass should be anti-correlated with GRB events. Here we carry out cross-correlations between GRBs and galaxy clusters from the Abell-Colwin-Olowin Catalog (Abell, 1958).

### 5.1 Results from literature

The present work is not the first to consider GRB-ACO correlations. Prior to 1997, before the first spectroscopic redshifts to GRBs were established, several workers have looked for correlations between GRBs and ACO clusters. Their motivation was that if GRBs are formed in galaxies, then the degree of their association with known populations of galaxy clusters will place constraints on their redshift distribution. All these studies tacitly assumed that weak lensing effect was negligibly small. Based on the current cosmological models, and our predictions in Section 4.3, their assumption was perfectly justified.

Some studies detected correlations between subsets of clusters and subsets of GRBs, others detected no correlations. For example, Kolatt & Piran (1996) claimed an association of 136 $`e1.6^{}`$ GRBs from the 3B Catalog and 3616 $`|b|30^{}`$ ACO clusters at separations $`\mathrm{\Theta }4^{}`$, at a significance level of 95%. Marani et al. (1997) reported correlation of 71 $`e1.685^{}`$ 3B GRBs with 185 nearby ACO clusters with $`R1`$ and $`D4`$ at $`2.93.5\sigma `$ level. They found even stronger correlations between 27 $`e<\mathrm{\hspace{0.17em}0.35}^{}`$ 3B GRBs and all 5250 ACO clusters, at $`3.54\sigma `$ level. Given these findings, they were surprised that 40 very well localized Inter Planetary Network (IPN) GRBs were not correlated with ACO clusters. Hurley et al. (1999) also used IPN positions but for a much larger sample, 157 GRBs from 4B Catalog. The average reduction in error area is $`50`$ for IPN compared to BATSE positions. Despite the accurate positions Hurley et al. (1999) did not find any correlations, for assumed ACO cluster radii of $`0.2^{}`$ or $`0.4^{}`$.

When GRB redshifts became available a few years ago, GRB-cluster cross-correlation work stopped. In view of our APM galaxy results in Section 4, it becomes interesting to revisit GRB–ACO correlations, especially since the Current BATSE Catalog contains about twice as many GRB events as were used in earlier studies, and error estimates for some GRB events have been recently revised (Paciesas et al., 1999).

### 5.2 Our results

Abell-Colwin-Olowin Catalog clusters are classified into distance and richness classes according to the standard Abell criteria (Abell, 1958). Since lensing should be sensitive to cluster distances we divide the clusters into three distance ranges: nearby, D=1-4, which have an average redshift of 0.06, medium distant, D=5, with an average $`z`$ of 0.11, and distant, D=6, with an average $`z=0.18`$. Distance class 7 is very incomplete; it has 8 clusters vs. about 2000 D=6 clusters, so we leave out D=7 clusters from the analysis. We further divide the clusters into poor ones, R=0-2, and rich ones, R$``$3. ACO catalog is incomplete close to the Galactic plane, so we mask out areas at $`|b|<30^{}`$. This leaves us with 3608 clusters and 1021 GRB events which we use to compute cross-correlations on a range of angular scales.

Fig. 8 shows GRB-ACO correlation functions between GRBs with $`e1^{}`$ 222This selection corresponds to the 74 GRB subset of the APM analysis, however, the number of GRBs here is actually 183. The difference arises mostly because in the APM analysis we did not use GRBs that were too close or too far from the plate center, whereas here there are no such constraints. and six subsets of ACO clusters. We used Landy & Szalay (1993) estimator, which has improved variance compared to the standard estimator. The two sets of solid lines in each panel are the 95% and 99% confidence limits derived from cross-correlations using 300 random realizations of the GRB catalog (compensated for the same declination dependent sky coverage as the real BATSE catalog). The dotted lines are the same correlations, but binned into narrower angular bins. (We do not plot confidence limits for these).

Out of sixty points plotted in six panels of Fig. 8 only one is significant at $`99\%`$: it states that GRBs with $`e1^{}`$ are anti-correlated with rich distant clusters on scales $`<20^{}`$. Of the six cluster subsets, the rich distant ones are most likely to act as weak lenses for background GRB, so this result supports our earlier GRB-APM findings. Whereas in the APM study the angular scales we could probe were limited by the plate size, here we are not restricted in that regard. We find that the full scale of anti-correlations corresponds to about 150$`h^1`$Mpc at the redshift of the clusters, which is the scale where cluster power spectrum is observed to peak: Retzlaff et al. (1998) find peak $`k=2\pi /l0.05h`$ Mpc-1 for ACO clusters, while Tadros et al. (1998) find peak $`k0.03h`$ Mpc-1 for APM clusters.

Results with $`0.5^{}e1^{}`$ (corresponding to the 46 GRB subset of the APM analysis) are similar to the ones in Fig. 8, however, the significance of anti-correlations with rich distant clusters drops to $`97\%`$. Correlations with GRBs unrestricted in positional error show no signal.

## 6 Discussion

Here we explore ways of accommodating GRB-APM anti-correlations, and earlier reports of QSO-galaxy correlations (Williams & Irwin, 1998; Norman & Williams, 2000; Bartsch et al., 1997). There are three possible avenues towards a resolution of the high amplitude of (anti-)correlations; these concern the sources (GRBs or QSOs), the lensing process, and the lenses (the mass distribution), respectively. We discuss these in turn. We assume that in the weak lensing regime the effects of these three sets of possibilities add up ‘linearly’, and so they can be considered separately.

### 6.1 Sources

Suppose our assumptions about GRB number distribution in fluence and error (Fig. 1) are incorrect, i.e. our assumed conversion between $`M`$ and $`q`$—exponent $`\beta `$ in eq. 4—is wrong. This could arise, for example, due to shot noise given small number of GRBs. The extreme value that $`\beta `$ can take is $`1`$, which is when there are no more sources beyond the limit of the survey, and area dilution of lensing reduces the sky density of sources by $`M`$. Even this extreme case does not come close to reproducing GRB-APM anti-correlations.

A similar solution to the problem was considered in Williams & Irwin (1998), who observed positive correlations between bright, optically selected QSOs and APM galaxies. The QSO number counts in the relevant redshift range gave $`q=M^{1.75}`$, whereas $`q=M^{20}`$ would be required to reproduce observations. That explanation was ruled out as unlikely.

### 6.2 Lensing process

Is it possible that propagation of light through an inhomogeneous Universe is not well modeled by the standard lensing equation? In terms of eq. 4 that would mean that the expression in round parentheses is incorrect. This was considered in Williams (2000), who explored the effect of the second order term in the lensing equation. Including that term had interesting consequences: magnification was increased, but only by 10%, and only in the mass distribution scenarios which were far from Gaussian random fields on large spatial scales. So this does not appear to be a promising avenue for the resolution of the problem. Furthermore, propagation of light through a numerically simulated Universe of popular cosmological models was done using the full geodesic equation of motion (Van Waerbeke et al., 2001; Tomita et al., 1999), and the results were found to be in good agreement with the predictions based on the standard lensing equation.

### 6.3 Lenses—constant biasing

The last set of possibilities is that the mass fluctuations are quite different from the observed projected galaxy density fluctuations. To reproduce the GRB-APM anti-correlations, $`\kappa _0/b`$ in eq. 4, would have to be increased by a factor of 10-20. Increasing $`\kappa _0`$ by that large a factor is ruled out (Section 4.3). On the other hand, a constant biasing of factor of 0.1-0.05 is ruled out by dynamical measurements on cluster and supercluster scales.

### 6.4 Lenses—density dependent biasing

Here we propose another variant of the third set of possibilities; we relax the requirement that $`b`$ is constant as a function of $`\delta \sigma /\sigma `$. This is a toy model only. We keep $`\kappa _0=0.028`$, which is appropriate for the average optical depth of the APM galaxies. Within this slab of APM galaxies, the projected galaxy density is assumed to be a monotonic, but not linear tracer of the projected mass density.

The shape of the distribution of APM galaxy counts-in-cells on $`\theta =1.5^{}`$ scales is fixed by observations. To maximize lensing effects we chose a skewed shape for $`p(\delta \sigma /\sigma )`$: let the distribution of projected mass densities averaged over circles of $`\theta =1.5^{}`$ have the shape of a half-Gaussian, with the sharp cut-off coinciding with $`\delta \sigma /\sigma =1`$, and the tail extending to positive values of $`\delta \sigma /\sigma `$. In terms of the optical depth the cutoff is at $`\kappa _0`$ (see insert in Fig. 6); the corresponding lines of sight are ‘empty beams’ and produce maximum possible demagnification of the sources. The width of the half-Gaussian is adjusted such that the average $`\delta \sigma /\sigma `$, and hence the average $`\kappa `$ are zero.

Since $`p(\delta \sigma /\sigma )`$ is highly asymmetric, while the observed $`p([n_{gal,D}/n_{gal,R}]1)`$ is symmetric, we must adjust the biasing function to map the former distribution onto the latter. Fig. 7 plots one possible version of such a biasing function, and the corresponding distribution of projected galaxy densities is shown as the dashed line in Fig. 6. It is symmetric, and nearly zero-centered, i.e. similar to the observed APM distribution (the solid histogram in this figure is the same as the solid histogram in Fig. 5). In short, the biasing function takes the very skewed projected mass distribution (insert in Fig. 6) and maps it onto a symmetric projected galaxy distribution (dashed histogram in Fig. 6).

Within this toy model, the average galaxy density for the 46 and 74 synthetic GRB subsets is about 0.97, and densities $`0.883`$ ($`0.915`$) for the 46 (74) subset occur in 9% (13%) of the cases. Thus the model can account for the observations. Furthermore, this model also reproduces the observed correlations reported in Williams & Irwin (1998): if the slope of QSO number counts is $`\alpha =1.1`$ (as quoted in that study), our model gives 1.018 for the average normalized galaxy foreground density in circles $`\theta =1.5^{}`$ around bright optically selected QSOs, a value consistent with Fig. 6 of that study ($`1.01`$).

Our model $`p(\kappa )`$ distribution and the resulting biasing function are very different from what is expected in currently accepted Universe models, with initial density perturbations specified by Harrison-Zel’dovich-type spectrum, and matter and cosmological constant (or the like) contributing comparably to the net zero curvature. In such models $`p(\kappa )`$ is Gaussian on large scales, and the biasing factor is close to unity for all values of density excess. For example, Fig. 13 of Jain et al. (2000) show that the symmetric shape for $`p(\kappa )`$ is already attained at smoothing scale of $`8^{}`$, and is symmetric for all larger scales.

Our model requires that there are lines of sight with projected radii $`20h^1`$ Mpc and extending from $`z0.1`$ to $`0.4`$, or $`600h^1`$ Mpc in terms of proper length, that are nearly devoid of mass. At the same time, these lines of sight are not nearly as empty in terms of galaxies. Within standard theoretical framework, where primordial fluctuation power spectrum and subsequent gravitational instability are solely responsible for structure on large scales, a dramatic redistribution of baryonic and dark matter, such as implied by our model, is not possible. Note however, that even though our model is astrophysically implausible, it does not violate the physical constraint that densities must remain positive everywhere: the sharp cut-off in $`p(\kappa )`$ corresponds to $`\kappa _0`$, the minimum possible lensing optical depth presented by mass between $`z=0.1`$ and $`0.4`$ to sources at $`z_s2.5`$, in a flat $`\mathrm{\Omega }=0.3`$ universe model.

Since baryonic and dark matter are unlikely to be distributed differently on large scales, we propose an alternative, speculative scenario. Suppose $`70\%`$ of the present day closure density is contributed by some dark energy that, because of its equation of state became dynamically important only recently, say $`z<\mathrm{\hspace{0.17em}1}`$. Until $`z1`$ structure formation proceeded according to the standard picture, with galaxies tracing the total mass, and with galaxy and mass distributions looking Gaussian on large scales. Suppose also that the dark energy can cluster on sub-horizon scales, or, more specifically, on scales as small as $`100200h^1`$ Mpc. Then at late epochs the dark energy will start contributing to the potential fluctuations that are already defined by dark and baryonic matter, on those scales. If these fluctuations develop quickly, then matter (light and dark), hampered by inertia, would not be able to respond quickly. Thus, power spectra obtained from galaxy surveys would not betray anything unusual. On the other hand, light rays from distant sources will traverse the fluctuations at $`c`$, thereby ‘capturing’ the total amplitude of the fluctuations. So, lensing is able to probe the full extent of gravitational potential wells, and if these are deep, then significant source-lens (anti-)correlations would result.

In this picture, $`p(\kappa )`$ distribution introduced earlier would refer not to the dark matter, but mainly to the projected clumping of dark energy. Assuming 3D clustering scale of $`100200h^1`$ Mpc, lines of sight $`600h^1`$ Mpc long that are devoid of dark energy are possible. Given the speculative nature of the proposed scenario, we will not develop it any further. Instead, we ask if the qualitative features of our model are compatible with the existing observations. By qualitative features we mean the relation between fractional galaxy excess and fractional dark energy excess (Fig. 7), with biasing now redefined in terms of these two quantities. Specifically, do observations rule out strong anti-biasing for regions with small $`|\delta \sigma /\sigma |`$, biasing for regions with $`\delta \sigma /\sigma >`$few, and moderate anti-biasing or $`b1`$ for all other regions. In 3D, this general behavior of the biasing function will be same as in projection.

Cosmic Microwave Background. According to the proposed scheme the dynamical evolution at very high redshifts is same as in the standard cosmological scenarios, so the primary CMB fluctuations would remain unaltered. The dark energy fluctuations grow at $`z<\mathrm{\hspace{0.17em}1}`$, and are therefore expected to contribute to CMB as late integrated Sachs-Wolfe effect, at $`\mathrm{}180^{}/\theta 100200`$. This location coincides with that of the first acoustic peak, and so the ISW signature will be either masked by the peak, or enhance the peak’s amplitude compared to that of the secondary peak. This is not ruled out by CMB observations.

Galaxy redshift surveys. Because of the late emergence of dark energy fluctuations, and galaxies’ sluggishness in responding to these, the observed galaxy distribution at redshifts $`<\mathrm{\hspace{0.17em}0.5}1`$ should be similar to what one would expect in standard cosmological models.

Cosmic velocity flows. Regions with $`|\delta \sigma /\sigma |<\mathrm{\hspace{0.17em}1}`$ would show significant anti-biasing according to our model. Such regions within a few hundred Megaparsecs around us are studied by cosmic velocity flows, and in the linear regime yield values for $`\mathrm{\Omega }_m^{0.6}/b`$ (Dekel, 1994). There are three basic methods of reconstructing 3D mass density fields from data: using density-density comparison, velocity-velocity comparison, or redshift-space distortions. These methods often produce discrepant results for $`\mathrm{\Omega }_m^{0.6}/b`$, with density-density giving consistently higher values than velocity-velocity comparison (Berlind et al., 2001). For example, IRAS 1.2Jy + Mark III data set yields $`\mathrm{\Omega }_m^{0.6}/b=0.89\pm 0.12`$ using density-density technique (Sigad et al., 1998), while $`\mathrm{\Omega }_m^{0.6}/b=0.50\pm 0.04`$ is obtained using velocity-velocity comparison (Willick & Strauss, 1998). Other data sets produce similar results. Berlind et al. (2001) argue that no reasonable biasing scheme would generate such different $`\mathrm{\Omega }_m^{0.6}/b`$ emerging from different methods, and ascribe the discrepancy to errors. This may well be the case, however, the discrepancy could also arise if biasing is a strong function of total density excess, as proposed by our toy model. In that case, density-density comparison would yield most reliable results, because the method makes no assumptions about the underlying mass density field, but derives it from observed velocities. Density-density comparisons currently yield $`\mathrm{\Omega }_m^{0.6}/b1`$ (Sigad et al., 1998); for $`\mathrm{\Omega }_m0.3`$ this would imply $`b0.5`$, not too different from the average biasing factor of our proposed model, which is 0.62. Thus our model is not necessarily in conflict with observed cosmic velocity flows.

Cosmic weak shear lensing. Another important observational test of our model is provided by the recent measurements of weak shear (Bacon et al., 2000; Hoekstra et al., 2002; Kaiser et al., 2000; Maoli et al., 2001; Van Waerbeke et al., 2000; Wittman et al., 2000), which are consistent with $`\mathrm{\Omega }_m=0.3`$ flat cosmological models, $`b1`$333To be more precise, Mellier et al. (2001) derive an empirical fit: $`\mathrm{\Omega }_m^{0.47}\sigma _80.59_{0.03}^{+0.03}`$., and Gaussian mass distribution on large scales. Even though the total mass and dark energy distribution in our model is far from Gaussian, our model and weak shear observations may not be discrepant; we list three possibilities: (i) The scales of weak shear and GRB-APM anti-correlations are different by about a factor of 3: the former apply to $`\mathrm{\Theta }<30^{}`$, while the latter are on scales $`<1.5^{}`$, thus the two techniques probe different physical scales. (ii) Observational factors may play a role. The redshifts of individual galaxies in the weak shear studies are not known, and so the lenses cannot be cleanly separated from the sources. Because of that, weak shear method may suffer from a ‘signal dilution’ effect: a foreground overdense patch will produce coherent shearing in background galaxies, but will also dilute the shear signal because of the increased number of foreground unsheared galaxies. (iii) In models with Gaussian distributed matter rms values of convergence and shear are the same. Suppose the distribution of dark energy trapped in potential wells of large overdense regions is top-hat-like. In that case the rms of convergence can be a few times larger than rms of shear, when averaged over large portions of the Universe. Thus the current weak shear observations may be underestimating the total amplitude of fluctuations.

We conclude that the qualitative features of our proposed model are not strongly ruled out by observations.

## 7 Summary and Discussion

As discussed in the Introduction, bright optically-selected and radio-selected QSOs are known to be correlated with foreground galaxies on large angular scales. The statistical significance of any one of these studies is generally $`23\sigma `$, and can be dismissed as a chance occurrence, but collectively these studies imply that the correlations are real, in which case the most likely explanation is weak magnification lensing. The purpose of this work was to see if cosmologically distributed sources other than QSOs are affected by weak lensing.

We started by computing the expected magnification bias of GRB events due to weak lensing, and found that most GRB subsets limited by statistical positional error should show negative magnification bias, i.e. should be anti-correlated with the foreground mass. Using APM galaxies to trace the mass we then looked for any correlations on the sky between APM galaxies and GRBs. Because the projected number density of GRBs is low, about 1 per APM plate, we did not use cross-correlation analysis of counting pairs, but instead counted galaxies in circles around individual GRBs. We found that GRBs with small positional errors are preferentially located in the directions where APM galaxies show $`10\%`$ deficits on degree angular scales. In particular, a subset of 46 GRBs with $`0.5^{}e\mathrm{\hspace{0.17em}1}^{}`$ and a subset of 74 GRBs with $`0^{}e\mathrm{\hspace{0.17em}1}^{}`$ have 12% and 8% less galaxies in circles of $`1.5^{}`$ radii around them (or $`20h^1`$ Mpc at the redshift of the typical lenses), than expected on average. These are significant at 99.7% and 97.6% confidence levels. Whereas these significance levels are less than overwhelming, what makes the finding especially interesting is that the GRBs are anti-correlated with the APM galaxies, as expected from the magnification bias (and opposite to what is expected if physical associations are significant). This is the first reported case of weak lensing induced anti-correlations on degree angular scales.

To test the observed anti-correlation of GRBs with foreground mass, we carried out a cross-correlation analysis between GRB subsets classed by positional error and Abell-Corwin-Olowin galaxy clusters of three distance ranges, D=1-4, D=5, and D=6, and two richness subsets, R=0-2, and R$`3`$. Only one combination of GRBs and clusters showed significant results: $`0^{}e\mathrm{\hspace{0.17em}1}^{}`$ GRBs and rich distant ACO clusters are anti-correlated on $`<20^{}`$ scales, consistent with predictions in Section 3, and observed anti-correlations with APM galaxies.

As was the case with earlier reports of QSO-galaxy correlations, present anti-correlations with GRBs cannot be accounted for by weak lensing using a standard model of mass distribution and moderate biasing on large scales; the discrepancy is rather severe, one needs to increase the lensing optical depth by 10, or reduce the biasing factor to $`b0.1`$, which is not a viable option. Since both the GRB-APM and GRB-ACO anti-correlation signals are of 2-3 $`\sigma `$ significance, they could well be statistical flukes. However, if they are not, and if the signal is due to weak lensing, then there are three types of resolutions: those that have to do with sources, lensing process, or lenses. We argue that the first two are inadequate, leaving us with the third.

We propose a speculative scenario which assumes that $`70\%`$ of closure density is contributed by dark energy that can clump on significantly sub-horizon scales, $`\mathrm{\hspace{0.17em}100}200h^1`$ Mpc, and that this clumping developed recently, $`z<\mathrm{\hspace{0.17em}1}`$. Then at present epochs there should be significant fluctuations in the dark energy component. However, the amplitude of baryonic and dark matter fluctuations would remain relatively unchanged, because it would take time for matter to be accelerated and displaced by significant amounts on large spatial scales. In this scenario, the galaxy distribution would trace the underlying potential wells but would severely underestimate the amplitude of fluctuations on $`>100h^1`$ Mpc scales. Dynamical measurements would do better, since velocities ‘respond’ quicker to acceleration than displacements. Lensing, which relies on light traversing the fluctuations at $`c`$, is the best tool to probe the full extent of these fluctuations. In this scheme, weak lensing induced GRB-APM anti-correlations would be strong because the total amplitude of gravitational potential fluctuations are substantially larger than galaxies would have us believe.

In principle, weak shear lensing should be able to detect the true extent of these fluctuations as well. However, weak shear observations support the standard cosmological model with cold dark matter and unclustered dark energy component, and seem to be in conflict with our proposed scenario. Current observations of microwave background anisotropy, galaxy redshift surveys and cosmic velocity flows are not strongly incompatible with our hypothetical scenario.