text:
# Networks with desired correlations

## Abstract

We discuss a simple method of constructing correlated random networks, which was recently proposed by M. Boguñá and R. Pastor-Satorras (cond-mat/0306072). The result of this construction procedure is a sparse network whose degree–degree distribution asymptotically approaches a given function at large degrees. We argue that this convergence is possible if the desired function is sufficiently slowly decreasing.

Recently, Marián Boguñá and Romualdo Pastor-Satorras proposed a natural way of constructing correlated networks with desired correlations of the degrees of the nearest neighbors . This practical algorithm is based on the idea that previously was used in Refs.  for constructing uncorrelated graphs and in Ref.  for building correlated networks (see also Ref. ). In all these constructions, (i) some weights (fitnesses, desired degrees, etc.) are ascribed to vertices, and (ii) each pair of vertices is connected with probability which depend on these weights.

In the first version of this work, I reproduced this algorithm without knowing that it had already been proposed. So, in the present, corrected and shortened and more methodical version I mostly discuss the range of validity of the algorithm and features of networks, generated by this method.

The correlations between degrees of the nearest neighbors in a graph are naturally described by the joint distribution of the degrees of end vertices of an edge of the graph, $`P(k,k^{})`$, $`P(k,k^{})=P(k^{},k)`$, $`_{k,k^{}}P(k,k^{})=1`$. The joint degree–degree distribution determines the degree distribution $`P(k)`$ of a network:

$$\underset{k^{}}{}P(k,k^{})=\frac{kP(k)}{\overline{k}}.$$
(1)
Consequently,

$$\overline{k}=\left[\underset{k,k^{}}{}\frac{P(k,k^{})}{k}\right]^1$$
(2)
and $`k^n=\overline{k}_{k,k^{}}k^{n1}P(k,k^{})`$.

The algorithm of Boguñá and Pastor-Satorras generates sparse random networks with desired degree–degree correlations. Suppose one wishes to obtain an ensemble of graphs with a desired joint distribution of the degrees of the nearest neighbors, $`P(q,q^{})`$. One should assume that $`P(q,q^{})`$ decreases with $`q`$ and $`q^{}`$ sufficiently slowly. Let the number $`N`$ of vertices in each graph of the ensemble be large and fixed.

The procedure  is as follows:

(i) Create $`N`$ vertices with a sequence of weights $`\{q_i\}`$, $`i=1,\mathrm{},N`$ independently sampled from the distribution

$$P(q)=\overline{q}\underset{q^{}}{}P(q,q^{})/q,$$
(3)
where $`\overline{q}=[_{q,q^{}}P(q,q^{})/q]^1`$.

(ii) Put a link between $`i`$ and $`j`$ vertices with probability

$$p_{ij}=p(q_i,q_j)=\frac{\overline{q}}{N}\frac{P(q_i,q_j)}{P(q_i)P(q_j)}.$$
(4)
In this model, (i) the average degree of an $`i`$-th vertex coincides with $`\overline{q}_i`$, and (ii) the distribution of degrees of an individual vertex is a relatively narrow function at large degrees. This can be proved by calculating the degree distribution of a vertex with weight $`q`$ (see Ref. ). Here, alternatively, we simply find the first two moments, $`\overline{k}_i`$ and $`k_i^2`$ for the degree of an $`i`$-th vertex.

The statistical weights of graphs $`gG`$ in the resulting ensemble, written in terms of $`p_{ij}`$ (adjacency matrix elements are $`a_{ij}=0,1`$), are of a rather standard form

$$\mathrm{\Pi }(g)=\mathrm{\Pi }(\{p_{ij}\},\{a_{ij}\})\underset{i,j}{}\left(\frac{p_{ij}}{1p_{ij}}\right)^{a_{ij}}=\underset{i,j}{}s_{ij}^{a_{ij}}$$
(5)
(compare with the statistical weights of classical random graphs and Ref. ). Here, $`s_{ij}=p_{ij}/(1p_{ij})`$. The partition function of the ensemble is

$`Z(G)=Z(\{p_{ij}\})`$ $`={\displaystyle \underset{gG}{}}\mathrm{\Pi }(g)={\displaystyle \underset{a_{ij}}{}}{\displaystyle \underset{i,j}{}}\mathrm{\Pi }(\{p_{ij}\},\{a_{ij}\})`$ (7)
$`{\displaystyle \underset{i,j}{}}(1p_{ij})^1={\displaystyle \underset{i,j}{}}(1+s_{ij}).`$
So,

$`a_{ij}={\displaystyle \frac{\mathrm{ln}Z(\{s_{ij}\})}{\mathrm{ln}s_{ij}}}=p_{ij},`$ (8)
$`a_{ij}a_{i^{}j^{}}a_{ij}a_{i^{}j^{}}={\displaystyle \frac{}{\mathrm{ln}s_{i^{}j^{}}}}{\displaystyle \frac{}{\mathrm{ln}s_{ij}}}\mathrm{ln}Z(\{s_{ij}\})=`$ (9)
$`(p_{ij}p_{ij}^2)\delta _{ii^{}}\delta _{jj^{}},`$ (10)
etc. \[The averages are over the statistical ensemble: $`X(g)=Z(G)^1_{gG}\mathrm{\Pi }(g)X(g)`$.\] Consequently, all the moments are equal: $`a_{ij}^n=p_{ij}`$, $`n1`$. This, in fact, is clear, since $`a_{ij}`$ takes only two values $`0`$ and $`1`$ with the probabilities $`1p_{ij}`$ and $`p_{ij}`$, respectively. Then,

$$\overline{k}_i=\underset{j}{}a_{ij}=\underset{j}{}p_{ij}=N\underset{q}{}P(q)p(q_i,q)=\overline{q}_i$$
(11)
\[Eq. (4) was used\], so the mean degree in the net $`\overline{k}=\overline{q}`$, and

$$k_i^2=\left(\underset{j}{}a_{ij}\right)^2=\underset{j}{}p_{ij}+\left(\underset{j}{}p_{ij}\right)^2\underset{j}{}p_{ij}^2.$$
(12)
Note that Eqs. (11) and (12) may be obtained by using more naive arguments. Suppose for brevity that a net is of three vertices, $`0`$, $`1`$, and $`2`$. Then $`\overline{k}_0=1[p_{01}(1p_{02})+(1p_{01})p_{02}]+2p_{01}p_{02}=p_{01}+p_{02}`$ \[compare with Eq. (11)\] and $`k_0^2=1^2[p_{01}(1p_{02})+(1p_{01})p_{02}]+2^2p_{01}p_{02}=p_{01}+p_{02}+2p_{01}p_{02}`$ \[compare with Eq. (12)\].

Using the expression (4), one can see that in a sparse network, the last term on the right-hand side of Eq. (12) is of the order of $`N^1`$ and so is negligible. Consequently, in the thermodynamic limit ($`N\mathrm{}`$),

$$k_i^2\overline{k}_i^2=\overline{k}_i=\overline{q}_i.$$
(13)
That is, the distribution of degrees of an individual vertex $`i`$ is indeed relatively narrow, if $`k_i`$ is large, and peaked at $`\overline{q_i}`$. One should take into account the fact that the number of edges in the net, which connect vertices with weights $`q`$ and $`q^{}`$, is

$`L(q,q^{})`$ $`=`$ $`{\displaystyle \frac{1}{2}}N^2P(q)P(q^{})a(q,q^{})=`$ (15)
$`{\displaystyle \frac{1}{2}}N\overline{q}P(q,q^{})=LP(q,q^{}),`$
where $`L`$ is the total number of edges \[note the symmetry factor $`1/2`$ in the second term of Eq. (15)\]. Then one can finally conclude that the resulting degree–degree distribution of the net is expressed in terms of a given function:

$$P_{\text{result}}(k,k^{})=P(k+\delta (k),k^{}+\delta ^{}(k^{})),$$
(16)
where the deviations, $`|\delta (k)|\sqrt{k},|\delta ^{}(k^{})|\sqrt{k^{}}`$, are relatively small at large degrees. The relation (16) shows that at large degrees the algorithm of Boguñá and Pastor-Satorras provides the degree–degree distribution $`P_{\text{result}}(k,k^{})P(k,k^{})`$, only if $`P(k,k^{})`$ decreases sufficiently slowly. The following arguments explain what this does mean.

We ask two (related) questions. (i) How should a function $`f(x)`$ behave at large $`x`$ to satisfy the condition: $`f(x+c\sqrt{x})/f(x)1`$ as $`x\mathrm{}`$? (ii) How should a function $`f(x)`$ behave at large $`x`$ to guarantee that $`[f(x+c\sqrt{x})+f(xc\sqrt{x})]/[2f(x)]1`$ as $`x\mathrm{}`$? Here, $`c=\text{const}`$.

At first sight, the second, more symmetric condition may demand less strong restriction on the behavior of $`f(x)`$. However, both the questions have the same answer. Indeed, the second and the third, and the higher terms in the series $`f(x+c\sqrt{x})/f(x)=1+c\sqrt{x}[df(x)/dx]/f(x)+c^2x[d^2f(x)/dx^2]/[2f(x)]+\mathrm{}`$ approach zero at large $`x`$ if the same condition is satisfied: $`f(x)`$ must decrease with $`x`$ slower than $`e^\sqrt{x}`$. So, if the given function $`P(k,k^{})`$ decreases slower than, say, $`e^{\sqrt{k}\sqrt{k^{}}}`$, then $`P_{\text{result}}(k,k^{})`$ asymptotically approaches $`P(k,k^{})`$ at large $`k`$ and $`k^{}`$ .

Moreover, this condition also guarantees that the resulting degree distribution $`P_{\text{result}}(k)`$ asymptotically approaches $`P(k)`$ at large degrees. On the other hand, in the region of small degrees, $`P_{\text{result}}(k,k^{})`$ deviates from desired $`P(k,k^{})`$, and $`P_{\text{result}}(k)`$ deviates from $`P(k)`$.

Similar arguments are also valid for other network constructions of this type : desired distributions must decrease sufficiently slowly. On the one hand, this implies a serious restriction on the range of degree (or degree–degree) distributions which can be reproduced in such a way. On the other hand, it is the slowly decreasing distributions that are most interesting.

Electronic address: sdorogov@fc.up.pt


format: text

text:
# The shallow slope of the 𝑧∼6 quasar luminosity function: limits from the lack of multiple image gravitational lenses

## 1 Introduction

The quasar luminosity function (LF) is the most basic property of the quasar population. At low redshifts several decades of study have yielded a well defined optical quasar LF with powerlaw slopes at both the faint and bright ends that do not evolve out to $`z3`$ (e.g. Boyle, Shanks & Peterson 1988; Hartwick & Schade 1990; Pei 1995a; Boyle et al. 2000). At higher redshifts only very bright quasars are currently observable at the magnitude limit of large surveys, and recent evidence suggests that the slope of their LF is significantly shallower than observed at $`z3`$ (Schmidt, Schneider, Gunn 1995; Fan et al. 2001a). There are now six quasars known with redshifts $`z5.8`$ (Fan et al. 2001b,2003). These very high redshift quasars provide important constraints for studies of structure formation (Turner 1990; Haiman & Loeb 2001) and reionisation (e.g. Madau, Haardt & Rees 1999; Wyithe & Loeb 2003). Determination of their LF (Fan et al. 2003) is critical if we are to address these issues.

As a population of sources for gravitational lensing the $`z6`$ quasars are unique, being the only sample where the gravitational lensing probability may be of order unity (Wyithe & Loeb 2002a,b; Comerford, Haiman & Schaye 2002). The high expected lensing rate arises through a large magnification bias which increases the fraction of gravitational lenses at a given flux level by drawing sources from the fainter, more numerous population into a flux limited sample. As a result the fraction of quasars in a sample that are multiply imaged by gravitational lenses is sensitive to the slope of the LF. Conversely, the observed multiple image lens fraction may be used to limit the unknown slope of a LF. This exercise was undertaken by Fan et al. (2003). They presented likelihood functions for $`\beta `$ given the absence of multiply imaged quasars in the $`z6`$ sample, and found the lack of lenses to be surprising at the $`90\%`$ level if $`\beta 3.5`$. This lensing constraint is consistent with their findings for $`\beta `$ through direct analysis of the luminosity distribution in the sample which yielded $`\beta 3.35`$ ($`90\%`$). However there is a complication. The lens fraction is linearly related (nearly) to the efficiency of the lens population, which is proportional to the expectation value of the velocity dispersion to the fourth power and to the space density of galaxies. Moreover, the lensing rate requires extrapolation of the local galaxy population to higher redshifts, and is also sensitive to cosmology (e.g. Kochanek 1996). Thus the unbiased lensing cross-section is quite uncertain.

While high resolution imaging data (Fan et al. 2001b,2003) shows none of the six $`z6`$ quasars to be multiply imaged, galaxies have been detected near the line of sight to two quasars (Shioya et al. 2002; White et al. 2003). This is puzzling because magnification bias should result in highly magnified multiply imaged sources being over represented among a population of quasars whose images are located near foreground galaxies. The effect becomes larger as the quasar LF becomes steeper. In this paper we present a Bayesian analysis that employs information on a posteriori alignments of the two lensed quasars. This statistic is much less sensitive to the uncertainties in the lens cross-section than the fraction of multiply imaged quasars, and we show that it produces a tighter limit on $`\beta `$.

The paper is set out as follows. In § 2 we compute the gravitational lens cross-section in light of the recently measured velocity function of galaxies, the probability of multiple imaging for different LFs and the limits on the quasar LF that result from the observed lack of multiply imaged quasars. We then discuss a Bayesian approach to computing the lens fraction in § 3. In § 4 we describe the two high redshift quasars thought to be magnified by gravitational lensing. The probability of getting a multiply imaged quasar within a sub-sample of quasars observed to be near a lens galaxy is discussed in § 5. These probabilities are used to compute likelihood functions for the fraction of lensed singly imaged quasars, and to derive limits on $`\beta `$. In § 6 we use the redshifts of the lens galaxies to constrain simple parametric models for the evolution of the lens galaxy population. We then combine these results with the limits on $`\beta `$ to estimate the expected multiple image fraction for $`z6`$ quasars in § 7. Finally in § 8 we discuss the implications of possible multiple imaging in SDSS J1148-5251 for $`\beta `$ before presenting our conclusions in § 9. Where required we assume the most recent cosmological parameters obtained through fits to WMAP data (Spergel et al. 2003). These include density parameters of $`\mathrm{\Omega }_m=0.27`$ in matter, $`\mathrm{\Omega }_b=0.044`$ in baryons, $`\mathrm{\Omega }_\mathrm{\Lambda }=0.73`$ in a cosmological constant, and a Hubble constant of $`H_0=71\mathrm{km}\mathrm{s}^1\mathrm{Mpc}^1`$.

## 2 The gravitational lens cross-section and multiple imaging rate

The probability that a Spherical Singular Isothermal (SIS) galaxy will lens a back-ground source is proportional to the 4th power of its velocity dispersion $`\sigma `$. Thus what is required to compute the cross-section for gravitational lensing is the *velocity function*. Until recently the velocity function for early type galaxies (the dominant lens population, Kochanek 1996) had to be computed through combination of a galaxy LF with the Faber-Jackson (1976) relation. However this procedure ignores the intrinsic scatter in the Faber-Jackson (1976) relation and is an unreliable method (e.g. Kochanek 1993). A more reliable representation is now possible using the measured velocity dispersion function of early type galaxies. Sheth et al. (2003) presented the measured velocity dispersion function for early type galaxies. They suggested an analytic fit of the form

$$\varphi (\sigma )d\sigma =\varphi _{}\left(\frac{\sigma }{\sigma _{}}\right)^\alpha \frac{\mathrm{exp}[(\sigma /\sigma _{})^\beta ]}{\mathrm{\Gamma }(\alpha /\beta )}\beta \frac{d\sigma }{\sigma },$$
(1)
where $`\varphi _{}`$ is the number density of galaxies and $`\sigma _{}`$ is a characteristic velocity dispersion. Sheth et al. (2003) found that the parameters $`\sigma _{}`$, $`\alpha `$ and $`\beta `$ are strongly correlated with one another. From Sheth we take $`\varphi _{}=(2.0\pm 0.1)\times 10^3`$ ($`h_{70}^1`$Mpc)<sup>-3</sup>, $`\alpha =6.5\pm 1.0`$, $`\beta =(14.75/\alpha )^{0.8}`$ and $`\sigma _{}=161\mathrm{\Gamma }(\alpha /\beta )/\mathrm{\Gamma }[(\alpha +1)/\beta ]\mathrm{km}\mathrm{s}^1`$.

Given the Einstein Radius (ER) for an SIS

$$\xi _0=4\pi \left(\frac{\sigma }{c}\right)^2\frac{D_\mathrm{d}D_{\mathrm{ds}}}{D_\mathrm{s}}$$
(2)
and a constant co-moving density of galaxies, the multiply imaged gravitational lens cross-section for sources at $`z_\mathrm{s}`$ is

$$\tau (z_\mathrm{s})=_0^{z_\mathrm{s}}𝑑z_0^{\mathrm{}}𝑑\sigma (1+z)^3\varphi (\sigma )\frac{cdt}{dz}\pi \xi _0^2.$$
(3)
We found the mean and twice the variance for a set of $`\tau `$ computed assuming the parameters $`\varphi _{}`$ and $`\alpha `$ to be distributed as Gaussian within their quoted uncertainties. This procedure yields $`\tau (z_\mathrm{s}=6)=(2.5\pm 0.25)\times 10^3`$. The statistical uncertainty of 10% is significantly smaller than is obtained through use of the Faber-Jackson (1976) relation and a LF. The value of $`\tau (z_\mathrm{s}=6)`$ obtained from the velocity function is a factor of $`3`$ smaller than obtained in estimates of the lens fraction in the very high redshift quasar samples (Wyithe & Loeb 2002a,b; Comerford, Haiman & Schaye 2002). While the implied lens fraction for these samples is still expected to be an order of magnitude higher than in lower redshift samples, it may result in less than one lens being expected in the current sample, hence limiting the use of the lens fraction for constraining the slope of the LF (Wyithe & Loeb 2002b; Comerford, Haiman & Schaye 2002; Fan et al. 2003).

As shown above the use of a measured velocity function reduces the statistical uncertainty in $`\tau `$ to 10%. However the remaining systematic error due to the uncertain evolution in the lens population (e.g. Keeton 2002) makes limits based on the lens fraction uncertain. For example, if the co-moving density of galaxies drops in proportion to $`(1+z)^\gamma `$ (density evolution) then for $`\gamma =1`$ we find $`\tau (z_\mathrm{s}=6)=(0.9\pm 0.09)\times 10^3`$, while for $`\gamma =2`$ the value drops to $`\tau (z_\mathrm{s}=6)=(0.4\pm 0.04)\times 10^3`$. Unless otherwise specified we assume a constant co-moving density of galaxies ($`\gamma =0`$). We consider only spherical lenses in this paper. Previous studies (Kochanek & Blandford 1987; Blandford & Kochanek 1987) have found that the introduction of ellipticities $`0.2`$ into nearly singular profiles has little effect on the lensing cross-section and image magnification. However the strong magnification bias will favour a high fraction of 4-image lenses (Rusin & Tegmark 2001). Finally we note that the $`i`$-band dropout quasars are selected independent of morphology, and so do not select against lenses, though the possibility that the lens galaxy itself will prevent detection should be accounted for (Wyithe & Loeb 2002b).

With these points in mind we compute the a-priori probability of multiple imaging in a sample of six quasars with $`z6`$ for different LFs. A similar calculation has already been performed by Fan et al. (2003), though we add the improvement of a more accurate $`\tau `$ computed from the velocity function, extend the calculation to account for constraints from the luminosity distribution (Fan et al. 2003), and provide a posteriori limits on $`\beta `$. We also modify $`\tau `$ to include only those lens galaxies which would not have contaminated the $`i`$band dropout selection of the quasars as discussed in Wyithe & Loeb (2002b). The probability that a quasar with luminosity $`L`$ and magnification bias $`B`$ will be multiply imaged by a foreground galaxy is

$$p_{\mathrm{lens}}(L)\frac{B(L)\tau }{B(L)\tau +(1\tau )},$$
(4)
where $`p_{\mathrm{lens}}`$ is slightly overestimated because we have not included any magnification bias for single image quasars. The sum of magnifications of multiple images ($`\mu `$) formed by an SIS has a probability distribution of the form

$$\frac{dP_\mathrm{m}}{d\mu }=\frac{8}{\mu ^3}\text{for}\mu 2,$$
(5)
resulting in a magnification bias $`B(L)`$ for a SIS and a LF $`\mathrm{\Phi }(L)`$ of

$$B(L)=_2^{\mathrm{}}\frac{d\mu }{\mu }\frac{8}{\mu ^3}\frac{\mathrm{\Phi }(L/\mu )}{\mathrm{\Phi }(L)}.$$
(6)
The probability given a LF $`\mathrm{\Phi }`$, that in a sample of six quasars at $`z6`$, none will be lensed is

$$p_{\mathrm{nolens}}=\underset{i=1,N_\mathrm{q}}{}\left(1p_{\mathrm{lens}}(L_i)\right),$$
(7)
where $`N_\mathrm{q}=6`$ and the $`L_i`$ are the luminosities of the $`N_\mathrm{q}`$ quasars.

### 2.1 lens fractions for double powerlaw luminosity functions

A successful fit to the low redshift ($`z2`$) quasar LF is the double powerlaw (e.g. Boyle et al. 2000)

$$\mathrm{\Phi }(L)=\frac{\mathrm{\Phi }_0}{L^\alpha +L^\beta },$$
(8)
where $`\alpha `$ and $`\beta `$ (note we have defined $`\alpha `$ and $`\beta `$ to be negative) are the slopes of the faint and bright ends of the quasar LF respectively, and we have expressed the luminosity $`L`$ in units of the characteristic break luminosity. Use of the double power-law form for $`\mathrm{\Phi }(L)`$ implies that we must specify two additional parameters before deriving limits on $`\beta `$. Firstly, since high magnifications will draw quasars that are fainter than the break into the sample, the magnification bias will depend on $`\alpha `$. Moreover use of the double powerlaw requires specification of the quasar luminosity with respect to the break. The $`z6`$ quasars range in luminosity from $`M_{1450}=27.1527.90`$ (Fan et al. 2001b,2003). In this paper we specify the luminosity of the quasars relative to the LF in terms of the difference between the magnitude of the faintest quasar and that of the break ($`\mathrm{\Delta }M=2.5\mathrm{log}_{10}L`$). For example if the break were at $`M_{1450}=26.15`$ then $`\mathrm{\Delta }M=1`$.

The resulting probability $`p_{\mathrm{nolens}}`$ that in the sample of six $`z6`$ quasars, none will be multiply imaged (given $`\alpha =1.75`$, $`\mathrm{\Delta }M=4`$) is plotted as a function of $`\beta `$ in the left hand panels of figure 1 (thick dashed grey line). The lack of multiply imaged quasars is only surprising at the 50% level for values of $`\beta <3.3`$. The probability $`p_{\mathrm{nolens}}`$ is also tabulated in Table 1 for various values of $`\alpha `$, $`\beta `$ and $`\mathrm{\Delta }M`$. The lack of multiply imaged quasars is only surprising at the 10% level if $`\beta 4`$ and $`\mathrm{\Delta }M3`$ or if $`\beta 3.75`$ and $`\mathrm{\Delta }M4`$. We note that these probabilities (and those in figure 9 of Fan et al. 2003) represent a likelihood function for $`\beta `$ rather than direct to limits on $`\beta `$. We now turn to computation of these limits.

The posterior probability for $`\beta `$ is

$$\frac{dP}{d\beta }|_{\mathrm{N}_{\mathrm{mult}}=0}=Np_{\mathrm{nolens}}\frac{dP_{\mathrm{prior}}}{d\beta },$$
(9)
where $`\frac{dP_{\mathrm{prior}}}{d\beta }`$ is the prior probability for $`\beta `$, and N is a normalising constant. We assume that the prior probability for the slope is flat between two bounds $`\beta _{\mathrm{min}}`$ and $`\beta _{\mathrm{max}}`$, hence

$$\frac{dP_{\mathrm{prior}}}{d\beta }=(\beta _{max}\beta _{\mathrm{min}})^1.$$
(10)
The absence of multiply imaged quasars in the sample is surprising if $`\beta `$ is small because the multiple image magnification bias is large in that case, but is less surprising as $`\beta `$ is increased. Hence the likelihood function $`p_{\mathrm{nolens}}`$ is increasing with $`\beta `$ so that the fraction of multiply imaged quasars carries no information on its upper limit. As a result, the confidence with which a small value of $`\beta `$ can be excluded given the lensing observations alone is sensitive to $`\beta _{\mathrm{max}}`$. This dependence implies that a small probability $`p`$ given some $`\beta `$ that a sample of quasars would contain no lenses does not translate to a limit on $`\beta `$ at the $`100(1p)\%`$ level. A second, independent constraint that limits the upper bound on $`\beta `$ and hence lowers the dependence on $`\frac{dP_{\mathrm{prior}}}{d\beta }`$ is required.

### 2.2 the addition of constraints from the luminosity distribution

Fan et al. (2003) have derived the constraint $`\beta \pm \mathrm{\Delta }\beta `$ where $`\beta =2.3`$ and $`\mathrm{\Delta }\beta =0.8`$ from their analysis of the luminosity distribution of the $`z6`$ quasars. We take this bound to be the 1-$`\sigma `$ level of a Gaussian likelihood function for $`\beta `$. We have plotted this likelihood function for $`\beta `$ (normalised to a maximum of 1), as well as the differential and cumulative probability distributions for $`\beta `$ in the in the left-hand, central and right-hand panels of figure 1 (solid grey lines).

We now have two constraints on $`\beta `$, one from the fraction of multiply imaged quasars and one from the distribution of luminosities. Assuming these constraints to be independent (this assumption is discussed below), we find a joint likelihood function, and hence posterior probability distribution for $`\beta `$

$$\frac{dP}{d\beta }|_{\mathrm{N}_{\mathrm{mult}}=0}=Np_{\mathrm{nolens}}\mathrm{exp}\left(\frac{1}{2}\left(\frac{\beta \overline{\beta }}{\mathrm{\Delta }\beta }\right)^2\right)\frac{dP_{\mathrm{prior}}}{d\beta }.$$
(11)
The likelihood function is now normalisable because the distribution of luminosities constrains large values of $`\beta `$ (Fan et al. 2003), and is quite insensitive to the prior $`\frac{dP}{d\beta }|_{\mathrm{N}_{\mathrm{mult}}=0}`$ as a result. We plot the posterior differential and cumulative probability distributions in the central and right hand panels of figure 1. These are shown by the dashed light lines, and should be compared to the posterior probability distributions based on the distribution of luminosities alone (solid light lines). The lensing constraint disfavours smaller values of $`\beta `$, resulting in a narrower probability distribution for $`\beta `$. The most likely value for the slope is $`\beta 2.2`$. The addition of the constraint from the fraction of multiply imaged lensed quasars improves on the limits obtained by Fan et al. (2003) from the distribution of luminosities alone. For this choice of $`\mathrm{\Delta }M`$ and $`\alpha `$, the lack of multiply imaged quasars rules out $`\beta 3.0`$ at the 90% level. We also construct posterior cumulative probabilities for various values of $`\alpha `$ and $`\mathrm{\Delta }M`$ and plot the results in figure 2 (thick dashed light lines). The addition of lensing constraints significantly improve the LF limits provided that the quasars are not too close to the LF break. For a constant co-moving density of lens galaxies (upper two rows) we find $`\beta 3.23.0`$ (90%) for $`\mathrm{\Delta }M2`$.

It should be noted that the two constraints are not quite independent. In general gravitational lensing tends to flatten the slope of the quasar LF by drawing populous faint quasars into a bright quasar sample (e.g. Pei 1995b). However there are two reasons to think that this is not a problem within the very high redshift quasar sample. Firstly the objects that are lensed are not multiply imaged and so have magnifications smaller than $`2`$. Secondly, the average change of slope is $`\mathrm{\Delta }\beta 0.2`$ even in the most optimistic lensing scenario (where $`\beta 3.5`$) for the $`z6`$ quasars (Wyithe & Loeb 2002b).

## 3 Bayesian approach to computing the multiple image fraction

We may also use a Bayesian approach to compute the multiple image fraction. This will provide us with a natural framework within which to add additional a posteriori information on alignments of quasars with foreground galaxies in § 5. Consider sources with unlensed impact parameters (in units of the ER) $`y=x1`$ with associated magnifications $`\mu `$. We write the likelihood per logarithm of $`x`$ of observing a singly imaged lensed quasar (including magnification bias)

$$L_{\mathrm{single}}=x(x1)\frac{1}{\mu _{\mathrm{single}}}\frac{\mathrm{\Phi }(L/\mu _{\mathrm{single}})}{\mathrm{\Phi }(L)},$$
(12)
where the factor $`(x1)`$ accounts for the additional solid angle available at large $`y`$, and $`\mu _{\mathrm{single}}=x/(x1)`$. This likelihood may be compared to the corresponding average likelihood of observing a multiply imaged quasar

$`L_{\mathrm{mult}}`$ $`=`$ $`{\displaystyle _0^1}𝑑yy{\displaystyle \frac{1}{\mu }}{\displaystyle \frac{\mathrm{\Phi }(L/\mu )}{\mathrm{\Phi }(L)}}`$ (13)
$`=`$ $`{\displaystyle \frac{1}{2}}{\displaystyle _2^{\mathrm{}}}{\displaystyle \frac{d\mu }{\mu }}{\displaystyle \frac{8}{\mu ^3}}{\displaystyle \frac{\mathrm{\Phi }(L/\mu )}{\mathrm{\Phi }(L)}}={\displaystyle \frac{B(L)}{2}}.`$ (14)
The likelihood that a quasar will be singly imaged at $`x`$ rather than multiply imaged is therefore

$$p_{\mathrm{single}}(\beta |x)=\frac{L_{\mathrm{single}}}{L_{\mathrm{single}}+L_{\mathrm{mult}}}.$$
(15)
We may also calculate the posterior probability that a quasar will not be lensed

$$p_{\mathrm{nolens},\mathrm{bayes}}=_2^{\mathrm{}}𝑑xp_{\mathrm{single}}(\beta |x)\frac{dP_{\mathrm{prior}}}{dx},$$
(16)
where $`\frac{dP_{\mathrm{prior}}}{dx}`$ is the prior probability for the $`x`$. This quantity is one minus the lens fraction and may be approximated using the usual formula for the lens fraction $`\tau B`$ or more accurate forms such as equation (4).

For the sample of $`z6`$ quasars, the posterior probability distribution for $`\beta `$ is therefore

$`{\displaystyle \frac{dP}{d\beta }}`$ $`=`$ $`N\mathrm{exp}\left({\displaystyle \frac{1}{2}}\left({\displaystyle \frac{\beta \overline{\beta }}{\mathrm{\Delta }\beta }}\right)^2\right){\displaystyle \frac{dP_{\mathrm{prior}}}{d\beta }}`$ (17)
$`\times `$ $`\left[{\displaystyle \underset{i=1,6}{}}\left({\displaystyle _2^{\mathrm{}}}𝑑xp_{\mathrm{single}}(\beta |x,L_i){\displaystyle \frac{dP_{\mathrm{prior}}}{dx}}\right)\right],`$
where the prior probability distribution for $`x`$ can be computed from the derivative of the Poisson probability that a source lies within a circle of radius $`x1`$ around a randomly positioned galaxy

$$\frac{dP_{\mathrm{prior}}}{dx}=2\tau (x1)e^{\tau (x1)^2}.$$
(18)
Equation (17) yields identical limits to those based on the multiple image fraction (dashed light lines in figures 1 and 2) as computed in the usual way from equation (4). Moreover, the magnification distribution for singly imaged sources is naturally normalised within the formalism, and hence the magnification bias of singly imaged sources is directly included in the calculation.

## 4 Two lensed $`z6`$ quasars

While none of the six $`z6`$ quasars discovered by Fan et al. (2001b,2003) in the Sloan Digital Sky Survey data are multiply imaged, two have close alignment with a foreground galaxy, implying that they are moderately magnified. We therefore refer to these quasars as lensed, though neither is multiply imaged.

* SDSS J1044-0125 at $`z=5.74`$: Shioya et al. (2003) have reported a faint foreground galaxy with $`z1.52.5`$ at a separation of $`\theta =1.9^{\prime \prime }`$. Shioya et al. (2003) estimate the velocity dispersion to be $`\sigma 140280`$km/sec for this redshift range. A second image would be detectable if $`\sigma 220`$km/sec. For a SIS galaxy this implies that the magnification of the image could be as high as $`\mu =2`$, and that the image is located at $`x210`$ ER (where $`x10`$ER corresponds to a $`\sigma 140`$km/sec SIS at $`z2.5`$).
* SDSS J1148-5251 at $`z=6.37`$: This is one of two quasars known to exhibit Gunn-Peterson absorption troughs. White, Becker, Fan & Strauss (2003) present spectra showing emission features in the Ly $`\beta `$ trough which they interpret as being Ly $`\alpha `$ emission from a foreground galaxy at $`z4.9`$. This is a likely scenario if gravitational lensing is important and was predicted by Wyithe & Loeb (2002b). Indeed, the smaller than expected Stromgren sphere implies that this quasar has been magnified by the intervening galaxy (White et al. 2003). While the alignment is presumably high, we do not know the degree of alignment with the foreground galaxy, and hence we cannot know $`x`$, though is the quasar is not multiply imaged it must be larger than 2, and is probably comparable to the case SDSS J1044-0125. White et al. (2003) do however provide an estimate of $`\sigma 250`$km/sec for the foreground galaxy from the velocity structure seen in the C IV absorption system. For an SIS at $`z=4.9`$ lensing a source at $`z=6.37`$ with a separation $`\theta `$ we have

$$x=4\left(\frac{\theta }{1^{\prime \prime }}\right)\left(\frac{\sigma }{250\text{km/sec}}\right)^2\text{ER}.$$
(19)

## 5 improved limits on $`\beta `$ from close alignments with foreground galaxies

The Bayesian approach (§ 3) to computing limits on $`\beta `$ from the multiple image fraction allows us to include the a posteriori information on the alignments of the two lensed the $`z6`$ quasars.

For illustration, we begin with a hypothetical sample of 6 quasars with $`\mathrm{\Delta }M=4`$, $`\alpha =1.75`$ and (lensed) impact parameter $`x`$. The relative likelihoods for different values of $`\beta `$ given the lack of multiple images are $`\left[p_{\mathrm{single}}(\beta |x)\right]^{N_\mathrm{q}}`$ as specified in equation (15) where $`N_\mathrm{q}=6`$. These likelihoods are shown in the left hand panel of figure 1 for values of $`x`$ ranging from 5 to 20 (dark lines). Smaller values of $`\beta `$ are strongly disfavoured, particularly if $`x`$ is not too large. Next we find find the joint likelihood function and hence a posterior probability distribution for $`\beta `$ given a common impact parameter $`x`$ for six quasars that are not multiply imaged

$$\frac{dP}{d\beta }|_x=N[p_{\mathrm{single}}]^{N_\mathrm{q}}\mathrm{exp}\left(\frac{1}{2}\left(\frac{\beta \overline{\beta }}{\mathrm{\Delta }\beta }\right)^2\right)\frac{dP_{\mathrm{prior}}}{d\beta }.$$
(20)
The resulting posterior differential and cumulative probability distributions are shown in the central and right hand panels of figure 1 (dark lines). The most likely value is near $`\beta 2`$ and the cumulative distributions suggest that $`\beta 3.12.7`$ at the 90% level where the systematic dependence is on $`x`$. Thus any additional information about close alignments produces constraints that may significantly tighten the lower limits on the slope of the $`z6`$ quasar LF, with improvements in the limit that are greater than 0.5 units in $`\beta `$ for cases where the alignment is high.

The above example suggests strong dependence of the limits derived for $`\beta `$ on the value of $`x`$. For the two lensed SDSS quasars discussed in § 4, there are observational limits on $`x`$ in the form

$$\frac{dP}{dx}=L_x\frac{dP_{\mathrm{prior}}}{dx},$$
(21)
where $`L_x`$ is the likelihood for $`x`$ given the observations of $`\sigma `$ and $`z`$ for the lens galaxy, and $`\frac{dP_{\mathrm{prior}}}{dx}`$ is the prior probability for $`x`$. Given the relation $`x=x(\sigma ,z)`$ the likelihood $`L_x`$ is

$$L_x=L_\sigma L_z\frac{x}{\sigma }\frac{x}{z}\frac{dP_{\mathrm{prior}}}{d\sigma }\frac{dP_{\mathrm{prior}}}{dz}.$$
(22)
For SDSS J1044-0215 Shioya et al. (2003) find the majority of the dependence in the likelihoods for $`\sigma `$ and $`z`$ to be systematic, while there is no information in this regard for SDSS J1044-0215 (White et al. 2003). We assume flat distributions $`\frac{dP}{dx}`$ with limits of $`2<x<10`$ for the two lensed quasars. For the other 4 quasars we assume the prior probability distribution for $`x`$ (equation 18).

The posterior differential probability distribution for $`\beta `$ then becomes

$`{\displaystyle \frac{dP}{d\beta }}`$ $`=`$ $`N\mathrm{exp}\left({\displaystyle \frac{1}{2}}\left({\displaystyle \frac{\beta \overline{\beta }}{\mathrm{\Delta }\beta }}\right)^2\right){\displaystyle \frac{dP_{\mathrm{prior}}}{d\beta }}`$ (23)
$`\times `$ $`\left[{\displaystyle \underset{i=1,2}{}}\left({\displaystyle _2^{10}}𝑑x_ip_{\mathrm{single}}(\beta |x,L_i)\right)\right]`$
$`\times `$ $`\left[{\displaystyle \underset{i=3,6}{}}\left({\displaystyle _2^{\mathrm{}}}𝑑xp_{\mathrm{single}}(\beta |x,L_i){\displaystyle \frac{dP_{\mathrm{prior}}}{dx}}\right)\right].`$
The resulting cumulative probability distributions are shown in figure 2. The limits on $`\beta `$ are significantly tighter than those obtained from the distribution of luminosities alone except in cases where $`\alpha `$ is large (shallow faint end slope) and $`\mathrm{\Delta }M`$ is small (so that magnifications associated with multiple images tend to draw quasars with $`\mathrm{\Delta }M<0`$ into the sample). In addition the limits are tighter than those obtained through consideration of the luminosity distribution and multiple image fraction (dashed light lines). For a constant co-moving density of lens galaxies and $`\mathrm{\Delta }M2`$ we find $`\beta 3.13.0`$, while for $`\mathrm{\Delta }M4`$ we obtain $`\beta 3.02.9`$ (both with 90% confidence). Thus the tightest limits come from the inclusion of the a posteriori information that two of the quasars have close alignment with foreground galaxies.

A second important point regarding the limit on $`\beta `$ provided by equation (23) is that unlike the limit from the multiple image fraction, it is nearly independent of the value of $`\tau `$. To demonstrate this independence we have computed constraints on $`\beta `$ (shown in the lower two panels of figure 2) that assume a dependence in the co-moving density of galaxies of $`(1+z)^\gamma `$ where $`\gamma =2`$ (resulting in $`\tau =0.0004`$). These limits may be compared with results that assume a constant co-moving density of lens galaxies (in the upper two panels of figure 2). The limits obtained from the multiple image fraction are much weaker if $`\gamma =2`$. On the other hand the limits that use a posteriori observations of the quasar-lens galaxy alignment are quite insensitive to $`\gamma `$. The reason is that the role of $`\tau `$ is replaced by $`\frac{dP}{dx}`$ for the two quasars which provide the largest contribution to the likelihood change between large and small values of $`\beta `$.

### 5.1 a posteriori choice of statistic

We have computed limits on the value of $`\beta `$ using two lensing based constraints, and a posteriori chosen the better one. This practice becomes unfair if a large number of different constraints are available where each produces a different limit. In the situation described we have two different lensing constraints. However the second constraint utilises additional rather than different information. Thus we are justified in choosing it a posteriori.

## 6 Limits on galaxy evolution from lens galaxy redshifts

As noted by White et al. (2003) the candidate lens galaxy in the system SDSS J1148-5251 is found at an improbably high redshift, which could provide an argument against the lens hypothesis in this system. This is quantified in figure 3, where the upper curves show the differential cross-section for a source at $`z=6.37`$ lensed by a constant co-moving density of galaxies (left panel), and the fraction of the total cross-section that is found at a redshift larger than $`z`$ (right panel). We see that the prior probability of finding a lens at $`z4.94`$ among two lensed quasars is only $`0.01`$. However the probability of finding a lens at high redshift is an a-posteriori statistic, i.e. we have chosen one of a possible number of a-priori unlikely events after the observation has already been made. Moreover the selection of lenses within the sample is not uniform in redshift. In particular, since the galaxy in front of SDSS J1148-5251 was identified spectroscopically via its Ly-$`\alpha `$ emission line, it would be more easily identified at high redshift. In addition, the cross-sections plotted in figure 3 refer to multiple imaging so that the low probability for a high redshift lens results from the small size of the Einstein Ring radius at high redshift. In contrast, the redshift distribution of galaxies (with velocity dispersions larger than $`\sigma `$) that lie within some specified angular separation along the line of sight to a back-ground quasar implies that there would be 1 chance in $`3`$ of finding such an aligned galaxy at $`z>4.96`$. The role of magnification bias in making lensing a likely scenario probably results in an a-posteriori probability of finding a high redshift fore-ground galaxy that lies somewhere between these two extremes.

The redshift distribution of gravitational lenses may be used to constrain evolution in the lens galaxy population (Kochanek 1992; Ofek, Rix & Maoz 2003). While the absolute probability of having observed a high redshift lens is difficult to quantify, we may more easily discuss the relative likelihoods of observing a high redshift lens as a function of lens galaxy population. For definiteness we consider two parameterisations for the evolution of the lens galaxy population: firstly evolution of the characteristic velocity $`\sigma _{}(z)=\sigma _{}(1+z)^\gamma `$, and secondly evolution of the characteristic density $`\varphi _{}=\varphi _{}(1+z)^\gamma `$, which we term velocity and density evolution respectively. The left hand panels of figure 3 demonstrate the effect on the lens cross-section of varying $`\gamma `$. Values of $`\gamma `$ that differ from 0 (constant co-moving evolution) result in a lens population that is truncated at high redshift. This effect may also be seen in the central panels of figure 3 where we have plotted the fraction of cross-section at redshifts larger than $`z`$. Values of $`\gamma >0`$ are disfavoured by the existence of a lensing galaxy at $`z=4.96`$ among a sample of only two lens galaxies.

To quantify this statement we construct a likelihood function for $`\gamma `$ from the product of the normalised probabilities for the lens redshifts. The likelihood should include constants ($`s_i`$) to account for the relative detectabilities of the two lenses (the two galaxies were discovered separately via different techniques), though the limits on galaxy evolution are independent of these since the constants are independent of the evolution. The likelihood function is

$$L_\gamma =\underset{i=1}{\overset{2}{}}s_i_2^{\mathrm{}}𝑑x\frac{dP_i}{dx}(x1)^2\frac{d\tau _i}{dz}B_i\underset{i=1}{\overset{2}{}}\frac{d\tau _i}{dz},$$
(24)
where the $`\frac{d\tau }{dz}`$ are differential cross-sections evaluated at the lens and source redshifts, the integrals over the distributions $`\frac{dP_i}{dx}(x1)^2`$ account for the relative alignments of the quasar and galaxy, and the $`B_i`$ are the magnification biases. The relative likelihood is dependent only on the product of the differential cross-sections. Note that the likelihood (equation 24) is only applicable if the magnification bias has aided in selection of the quasar, so that the source may be considered lensed. This is a caveat to the constraints imposed on $`\gamma `$ in this section. However the quasar does appear to be magnified, as evidenced by the smaller than expected Stromgren sphere (White et al. 2003).

In the right hand panels of figure 3 we plot the posterior cumulative probability for $`\gamma `$

$$P(<\gamma )=^\gamma NL_\gamma ^{}\frac{dP_{\mathrm{prior}}}{d\gamma ^{}},$$
(25)
assuming a flat prior probability distribution for $`\gamma `$ at values greater than $`0`$. By excluding the possibility of $`\gamma <0`$ we are assuming that the lens galaxy population increases monotonically in time as expected in hierarchical merging scenarios. This choice also leads to more conservative limits on $`\gamma `$. We find $`\gamma 0.4`$ and $`\gamma 1.6`$ at the 90% level assuming velocity and density evolution respectively. The possible presence of a lens galaxy at such a high redshift therefore offers an opportunity to constrain the (mass selected) co-moving density of massive galaxies to be close to constant out to high redshifts. This result is consistent with the study of Ofek, Rix & Maoz (2003) who performed a detailed study on a large sample (15) of multiple image lenses at $`z12`$.

A constant co-moving density of lens galaxies out to $`z5`$ may not be surprising in light of the Press-Schechter (1976) prediction for the velocity function of dark-matter halos (number per cubic comoving Mpc per unit velocity). Taking the circular velocity $`v_{\mathrm{vir}}`$ to equal the virial velocity of an SIS dark matter halo with mass $`M`$ (Barkana & Loeb 2001), and a velocity dispersion $`\sigma =v_{\mathrm{vir}}/\sqrt{2}`$ we find

$$\frac{dn}{d\sigma }=M\frac{dn}{dM}\frac{\sqrt{2}}{3v_{\mathrm{vir}}},$$
(26)
where $`\frac{dn}{dM}`$ is the Press-Schechter (1976) mass function. The resulting velocity function of dark matter halos is plotted in figure 4 at a series of redshifts. Note that near velocity dispersions of $`\sigma 200`$km sec<sup>-1</sup>, which dominate the lens cross-section, there is little evolution in $`\frac{dn}{d\sigma }`$ from $`z1`$ to $`z5`$. If massive galaxies occupied dark matter halos in the past as they do today, we would therefore expect little evolution in the lens population, even out to large redshifts.

## 7 What is the multiple image lensing rate for $`z6`$ quasars?

We may combine the information obtained for $`\beta `$ and $`\gamma `$ and estimate the expected multiple imaging rate for $`z6`$ quasars. Since our earlier limits on $`\beta `$ are independent of $`\tau `$, while the limits derived for $`\gamma `$ are insensitive to the magnification bias, the constraints placed on $`\beta `$ and $`\gamma `$ in § 5 and § 6 are independent. In figure 5 we have plotted the resulting joint probability function (dark contours); the contours shown are at 61, 26, 14 and 3.6% of the peak value, corresponding to the 1, 2, 3 and 4 $`\sigma `$ levels of a Gaussian distribution. The upper two and lower two rows of figure 5 correspond to velocity and density evolution respectively and in each case results are shown for a double powerlaw LF with various values of $`\alpha `$ and $`\mathrm{\Delta }M`$. The figure shows that the preferred values are found near near $`\gamma =0`$ and $`\beta =2.1`$. We also show contours of the multiple image fraction (light contours); the solid, dashed, dot-dashed and dotted contours correspond to lens fractions of 0.1, 0.03, 0.01 and 0.003. We find that the multiple image fraction should be $`13\%`$. This value is lower than previous estimates due to constraints on the shallow luminosity function, and unfortunately implies that a $`z6`$ lens may not be found among the complete sample of SDSS $`z6`$ quasars.

## 8 What if SDSS J1148-5251 were multiply imaged?

While current high resolution imaging suggests that all of the $`z6`$ quasars are point sources, White et al. (2003) note the possibility that SDSS J1148-5251 is multiply imaged cannot be ruled out by current observations due to the small angular diameter of the Einstein ring in this system; $`\mathrm{\Delta }\theta 0.3^{\prime \prime }`$ for a $`\sigma =250`$km/sec galaxy at $`z=4.94`$ lensing a quasar at $`z=6.37`$. For comparison the $`z6`$ quasars have been imaged at a resolution of $`0.4^{\prime \prime }`$ (Fan et al. 2003). It is therefore possible that SDSS J1148-5251 is multiply imaged but appears as a point source, though we note that this is an unlikely scenario since with $`4^{\prime \prime }`$ seeing, a double with a $`3^{\prime \prime }`$ separation should be recognisable if the flux ratio is smaller than $`10:1`$ (Chris Kochanek 2003, private communication). Multiple imaging of SDSS J1148-5251 would have important implications for the study of the $`z6`$ LF, invalidating the constraints on $`\beta `$ obtained in § 5. We have therefore computed the limits imposed on $`\beta `$ by the observation of one multiply imaged source (SDSS J1148-5251) and one singly imaged source with high alignment (SDSS J1044-0125) among a sample of six quasars at $`z6`$.

The likelihood function for this scenario may be written as the product of the probability that a source is lensed with the likelihood that the remainder are singly imaged. The observation that one of the quasars is multiply imaged constrains large values of $`\beta `$, while the observation of high alignment without multiple imaging limits small values of $`\beta `$ as discussed in previous sections. The maximum of the combined likelihood function lies in the lower end of the range specified by the luminosity distribution (Fan et al. 2003). The posterior probability function for the sample of $`z6`$ quasars may be written

$`{\displaystyle \frac{dP}{d\beta }}`$ $`=`$ $`N\mathrm{exp}\left({\displaystyle \frac{1}{2}}\left({\displaystyle \frac{\beta \overline{\beta }}{\mathrm{\Delta }\beta }}\right)^2\right){\displaystyle \frac{dP_{\mathrm{prior}}}{d\beta }}`$ (28)
$`\times p_{\mathrm{lens}}(\beta ,L_1)\times \left[{\displaystyle _2^{10}}𝑑xp_{\mathrm{single}}(\beta |x,L_2)\right]`$
$`\times \left[{\displaystyle \underset{i=3}{\overset{6}{}}}\left({\displaystyle _2^{\mathrm{}}}𝑑xp_{\mathrm{single}}(\beta |x,L_i){\displaystyle \frac{dP_{\mathrm{prior}}}{dx}}\right)\right]`$
and is plotted in figure 6 for various values of $`\mathrm{\Delta }M`$ and $`\alpha `$. We find the preferred value in this case to be $`\beta 3`$. As noted in § 2.1 it is surprising that one quasar in the sample would be lensed, but not at a highly significant level. The observation of one lensed lensed quasar therefore prefers smaller values of $`\beta `$, for which the magnification bias is larger, and also slightly tightens the allowed range for $`\beta `$. In summary, if SDSS J1148-5251 were multiply imaged the preferred value for the slope would be $`\beta 3`$, which is ruled out at the 90% level if the quasar is singly imaged but with a high alignment. From figure 5 we see that $`\beta 3`$ implies a lens fraction of $`310\%`$ rather than $`0.31\%`$, which is more consistent with previous estimates. This underlines the importance of determining whether this quasar is multiply imaged or merely lensed.

## 9 conclusion

From their analysis of the luminosity distribution of quasars at $`z6`$ Fan et al. (2003) determined a slope for the quasar LF of $`\beta 3.3`$ (90%). This slope is consistent with the value found for the slope of the LF at $`z4.3`$, but is not consistent with the slope of the LF of bright quasars at $`z3`$. It is also possible to constrain the slope of the LF using the fraction of multiply imaged lensed quasars. Fan et al. (2003) computed the probability of obtaining a lens fraction of zero as a function of $`\beta `$ and found that the constraints were similar to those of the LF. We have performed a Bayesian analysis including both of these (nearly) independent constraints, yielding the result that at 90% confidence $`\beta 3.33.0`$ provided that the quasars are at least 2 magnitudes brighter than the unknown position of a break in a double powerlaw LF. The systematic dependence in the constraint is due to the unknown slope of the LF at fainter luminosities, the luminosity of the break and the uncertain evolution in the lens galaxy population.

While inclusion of constraints from the multiple image fraction somewhat improves the limits on $`\beta `$, we have shown that the additional information from observations that neither of the two quasars which lie near to the line of sight to foreground galaxies (and which are therefore lensed) are multiply imaged provides a stronger lensing based constraint on the slope $`\beta `$. We find that for a double powerlaw LF $`\beta 3.12.9`$ with 90% confidence. Unlike the constraint that uses only the multiple image fraction, this limit is nearly independent of evolution in the lens population, and adds further evidence of a trend to shallower LF slopes at large redshifts.

We also find that the existence of a lens galaxy at $`z5`$ in a sample of two lenses constrains the evolution in the massive galaxy population to be close to that of constant co-moving density (provided that the quasar behind the $`z5`$ galaxy is subject to magnification bias). This lack of evolution in the lens population is consistent with the lack of redshift evolution in the velocity function of dark-matter halos (for velocity dispersions near 200km sec<sup>-1</sup>) as predicted by the Press-Schechter formalism.

Finally the constraints on the quasar luminosity function and lens population have been used to compute an improved estimate for the expected $`z6`$ multiple image lens fraction of $`13\%`$. This value is lower than previous estimates due to the tight constraints on the slope of the LF.

## Acknowledgements

I would like to thank Chris Kochanek for helpful comments during the course of this work.


text:
# Single-Spin Beam Asymmetry in Semi-Exclusive Deep-Inelastic Electroproduction

## Abstract

Recent measurements from Jefferson Lab show significant beam single spin asymmetries in deep inelastic scattering. The asymmetry is due to interference of longitudinal and transverse photoabsorption amplitudes which have different phases induced by the final-state interaction between the struck quark and the target spectators. We developed a dynamical model for a single-spin beam asymmetry in deep-inelastic scattering. Our results are consistent with the experimentally observed magnitude of this effect. We conclude that similar mechanisms involving quark orbital angular momentum (‘Sivers effect’) are responsible for both target and beam single-spin asymmetries.

Single–Spin Asymmetries (SSA) in semi–exclusive electroproduction reactions provide important new information about the nucleon spin structure sivers91 ; collins92 ; jjmt . Recent successful measurements of SSA in leptoproduction hermes\_ssa ; CLAS02  and anticipated experiments at Jefferson Lab after the 12-GeV energy upgrate stimulate a lot of interest in this subject.

One can relate SSA to the products of parton distribution and fragmentation functions, see Ref.baronereview  for review. In the quark–parton model, the origins of SSA can be attributed both to the time–reversal odd quark distribution, often called the ‘Sivers effect’ sivers91 , and to the time–reversal odd fragmentation also known as ‘Collins effect’ collins92 . Recently, Brodsky and collaborators bhs  used a model to demonstate that the target SSA can be generated at a leading twist level without time–reversal odd fragmentation. It led to the revision collins ; bjy  of time–reversal arguments, validating phenomenology with time–reversal odd parton distributions.

It has been known for a long time that single–spin polarization observables in particle scattering are caused by the spin–orbit interaction mott ; nucl . SSA in deep–inelastic scattering should become a sensitive probe of the orbital angular momentum (OAM) of the participating quarks since SSA would simply vanish without contributions from both spin and OAM. It was demonstrated bhs  that indeed OAM of quarks is crucial for describing target SSA.

In this work, we extend the calculation of SSA to the case of a polarized electron beam using the model that Ref.bhs  used for the polarized target case. The model is reminiscent of model decomposition of the proton studied earlier in Ref.ezawa-drell-lee . We demonstrate that a mechanism similar to Ref.bhs  is also generating beam SSA of a magnitude compatible with experimental data. The considered mechanism does not require the Collins fragmentation function and therefore does not allow relation of the measured asymmetry to the quark distribution $`e(x)`$ introduced by Jaffe and Ji jjmt . In this respect, our result for the beam SSA differs from analyses of Refs.leveltmulders ; efremovgoeke .

We calculate the beam SSA $`A_{LU}^{\mathrm{sin}\varphi }`$ measured in the experiment and defined as

$$d\sigma (1+h_eA_{LU}^{\mathrm{sin}\varphi }\mathrm{sin}\varphi ),$$
(1)
where $`h_e=\pm 1`$ indicates the electron beam helicity, and $`\varphi `$ is the angle between the planes of electron scattering and hadron production. Following Ref.bhs , we describe the deep–inelastic electroproduction process with a model of a spin 1/2 proton of mass $`M`$ which is made from charged spin–1/2 and spin–0 constituents of masses $`m`$ and $`\lambda `$, respectively (Fig.1). During the scattering, the virtual photon is absorbed by an active spin–1/2 parton, while the scalar diquark is used to describe the spectator system. The momentum of the knocked–out quark is set equal to the observed hadronic momentum. The proton–quark–diquark vertex function should be a solution to a bound–state equation, which we here simplify as a scalar vertex with a coupling parameter $`g`$.

The leading–order diagram Fig.1a alone produces zero asymmetry $`A_{LU}^{\mathrm{sin}\varphi }`$, since its contribution to the scattering amplitude is purely real. The one–loop next–to–leading–order correction Fig.1b needs to be included in order to generate nonzero beam and target SSA. The loop diagram Fig.1b contribution is given by the formula

$$J_\mu ^{NLO}=ie_qgC_F4\pi \alpha _s\frac{d^4k}{(2\pi )^4}\frac{\overline{u}(q+r)\gamma _\nu (\overline{)}q+\overline{)}k+m)\gamma _\mu (\overline{)}k+m)(g_{\nu \tau })(2pkr)_\tau }{D_1D_2D_3D_4},$$
(2)
where $`C_F=4/3`$ and $`D_i`$ are denominators that come from the four propagators.

We calculate the contribution of the Fig.1 mechanism to the anti–symmetric part of the hadronic tensor defined by the imaginary part of the intereference between Fig.1a and Fig.1b amplitudes. Cutkosky rules are used to calculate the absorptive part of the loop contribution, replacing two of the denominators $`D_i`$ in Eq.(2) by the corresponding delta functions as follows:

$$\frac{1}{((q+k)^2m^2+iϵ)((pk)^2\lambda ^2+iϵ)}(2\pi i)^2\delta ((q+k)^2m^2)\delta ((pk)^2\lambda ^2),$$
(3)
where $`m`$ and $`\lambda `$ are the masses of the quark and diquark, respectively. As a result, the 4–dimensional integration over the loop momentum in Eq.(2) is reduced to a 2–dimensional angular integration, making the result safe from ultraviolet divergence. The beam asymmetry is generated due to interference between absorption of longitudinal and transverse virtual photons (LT-interference) which acquire different phases as a result of final–state gluon exchange. In the scaling limit, the asymmetry is proportional to $`r_T/\sqrt{Q^2}`$, where $`r_T`$ is a transverse component of the final quark momentum. The factor $`1/\sqrt{Q^2}`$ points to the twist–three origin of this observable. The calculations were also performed in the light–cone gauge, producing the same results.

The amplitude Fig.1b is infrared–divergent as the momentum of the exchanged gluon approaches zero. However, the infrared–divergent terms exactly cancel at the level of observable asymmetry. The same spin–dependent mechanism reduces relative contribution of soft gluons to the integral Eq.(2). We also note that we neglected the loop correction to the unpolarized cross section, otherwise we would have to add real gluon bremsstrahlung or introduce an infared cut–off to make the theory infrared–safe.

The results of our model for beam SSA is shown in Fig.2 for different beam energies as a function of Bjorken $`x_{Bj}`$. We used the values of the parameters $`m=`$ 0.3 GeV, $`\lambda `$= 0.8 GeV and $`\alpha _s`$= 0.3. The value of the kinematic variable $`y=`$ 0.7 and beam energies of 4.25 and 6 GeV are chosen to match the experimental conditions CLAS02 . One can see that the prediction for beam SSA is consistent with its observed magnitude. For higher beam energies of 27 GeV relevant to the HERMES experiment hermes\_ssa  and fixed values of the kinematic variables $`x_{Bj},y`$, the beam asymmetry is supressed by about a factor of two, following the $`1/\sqrt{Q^2}`$ dependence. Our model predicts that the asymmetry changes its sign around $`x=`$ 0.5, where no experimental data are available. Thus extension of experimental measurements to higher values of $`x`$ can test validity of the proposed mechanism. Note that beam SSA calculated in a simplified model with scalar partons is approximately constant as a function of $`x_{Bj}`$ and equals about +0.08 for the 6-GeV kinematics of Fig.2.

Let us discuss the relation of the obtained results to the target SSA calculations of Ref.bhs . Target SSA are due to the T-odd observable $`\stackrel{}{S}_p\stackrel{}{q}\times \stackrel{}{p}_h`$ made from the proton polarization vector $`\stackrel{}{S}_p`$, the photon momentum $`\stackrel{}{q}`$, and the outgoing hadron’s momentum $`\stackrel{}{p}_h`$. The essential ingredient of the model bhs  is quark OAM in the initial nucleon state, so that the struck parton may have helicities both parallel and antiparallel to the nucleon helicity. Interference between transitions to different helicity states, supplemented by the phase difference due to final gluon exchange is the mechanism that generates the target SSA. Very similarly, beam SSA is due to the correlation $`\stackrel{}{S}_\gamma \stackrel{}{q}\times \stackrel{}{p}_h`$, where $`\stackrel{}{S}_\gamma `$ is the virtual photon polarization. The beam SSA arises due to intereference between longitudinal and transverse photoabsorption amplitudes, with the transverse momentum and phase differences generated by the same gluon–exchange mechanism as in the polarized target case. The quark OAM now unambiguosly enters at the photon–parton level and the result does not depend on the quark OAM contribution to the nucleon light–cone wave function. Therefore, measurements of beam SSA are needed for adequate interpretation of target SSA measurements within a parton model. By analogy with the target asymmetry and in absence of T–odd fragmentation, we can view the considered mechanism as a photon Sivers effect, emphasizing the importance of this effect in polarized deep–inelastic scattering.

We acknowledge useful discussions with H. Avakian, S. Brodsky and D.S. Hwang. This work was supported by the US Department of Energy under contract DE-AC05-84ER40150. In addition, C.E.C. thanks NSF for support under grant PHY-0245056.


text:
# Non-linear resonance in the accretion disk of a millisecond pulsar

## 1 Introduction

A millisecond pulsar that accretes matter from a binary stellar companion provides a unique probe of the physics of accretion flow a few Schwarzschild radii away from the neutron star. In contrast with the case of the slowly rotating steady X-ray pulsars, the magnetic field of the (X-ray transient) millisecond pulsar is too weak to prevent the accreting gas from orbiting the neutron star very close to its surface. Hence, in such a pulsar oscillations in the inner parts of the disk may be observed at their characteristic frequency of several hundred Hertz. The pulsar is expected to disturb the accretion disk at its spin frequency. That it does, is demonstrated by the discovery in the 2.5 ms accreting pulsar of a frequency difference between the two QPOs equal to one-half the pulsar spin rate.

We point out that the two QPO frequencies themselves may be commensurable with the spin frequency and with each other. The two simultaneously observed QPO frequencies are close to 500 Hz and 700 Hz, i.e., they may be in a 5:7 ratio. Twin QPOs have been proposed to be a manifestation of non-linear resonance in relativistic accretion disks (Kluźniak and Abramowicz 2001), and the resulting suggestion that twin QPO frequencies should be in the ratio of small integers has subsequently been confirmed for a number of black hole transients (Abramowicz and Kluźniak 2001; Remillard et al. 2002; Kluźniak and Abramowicz 2002; McClintock and Remillard 2003) and for the bright steady source Sco X-1 (Abramowicz et al. 2003a).

The two QPOs discovered in the 2.5 ms accreting pulsar (Wijnands et al. 2003) are similar to QPOs previously detected in many unpulsed X-ray sources, in which an accretion disk is thought to feed matter and angular momentum to the neutron star or black hole (van der Klis et al. 2000; Lewin et al. 1995). Such QPOs in black hole systems were intrerpreted as oscillations of the accretion disk (Wagoner 1999; Kato 2001), but a direct proof of this was lacking, and for neutron star systems other models had been proposed (van der Klis et al. 2000). SAX J1808.4-3658 is the first kHz QPO source in which the rotational period of the neutron star has been unambiguously determined, and the first in which the frequency difference of the twin QPOs is known to be a subharmonic of the stellar rotation rate. For the first time, we have direct evidence that kHz QPOs are caused by disk oscillations.

## 2 Non-linear resonance and epicyclic frequencies

The presence of a sub-harmonic frequency is a clear signature of non-linear resonance (Nayfeh and Mook 1979). Theodor von Karman observed in 1940 that high-frequency vibrations of an airplane engine excite lower-frequency resonances in the airframe (von Karman 1940; Nayfeh and Mook 1979). A tragic example is known when a subharmonic resonance was excited in an airplane wing, which in turn excited a subharmonic resonance in the rudder, at 1/2 the wing eigenfrequency (Lefschetz 1956; Nayfeh and Mook 1979). In the present case, the 2.5 ms pulsar SAX J1808.4-3658 plays the role of the engine, while the accretion disk (whose response is observed as the quasi-periodic oscillations) is the airframe with its own set of known eigenfrequencies. The fact that the QPO frequency difference coincides with one half the known pulsar spin frequency,

$$(694\pm 4)\mathrm{Hz}(499\pm 4)\mathrm{Hz}=(195\pm 6)\mathrm{Hz}\frac{1}{2}(401\mathrm{Hz}),$$

follows directly from the nature of the system. The structure of the accretion disk is determined by the non-linear equations of hydrodynamics, and the effective gravitational potential (in which the accreting gas moves) is not harmonic, so a non-linear response is expected. Below, we discuss in detail how this comes about for an accretion disk in the space-time of a rotating neutron-star.

The motion of a test particle in nearly circular orbits close to the equatorial plane can be decomposed into three components, circular planar motion at the orbital frequency $`\mathrm{\Omega }2\pi \nu _{\mathrm{orb}}`$, harmonic radial motion at the radial epicyclic frequency $`\kappa 2\pi \nu _\mathrm{r}`$, and harmonic vertical motion at the meridional epicyclic frequency $`\zeta 2\pi \nu _{\mathrm{vert}}`$. In Newtonian theory of spherically symmetric gravitating bodies the three frequencies coincide, but in Einstein’s gravity $`\kappa <\mathrm{\Omega }`$. (The relativistic precession of the perihelion of Mercury occurs at the rate $`\mathrm{\Omega }\kappa `$.) The same three frequencies are important in a discussion of fluid motion about a gravitating body in general relativity (Kato 1998). The accretion disk is a body of hot gas that is supported against infall primarily by rotation and is nearly in hydrostatic equilibrium (Lewin et al. 1995; Kato 1998). Theories of accretion (Shakura and Sunyaev 1973; Jaroszyński et al. 1980) admit solutions in which the disk thickness is much smaller than its radial extent, but also solutions in which the disk has the geometry of a torus. Like other extended bodies in equilibrium, the disk is capable of motion in a variety of modes. In the linear regime, small oscillations of geometrically thin accretion disks, as well as waves in this body of fluid, have been extensively studied in the Kerr metric (Wagoner 1999; Kato 2001). The radial vibrations of two-dimensional models of geometrically thick disks (accretion torii) have been investigated numerically in the Schwarzschild metric (Rezzolla et al. 2003). Qualitatively, these black-hole metrics are similar to the metrics of neutron stars.

We have studied numerically the response to a transient external perturbation of an ideal-gas torus that is initially in equilibrium rotation about a (Schwarzschild) black hole or neutron star. We find that the torus performs harmonic oscillations both in the radial and vertical directions, even if the impulse imparted to it at the beginning of the computation is purely radial. Further, the vertical oscillations have variable amplitude. Both these effects speak of a non-linearity in the system.The vertical and radial displacements of the center of the torus (defined as the point of maximum pressure in a meridional cross-section) are shown as a function of time in Fig. 1. The radial displacement has been rescaled to fit the graph—in the simulation described here, the amplitude of radial motion is actually larger than the amplitude of vertical motions.

For the numerical simulation (Fig. 1), we use a Newtonian 2-d smooth particle hydrodynamics (Monaghan 1992, SPH) code (Lee and Ramirez-Ruiz 2002) simulating ideal gas with adiabatic index 4/3 in a pseudo-potential (Kluźniak and Lee 2002) $`\psi _{KL}(r)=[1\mathrm{exp}(r_{\mathrm{ms}}/r)]GM/r_{\mathrm{ms}}`$, which reproduces the Schwarzschild ratios of the orbital and epicyclic frequencies, $`\zeta /\mathrm{\Omega }=1`$, $`(\kappa /\zeta )^2=1r_{\mathrm{ms}}/r`$. Here, $`r_{\mathrm{ms}}=6M(G/c^2)`$ is the radius of the marginally stable orbit. A torus in equilibrium, with $`r_0=12.25M(G/c^2)`$, was perturbed at time $`t=0`$ by imparting to it a radial velocity field with magnitude proportional to $`\sqrt{r_{\mathrm{ms}}/r}`$. The frequencies of the resulting motions are found to be inversely proportional to the central mass $`M`$, whose value was adjusted to match the observed QPO frequencies.

The two frequencies observed in this numerical simulation coincide, within errors, with the known epicyclic frequencies, $`\kappa `$ and $`\zeta `$, at the equilibrium radial position, $`r_0`$, of the locus of maximum pressure, $`r_c`$. We conclude that a slender torus responds to a radial external perturbation with oscillations occuring at frequencies equal to the two epicyclic frequencies, $`\kappa (r_0)`$ and $`\zeta (r_0)`$, at that radius where the pressure of the torus in equilibrium is highest. An important implication is that these frequencies can vary if the torus varies with time as it accumulates less or more of the matter flowing through it.

For the radial perturbation applied in the simulation, the amplitude of vertical motion is smaller than that of the radial motion. Presumably the ratio of the amplitudes depends on the perturbation applied. Quantifying this requires further study. A tilted magnetic dipole may perturb the disk quite strongly in the vertical direction. In the remainder of this section we assume that the perturbation is such that only one of the two motions is typically excited with a large amplitude, and that the other would not be as easily discernible in an X-ray observation, unless its initially much smaller amplitude were amplified by a resonance.

A pulsar is a rotating neutron star with a strong magnetic dipole not aligned with the rotation axis. As a result, in the 2.5 ms accreting pulsar (Wijnands and van der Klis 1998) the accretion disk suffers a periodic disturbance at the spin frequency $`\nu _1=401`$Hz. Observations indicate the presence of a QPO varying in frequency between 280 and 750 Hz (Wijnands et al. 2003). We interpret this as one of the two epicyclic frequencies, $`\nu _0`$, and expect that a second QPO will be present at the second epicyclic frequency, when a resonance occurs at frequency $`\nu _{\mathrm{res}}`$ between the corresponding disk oscillation and the disturbing forces.

Twin QPOs have been observed in this source on only one occasion, and it is not clear to us whether the single variable frequency $`\nu _0`$ present in the remaining observations corresponds to the upper or the lower QPO, i.e., to $`\nu _{\mathrm{vert}}`$ or to $`\nu _r`$ (but see below). A simple and general reason why the difference is equal to one half the spin frequency, $`\nu _{\mathrm{res}}\nu _0=\nu _1/2`$, and not the spin frequency, $`\nu _1`$, can be given if $`\nu _0=\nu _r`$. Then $`\nu _{\mathrm{res}}=\nu _{\mathrm{vert}}`$ at a particular radius, and the resonance would occur in an oscillator without quadratic anharmonicity terms. The restoring force in the vertical motion of the torus is obtained by taking the $`z`$ derivative of the effective potential expanded about the equatorial plane. Because the potential is symmetric under reflection in this plane, this expansion has only even powers of the cylindrical co-ordinate $`z`$, $`V(z)=V_0+(\zeta z)^2/2+\beta z^4/4+`$… The vertical motion of the torus can then be described as that of an oscillator with (angular) eigen-frequency $`\zeta `$, and a cubic anharmonicity in leading order, $`\beta z^3`$.

Consider resonance in an accretion disk which can be described as a damped oscillator with a cubic anharmonicity. In view of the very strong disturbance by the pulsar, the resonance may be driven by a combination of the pulsar spin frequency $`\nu _1`$ and another, much weaker, disturbance at frequency $`\nu _2`$, which we take to be a harmonic (i.e., an integer multiple) of the QPO frequency $`\nu _0`$. In the present case the frequencies of the two harmonic driving forces satisfy $`\nu _2>\nu _1`$. The oscillator may be in resonance with both of these forces at the same time (Nayfeh and Mook 1979, simultaneous resonance) in one of three cases:

a) $`\nu _{\mathrm{res}}2\nu _1\pm \nu _2`$, or

b) $`\nu _{\mathrm{res}}2\nu _2\pm \nu _1`$, or

c) $`\nu _{\mathrm{res}}(\nu _2\pm \nu _1)/2`$.

The single occasion when two QPOs were detected in the 2.5 ms pulsar corresponds to case c), with $`\nu _2=2\nu _0`$, i.e., $`\nu _{\mathrm{res}}\nu _0\pm \nu _1/2`$. Since for neutron stars in general relativity $`\nu _{\mathrm{vert}}>\nu _r`$, the upper sign corresponds to the case $`\nu _{\mathrm{res}}=\nu _{\mathrm{vert}}`$, $`\nu _0=\nu _r`$, the lower sign to the case $`\nu _{\mathrm{res}}=\nu _r`$, $`\nu _0=\nu _{\mathrm{vert}}`$. Indeed, the two observed QPO frequencies, are in this relation to the spin frequency of 401 Hz, $`\nu _{\mathrm{upper}}\nu _{\mathrm{lower}}=694\mathrm{Hz}499\mathrm{Hz}(401\mathrm{Hz})/2=\nu _1/2`$.

We note that the two QPOs cannot be separated by the spin frequency in this scheme. A separation between the two QPO frequencies of twice the spin frequency, ($`\nu _{\mathrm{res}}\nu _0=802`$ Hz, case a with $`\nu _0=\nu _2`$) is in principle also possible in this scheme if resonance occurs at the very inner edge of the disk, at frequencies of about 650 Hz and 1450 Hz, but it does not seem likely that the bulk of the torus could reside so close to the neutron star.

We have numerically computed realistic models of rotating neutron stars and their exterior metrics. Fig. 2 presents a typical result consistent with observations of the two QPOs in the 2.5 ms pulsar. The vertical and radial epicyclic frequencies for a 1.22 solar mass neutron star, modeled with equation of state FPS, and rotating at 401 Hz are illustrated as a function of the radius. The case c) resonance discussed above occurs when pressure in the torus has a maximum at radius $`r_0=11.27M(G/c^2)=20.6`$km. The stellar mass and resonance radius are only indicative of the true values for this pulsar, as the frequencies for a flattened torus of large radial extent may differ from $`\kappa (r_0)`$ and $`\zeta (r_0)`$ by several per cent (Rezzolla et al. 2003).

Wijnands et al. (2003) give phenomenological arguments in favor of the frequency being the upper QPO. In this case, in addition to the resonance described below, a resonance at the combination frequency $`\nu _r=\nu _{\mathrm{vert}}\nu _1`$ could also occur in principle if the radial oscillator with its quadratic anharmonicity were perturturbed simultaneously at the spin frequency and the frequency of the vertical oscillation. This would happen at $`r_0=16`$km for the same model of the neutron star, when the upper QPO frequency has the (as yet unobserved value) 1015 Hz and the lower frequency is 615 Hz (dashed vertical line in Fig. 2).

To compute the frequencies in Figs. 2 and 3, we constructed numerical models of rotating neutron stars and their exterior metrics using the relativistic code of Stergioulas and Friedman (20 1995), which solves Einstein’s field equations for arbitrarily large rotation rate in an integral form. Details of the numerical method and extensive accuracy tests can be found in (Stergioulas 2003). The results presented here are not very sensitive to the choice of equation of state (EOS) of neutron-star matter.

## 3 Parametric resonance

So far we have been discussing a resonance directly caused by an external disturbance. But when the frequencies $`\kappa (r)`$ and $`\zeta (r)`$ happen to be in particular ratios, large-amplitude motions of the torus may be excited spontaneously. This mechanism, and specifically parametric resonance between the two modes of oscillation, has been suggested (Kluźniak and Abramowicz 2002) as the origin of the high-frequency QPOs in black hole systems, where the frequency ratios 5:3 and 3:2 have been noted (Abramowicz and Kluźniak 2001; Remillard et al. 2002; Kluźniak and Abramowicz 2002; McClintock and Remillard 2003), as well as in the candidate neutron-star system Sco X-1 (Abramowicz et al. 2003a, b). The presence of subharmonics in some of the black hole systems has also been discussed (Remillard et al. 2002; Kluźniak and Abramowicz 2003). On the other hand, in several X-ray bursters the frequency difference between the twin QPOs has been reported to be close to the suspected spin frequency (or one half of it) of the neutron star (van der Klis et al. 2000). It is interesting to note that in addition to being separated by one-half the pulsar spin frequency, the two QPO frequencies in SAX J1808.4-3658 are approximately in a 7:5=1.4 ratio. The system may thus be the missing link between the QPOs in black hole systems and X-ray bursters. It is known that if the frequency difference of two modes is accidentally close to that of an external disturbance, a resonance involving the three frequencies may occur (Kato 1974). In neutron star systems this could favor the appearance of strong QPOs with a frequency difference selected by the perturbing neutron-star rotation, while in black hole disks only the spontaneous resonance would occur.

Parametric resonance between two oscillations is possible in accretion disks because the coupling of modes is non-linear in hydrodynamics. For thin disks, this can lead to excitation of two modes whose frequency ratio is closer to 1.4 than to 1.5. In the Schwarschild metric, the damping or excitation of g-modes with azimuthal dependence $`\mathrm{exp}(im\varphi )`$ has been considered by Kato (Kato 2003), who found that an $`m=0`$ mode and an $`m=1`$ mode will undergo resonant amplification when interacting with a non-rotating one-armed ($`m=1`$) stationary warp, assumed to be present in the accretion disk. We assume that the same modes are unstable also in a gravitational field that is not spherically symmetric, and compute their frequencies for the metric of rotating neutron stars. The sum of the two frequencies turns out to be $`\nu _{\mathrm{orb}}`$, and if they are identified with the two observed QPO frequencies this constrains the mass of the 2.5 ms pulsar to be less than 1.5 solar mass. The two mode frequencies, $`\nu _{\mathrm{lower}}`$ and $`\nu _{\mathrm{upper}}`$, as well as $`\nu _{\mathrm{orb}}`$, $`\nu _r`$, and $`\nu _{\mathrm{orb}}\nu _r`$ are shown in Fig. 3 for a 1.49 solar mass neutron star, spinning at a frequency of 406 Hz. The resonance occurs in the region where $`\nu _{\mathrm{lower}}<\nu _r`$, with the largest modulation of X-ray flux expected when the two frequencies are evaluated at the edge of the resonance zone, i.e., when they are $`\nu _{\mathrm{lower}}\nu _r`$ and $`\nu _{\mathrm{upper}}\nu _{\mathrm{orb}}\nu _r`$.

The two unstable modes in Fig. 3 were taken, following (Kato 2003), to be $`m=1`$, $`n=1`$, and $`m=0`$, $`n=1`$, g-modes excited by parametric resonance, which is mediated by a warp assumed to be present in the thin accretion disk. For a non-rotating warp their frequencies (Kato 2003) would be $`2\nu _{\mathrm{orb}}\sqrt{2}\nu _{\mathrm{vert}}`$ and $`\sqrt{2}\nu _{\mathrm{vert}}\nu _{\mathrm{orb}}`$, respectively, with their ratio tending to $`\sqrt{2}=1.41`$ in the limit of non-rotating stars, or in the limit of large radial distance from the star, $`r\mathrm{}`$. The warp can be identified with a c-mode of disk oscillations and is expected to rotate at the Lense-Thirring angular frequency $`2\pi (\nu _{\mathrm{orb}}\nu _{\mathrm{vert}})`$, evaluated at a certain characteristic radius, $`r_w`$ (Silbergleit et al. 2001). The two curves plotted in Fig. 3 are $`\nu _{\mathrm{lower}}(r)=(\sqrt{2}1)\nu _{\mathrm{vert}}(r)`$ and $`\nu _{\mathrm{upper}}(r)=\nu _{\mathrm{orb}}+(1\sqrt{2})\nu _{\mathrm{vert}}(r)`$. The actual mode frequencies coincide with these at $`r=r_w`$, and at arbitrary radius differ from these by $`[\nu _{\mathrm{orb}}(r_w)\nu _{\mathrm{vert}}(r_w)][\nu _{\mathrm{orb}}(r)\nu _{\mathrm{vert}}(r)]`$, i.e., by a few percent. For the present purposes, we take the quantity $`r_w`$ to be an unknown parameter.

## 4 Conclusions

In conclusion, we note again that it is not yet known whether the accretion disk resembles a torus or, to the contrary, is geometrically thin. The 2.5 ms pulsar is the first source in which the frequencies of the twin QPOs are at once clearly related to the neutron-star spin frequency and are in a 1.4 ratio. With one observation of twin frequencies in only one source with a known rotation rate (but an unknown mass), it is impossible to determine which of the two pairs of oscillation modes discussed here is responsible for the non-linear resonance. The ambiguity will be removed if twin QPOs are observed in one or more of the four other accreting pulsars (Markwardt and Swank 2003) known, where the spin periods are different.

It is worth noting that the fact that the frequency difference between the two observed kHz QPOs is equal to one half the spin frequency of the neutron star is the first clear indication that the commonly observed millisecond modulations of the X-ray flux in low mass X-ray binaries are caused by oscillatory motions of the accretion disk around the neutron star (or black hole in other systems).

The authors thank Luciano Rezzolla for pointing out to them the importance of the radial oscillations of a torus. WK wishes to thank Prof. Shin Mineshige for generous hospitality at the Yukawa Institute, and Prof. Jacques Paul for welcoming him at the Cargese School (IESC). MAA, WK, WHL and NS acknowledge the hospitality of SISSA. This work has been supported in part by the EU Programme Improving the Human Research Potential and the Socio-Economic Knowledge Base (Research Training Network Contract HPRN-CT-2000-00137), KBN grant 2P03D01424, the Greek GSRT Grant EPANM.43/2013555, by CONACYT, and by le Departement Sciences de l’Univers du CNRS.