text:
# Microlensing of Large Sources

## 1 Introduction

With the advent of sensitive large-area synoptic surveys of the sky, it may be possible in the near future to detect non-accreting compact objects in binaries via gravitational lensing during transit of the companion stars, either main-sequence or white dwarf (Maeder, 1973; Marsh, 2001; Agol, 2002; Beskin and Tuntsov, 2002; Sahu and Gilliland, 2002). The problem of modeling the lightcurve of a transiting star has been solved for a uniform source by Agol (2002), while a limb-darkened source requires numerical integration (Agol, 2002; Sahu and Gilliland, 2002). A related problem occurs in the search for planetary companions of stars via perturbations to microlensing lightcurves which requires numerical solution for the planetary perturbation. Here we prove a theorem that the solid angle added by a lens star is equal to $`2\pi `$ times the square of the angular Einstein radius for lensing of a source which has an angular size much larger than angular Einstein radius of the lens ($`\mathrm{\S }2`$). We analytically compute the transit lightcurve of a compact object in front of uniform star, including eclipse by the lens ($`\mathrm{\S }3`$), which we use to derive the lightcurve for a limb-darkened source star ($`\mathrm{\S }4`$), which we generalize to a lens with shear and convergence ($`\mathrm{\S }5`$) as well as the case of multiple lenses in front of a large source ($`\mathrm{\S }6`$), and finally we apply these formulae to several astrophysical problems ($`\mathrm{\S }6`$). Capitalized $`R`$’s, such as the Einstein radius, $`R_E=\left(4GMc^2(D_L^1D_S^1)\right)^{1/2}`$, are in units of angular size, while lower-case $`r`$’s, such the source size $`r_{}`$, are in units of $`R_E`$, $`r_{}=R_{}/R_E`$, where $`D_L`$ and $`D_S`$ are distances to the lens and source and $`M`$ is the lens mass.

## 2 Magnification of a large source

We show first that the magnification of a large uniform source is approximately constant when the lens is in front of and away from the edges of the source. The lens equation normalized to $`R_E`$ is $`y=x1/x`$ where $`y`$ is the source position and $`x`$ is the image position. In the case of a lens in front of a uniform source of radius $`R_{}`$ for which $`R_ER_{}`$ and for which the closest distance of the lens to the edge of the source is much smaller than the Einstein radius, then $`y`$ much larger than unity. The magnification is the area contained within the outer lensed image divided the unlensed area, so $`xy+1/y`$ is large as well (the area within the inner image should in principle be subtracted off, but it is of order $`y^2`$ so that it can be ignored). Thus, the magnification is given by

$$\mu _e=\frac{1}{A}_0^{2\pi }𝑑\varphi \frac{1}{2}x^2$$
(1)
where $`x`$ is the distance to the images located at a position angle $`\varphi `$ relative to the lens center (see Figure 1) and $`A`$ is the area of the source in units of $`R_E^2`$. Squaring the lens equation gives $`y^2=x^22+1/x^2`$, but since $`x`$ is large, $`x^2y^2+2`$. Thus

$$\mu _e=\frac{1}{A}_0^{2\pi }𝑑\varphi \frac{1}{2}(y^2+2).$$
(2)
Now, $`\frac{1}{2}y^2d\varphi `$ is simply the differential area of the unlensed source, so the integral of this term gives $`A`$. The remaining term is $`d\varphi `$ which means that the additional area due to lensing is constant per differential angle.

Figure 1 shows two shaded regions between the source and image boundaries subtended by the same angle $`\delta \varphi `$. These shaded regions represent the extra solid angle added by lensing. Despite their different shape, the two regions have the same area to order $`y^2`$. This is simply due to the fact that more distant points $`y`$ have smaller displacements $`\mathrm{\Delta }y1/y`$ so that the total area, $`dA=yd\varphi \mathrm{\Delta }yd\varphi `$ is simply equal to the subtended angle, independent of the distance to the boundary or shape of the boundary. Thus, the total of the extra area is simply $`2\pi `$ (in units of $`R_E^2`$) independent of the size or shape of the source or the position of the lens (as long as it is away from the edges). The maximum area added is equal to $`2\pi `$ since for sources small compared to an Einstein radius, $`r_{}1`$, the extra area is $`2\pi r_{}2\pi `$ (however, the magnification is larger for small sources since the unlensed flux is proportional to $`r_{}^2`$). So the total magnification for a large source then becomes

$$\mu _e=1+\frac{2\pi }{A}.$$
(3)
Maeder (1973) has derived a similar result for the special case of a circular source. So we have proved that the magnification of a lens in front of a large source is constant to order $`y_{min}^2`$, where $`y_{min}`$ is the shortest distance of the center of the lens to the edge of the source in units of $`R_E`$. This creates a lightcurve which looks like a top hat.

For a circular source, $`A=\pi r_{}^2=\pi (R_{}/R_E)^2`$, so $`\mu _e=1+2(R_E/R_{})^2`$. When the lens is outside of the source, one can show that there is no extra magnification (the area is the difference of the inner and outer images of the source, so the factor of 2 cancels when taking the difference), so to a good approximation the magnification is simply

$$\mu _e=1+\left[2(R_E/R_{})^2(R_L/R_{})^2\right]\mathrm{\Theta }(r_{}z),$$
(4)
where $`z`$ is the position of the lens with respect to the origin of the circular source and $`R_L`$ is the angular size of the lens star.

## 3 Magnification of a large source near the edge

To compute the magnification of an extended source near the edge, we cannot ignore the $`x^2`$ term in the squared lensing equation, so we must include the distortion of the edge of the star. In the case of a large circular star, $`R_{}R_E`$, the edge can be approximated as a straight line and the change in magnification is simply the extra area enclosed by the image of this line. Figure 2 shows the geometry of a line being lensed by a point-mass lens. As the lens approaches the line, it is distorted outward, curving around the Einstein ring. The inner image of the line grows inside the Einstein ring, filling half of the Einstein ring when the line crosses the origin. Then, as the line crosses to the other side, the pattern repeats itself as the line recedes from the center of the disk. This creates a symmetry about the point at which the lens sits on the limb of a large source. The magnification is given by

$$\delta \mu _e=\frac{1}{r_{}^2}\frac{2(zr_{})}{\pi r_{}^2k}\left[K(k)E(k)\right],$$
(5)
where $`k^1=\sqrt{1+(zr_{})^2/4}`$. This approximates the exact magnification quite well in the large-source limit. Unfortunately, this expression is complicated enough that we cannot integrate it over limb-darkening analytically.

Figure 3 compares equation (5) with the exact equations from Agol (2002). For smaller source size, the edge has a Napoleonic bicorne shape, while larger sources look more like a tophat. The straight-edge approximation improves as the source becomes larger, and approaches a step function as $`R_{}/R_E\mathrm{}`$.

When the lens has a finite size, $`r_L=R_L/R_E>0`$, the obscuration of the images by the lens gives

$`\delta \mu `$ $`=`$ $`{\displaystyle \frac{1}{r_{}^2}}+{\displaystyle \frac{B}{\pi r_{}^2}},`$ (6)
$`B`$ $`=`$ $`s_\zeta \left({\displaystyle \frac{\pi }{2}}\varphi _0\right)(1r_L^2){\displaystyle \frac{\pi r_L^2}{2}}+G\left({\displaystyle \frac{\pi }{2}}\right)s_rG\left(\varphi _0\right)\zeta \frac{1}{2}\mathrm{tan}\varphi _0\left[|\zeta |+s_r\sqrt{\zeta ^2+4\mathrm{cos}^2\varphi _0}\right],`$ (7)
$`G(\varphi )`$ $`=`$ $`{\displaystyle \frac{\zeta }{k}}\left[K(\varphi ,k)E(\varphi ,k)\right],`$ (8)
$`\varphi _0`$ $`=`$ $`\{\begin{array}{cc}0\hfill & |\zeta |r_L|1r_L^2|\hfill \\ \mathrm{cos}^1\left(\frac{|\zeta |r_L}{|1r_L^2|}\right)\hfill & |\zeta |r_L<|1r_L^2|,\hfill \end{array}`$ (11)
where $`\zeta =r_{}z`$, $`s_\zeta =sign(\zeta )`$, $`s_r=sign(r_L1)`$.

If an observed microlensing (or occultation) event of a uniform source has a maximum depth given by $`\delta \mu _{max}`$, and if the limb-darkening function is known, then the measured depth can be fit by any $`R_L=\left[2R_E^2\delta \mu _{max}R_{}^2\right]^{1/2}`$. However, the lightcurve near the edge of the source depends differently on $`R_E`$ and $`R_L`$. In Figure 4, we show lightcurves near the edge computed from equation 6 which have the same $`\delta \mu _{max}`$, but have different ratios of $`R_E/R_{}`$ and $`R_L/R_{}`$. The cases with larger Einstein radii have a larger magnification and dip near the edge than the pure occultation case. The peak in $`\delta \mu `$ corresponds to the point at which the image outside of the Einstein radius is first occulted, at which point $`zr_{}=r_L1/r_L`$ for $`r_L>1`$. At this point, the magnification is

$$\delta \mu _{peak}=\frac{1}{2r_{}^2}\left[1\frac{r_L^2r_L^2}{\pi }\left\{K\left(\frac{2}{r_L+r_L^1}\right)E\left(\frac{2}{r_L+r_L^1}\right)\right\}\right].$$
(12)
In the limit that $`r_{}r_L1`$, this reduces to $`\delta \mu _{peak}=4^1r_{}^2r_L^2`$, which becomes small rapidly. In principle this can be used to break the degeneracy between $`R_E`$ and $`R_L`$ since the maximum magnfication has a different dependence on these quantitities.

## 4 Lensing of a large limb-darkened source

When limb-darkening is included using the parameterization $`I(R)=1\gamma _1(1\mu )\gamma _2(1\mu )^2`$ and $`\mu =[1R^2]^{1/2}`$, the magnification outside of eclipse is decreased slightly due to the darker limb, while the dip within the eclipse disappears since the dip becomes larger as the lens covers the brighter regions of the star (Figure 4). These deviations may be useful for distinguishing a brown dwarf from a planetary transit, or a close white dwarf from a terrestrial planet.

In the case of a limb-darkened star, we must integrate the magnification for a uniform source over the limb-darkened profile of the star. The magnification of a limb-darkened source is approximately

$$\mu _l=1+F_0^1_0^R_{}2\pi 𝑑R\frac{d[\delta \mu _e(R)R^2]}{dR}I(R),$$
(13)
where $`F_0=2\pi _0^R_{}𝑑RRI(R)`$ is the total flux of the source star and $`\delta \mu _e(R)=\mu _e(R)1=2(R_E/R)^2\mathrm{\Theta }(R/R_Ez)`$. Since $`\delta \mu _eR^2`$ is simply proportional to the step function, its derivative is the delta function. Integrating over this delta function picks out the surface brightness at the position of the lens. We finally find that

$$\mu _l=1+\left[2\left(\frac{R_E}{R_{}}\right)^2\left(\frac{R_L}{R_{}}\right)^2\right]\frac{I(zR_E)}{I}\mathrm{\Theta }(R_{}zR_E),$$
(14)
where $`I=F_0/(\pi R_{}^2)`$ and we have added in occultation by the foreground star. It is remarkable that although the magnification involves integrating over the surface of the entire star, it ends up depending only on the surface brightness underneath the lens! This results from the fact that the microlensing amplification of a point source becomes infinite for an exactly aligned source, which may be approximated as a delta function.

Figure 5 shows a comparison of this approximation with exact lightcurves. The agreement between the approximate lightcurves and exact lightcurves away from the edges is astoundingly good, especially considering that the analytic approximation involves a great deal less computational time. This formula gives a very good approximation to the magnification of a source with arbitrary limb-darkening for $`R_E<0.1R_{}`$.

## 5 Magnification of a large source by a Chang-Refsdal lens

We now generalize the theorem to the case of a point-mass lens embedded in a field with constant shear and convergence. A shear field distorts circles into ellipses, while a convergence field simply expands the size of a source, and since lensing by a point mass creates extra area equal to $`2\pi `$ independent of source size or shape, one might expect that the constant magnification should still apply. The Chang-Refsdal lensing equation is

$$𝐲=𝐱\left(1\frac{1}{x^2}\kappa \right)+\gamma \left(\begin{array}{cc}1& 0\\ 0& 1\end{array}\right)𝐱,$$
(15)
where $`𝐱=(x_1,x_2)`$ and $`𝐲=(y_1,y_2)`$ are the image and source coordinates, respectively, $`\kappa `$ is the convergence or normalized surface density, and $`\gamma `$ is the shear (Schneider, Ehlers & Falco, 1992). The coordinates have been chosen such that the shear axes lie along the coordinate axes. Now, the area of the image is given by the closed integral $`A_x=\frac{1}{2}𝐱\times 𝑑𝐥_𝐱`$, where $`d𝐥_𝐱`$ is the line element for the contour described by $`𝐱`$, and similarly for the source replacing $`x`$ with $`y`$. One can use the lens equation (15) to find the zeroth order solution for $`𝐱`$ which amounts to neglecting bending by the lens

$`x_{01}`$ $`=`$ $`y_1\left(1\kappa +\gamma \right)^1,`$ (16)
$`x_{02}`$ $`=`$ $`y_2\left(1\kappa \gamma \right)^1.`$ (17)
Then, the bending angle can be approximated as $`𝐱/(x_{01}^2+x_{02}^2)`$, giving a first-order solution

$`x_1=y_1\left(1{\displaystyle \frac{1}{x_0^2}}\kappa +\gamma \right)^1,`$ (19)
$`x_2=y_2\left(1{\displaystyle \frac{1}{x_0^2}}\kappa \gamma \right)^1,`$ (20)
where $`x_0^2=x_{01}^2+x_{02}^2`$. Using this solution, one can compute the area of the lensed image, $`A_x`$. For a lens at the origin, the magnification is

$$\mu _{CR}(0)=\mu _0+\frac{2\mathrm{\Theta }(|1\kappa |\gamma )}{(1\kappa )r_{}^2},$$
(21)
where $`\mu _0=|(1\kappa )^2\gamma ^2|^1`$. This is a generalization of the formula in Gould and Gaucherel (1997) which was derived for $`\kappa =0`$. For a lens located off of the origin, it can be shown numerically that the magnification is approximately constant for a lens within a large source, and then approaches $`\mu _0`$ outside of the source, so the extra magnification can be multiplied by a step function as in the point-mass case, and the specific intensity can be included as well. This gives

$$\mu _{CR}(z)=\mu _0+\frac{I(zR_E)}{I}\frac{2\mathrm{\Theta }(|1\kappa |\gamma )}{(1\kappa )r_{}^2}\mathrm{\Theta }(r_{}z).$$
(22)
Figure 6 shows a comparison of lightcurves for a Chang-Refsdal lens passing in front of a source of size $`r_{}=10`$ for different values of $`\kappa `$ and $`\gamma `$. Note that if $`|1\kappa |<\gamma `$, then the lightcurve is unperturbed. The magnification differs strongly near the edges due to the presence of caustics, but otherwise does an excellent job of describing the microlensing.

## 6 Magnification of a large source by many lenses

We next generalize the theorem to the case of lensing by a distribution of point masses. In this case the normalized lens equation is

$$𝐲=𝐱\underset{i}{}m_i\frac{𝐱𝐱_𝐢}{|𝐱𝐱_𝐢|^2},$$
(23)
where $`m_i`$ is the fraction of the total mass in each lens ($`_im_i=1`$) and $`𝐱_𝐢`$ is the position of each lens. As in the single-lens case, we can square this equation, disregard terms of order $`x^2`$, and rearrange terms to give

$$y^2=x^222\underset{i}{}m_i\frac{𝐱_𝐢(𝐱𝐱_𝐢)}{|𝐱𝐱_𝐢|^2}.$$
(24)
This looks similar to the single-lens case except for the extra term which appears since each lens is not at the center of the coordinate system. One can approximate $`𝐱𝐲`$ and then carry out the integral over the last term. In the case of a circular source, this term gives zero. So, we end up deriving the same result as in the single-lens case, equation (4), the only difference is that the Einstein radius is determined by the sum of all masses in front of the lens; the condition $`R_ER_{}`$ must still hold.

Now consider a random star field with a surface number density of stars $`n_{}`$. The lensing optical depth is defined as $`\tau _{}=\pi R_E^2n_{}`$, where $`R_E`$ is the Einstein radius for each star. The mean number of stars in front of a large source is then $`N_{}=n_{}\pi R_{}^2`$, so if $`R_E=\sqrt{N_{}}R_ER_{}`$, or $`\tau _{}1`$, then we can apply the above result to find the mean magnification: $`\mu =1+2\tau _{}`$. Of course, a random distribution of stars will create Poisson fluctuations in the number of stars in front of the source. Thus, the standard deviation of the magnification for an ensemble of star fields is $`\sigma _\mu =2\sqrt{\tau _{}}R_E/R_{}`$. This result has already been derived by Refsdal and Stabell (1991) using a different computational technique.

## 7 Applications

The formulae derived here will be useful for discovering compact objects in eclipsing binaries. A measurement of lensing during a transit by a compact object of a companion star places a constraint on the lens mass (Sahu and Gilliland, 2002; Beskin and Tuntsov, 2002). Brown dwarfs (BD), white dwarfs (WD), black holes (BH), and neutron stars (NS) in binaries of semi-major axis $`a`$ each have Einstein radii, $`R_E=\sqrt{4GMc^2a}`$, much smaller than the size of a companion main-sequence star (MS) for $`a1`$ AU and $`M1M_{}`$, and thus the approximation of equation (14) will apply. Figure 8 shows an example lightcurve of a white dwarf orbiting a main sequence star for $`M_{WD}=0.7M_{}`$, $`M_{MS}=1M_{}`$, $`a=0.1`$ AU, $`t_{WD}=10`$ Gyr, $`\gamma _1=0.6,\gamma _2=0`$, and the observer in the equatorial plane of the binary. In the upper left panel, the analytic approximation (equation 14) is plotted as a dotted line, which is barely distinguishable from the exact lightcurve (solid line). The transit/eclipse duration is about 3 hours. There is slight tidal distortion of the main-sequence star with an amplitude of $`5\times 10^5`$ and sinusoidal dependence of twice the frequency of the binary. In this particular example, $`r_L=1`$ and $`r_S=88`$, or the Einstein radius is $`R_E=0.011R_{}`$. The V-band magnitude of the WD is $`M_{V,WD}=17`$, while the MS star is $`M_{V,MS}=4.8`$, which makes the eclipse of the white dwarf by the main-sequence star (at 160 hours) rather faint. The maximum deviation of the lightcurve is to be $`8\times 10^5`$, this is modified somewhat due to limb-darkening. To detect such a transit requires better than $`0.01`$ millimagnitude sensitivity over several hours, which is a requirement that should be met by the planned Kepler and Eddington missions.

A simple estimate of the size of the lensing signal can be made in the limit that the size of the lensing object can be neglected (usually the case for neutron stars, black holes, and white dwarfs with large semi-major axes). The extra flux caused by lensing is simply $`2r_{}^2F`$ (neglecting limb-darkening), where $`F`$ is the flux of the unlensed companion star. The duration of transits is of order $`2R_E/v`$, where $`v`$ is the relative velocity of the binary, while the number of transits scales as $`T_{obs}/P=T_{obs}v/(2\pi a)`$, where $`T_{obs}`$ is the total observing time (assumed to be much larger than the period of the binary and the monitoring is assumed to be continuous), $`a`$ is the semi-major axis of the binary (assumed to be circular). Then the total fluence added by the lensing is equal to the lensing flux times the transit duration times the number of transits, or

$$\delta E=2F\frac{R_E^2}{R_{}^2}\frac{2R_E}{v}\frac{T_{obs}v}{2\pi a}=\frac{8}{\pi }\frac{R_g}{R_{}}FT_{obs},$$
(25)
where $`R_g=GM_{CO}/c^2`$ is the gravitational radius of the lensing compact object. Remarkably, any dependence on the semi-major axis or mass of the main-sequence star has disappeared! This is due to the fact that more distant binaries have a larger flux signal and longer transit duration, but fewer transits during a given experiment duration. Since $`E=FT_{obs}`$ is the fluence of the unlensed star, the total signal is

$$\frac{\delta E}{E}=\frac{8}{\pi }\frac{R_g}{R_{}}=5\times 10^6\frac{m}{R_{}},$$
(26)
where $`m`$ is the mass of the lensing, compact object in solar masses and $`R_{}`$ is the radius of the main-sequence star in units of solar radii. This relation suggests a technique for measuring the mass of the eclipsing object: first, estimate the radius of the main-sequence companion via spectral or photometric modeling (the distance can be unknown), and then measure $`\delta E/E`$ to derive $`m=2\times 10^5R_,\delta E/E`$.

The foregoing formula is only true for an equatorial transit of a main-sequence star in a band which is not limb-darkened. Including limb-darkening and a non-equatorial transit, one finds

$$\frac{\delta E}{E}=\frac{8}{\pi }\frac{R_g}{R_{}}\frac{2xI_{line}}{I_{disk}},$$
(27)
where $`x=\sqrt{1b^2}`$, $`b`$ is the minimum impact parameter of the lens with respect to the center of the source star in units of $`R_{}`$, $`I_{line}`$ is the mean surface brightness of the source along the path of the lens, and $`I_{disk}`$ is the mean surface brightness of the source. For quadratic limb-darkening, we find

$$H(b)=\frac{2xI_{line}}{I_{disk}}=\frac{(1\gamma _1\gamma _2)x+\frac{\pi }{4}(\gamma _1+2\gamma _2)x^2\frac{2}{3}\gamma _2x^3}{\pi \left(1\frac{1}{3}\gamma _1\frac{1}{6}\gamma _2\right)}.$$
(28)
If $`R_{}`$ and $`\gamma _1,\gamma _2`$ are known from, say, modeling of the spectrum of the main-sequence star, then equation 27 has two unknowns, $`b`$ (or $`x`$) and $`M_{CO}`$. Another equation can be obtained from the ratio of the duration of the transit, $`t_{transit}`$, to the period of the orbit, $`P`$. Since $`P`$ depends on the semi-major axis, which depends on the masses of both stars, one finds

$$\sqrt{1b^2}R_{}=t_{transit}\left[2\pi GP^1\left(M_{MS}+M_{CO}\right)\right]^{1/3},$$
(29)
which constitutes another equation for $`b`$ and $`M_{CO}`$ assuming that $`M_{MS}`$ is known. The mass of the compact object may be found by eliminating $`M_{CO}`$ in equations 27 and 29, which yields a 6th order equation for $`x`$ which can be solved numerically for $`x`$ which may then be used in equation 27 to find $`M_{CO}`$. If the data are of sufficient quality then $`b`$ may be found by fitting the shape of the transit and then $`M_{CO}`$ may be found directly from equation 27 (this is described in Mandel & Agol (2002) for the occultation case). Uncertainties in $`M_{},R_{},\gamma _1,`$ and $`\gamma _2`$ will contribute to uncertainty in $`M_{CO}`$ and $`b`$. Another technique for finding $`b`$ relies on measuring the difference of the transit depth in different bands. Figure 7 shows the ratio of $`H(b)`$ for the V and I bands as a function of impact parameter for a solar-type star. The variation is nearly a factor of 3 for $`b`$ ranging between 0 and 1. Given this ratio, $`x`$ may be solved for using the quadratic equation in $`x`$ determined by equating the functions 28 for V and I with the observed value. The value of $`b`$ measured from this relation can then be used to estimate $`H(b)`$ and determine $`R_g/R_{}`$.

The transit detection probability is $`R_{}/a`$, from which we can estimate the fraction of binaries for which transits can be detected assuming the the required fluence sensitivity can be achieved. This is considered in detail for WD-MS binaries by Farmer and Agol (2003) who find that $`>10^2`$ such systems might be found with the Kepler and Eddington missions, and for WD-WD, WD-NS, WD-BH binaries by Beskin and Tuntsov (2002). Here we make some general considerations for the detection of NS-MS and BH-MS binaries. For the binary to not accrete, the main-sequence star must have a size smaller than the Roche-lobe size, $`R_L`$, by some factor $`f_R=R_{}/R_L`$. To a good approximation, the Roche-lobe is given by $`R_L=0.5af_M^{1/3}`$ where $`f_M=M_{MS}(M_{MS}+M_{CO})^1`$ and the $`CO`$ subscript stands for compact object (Pringle & Wade, 1985). Typically, binaries show a distribution of semi-major axes which is constant in $`\mathrm{ln}a`$, which we assume to be the case from $`a_{min}`$ to $`a_{max}`$ for BH-MS and NS-MS binaries (clearly, reality is likely to be more complicated due to the effects of kick-velocities during the formation of the compact object, e.g. in a supernova collapse), common-envelope evolution of the massive star, and tidal evolution due to distortion of the main-sequence star). Since $`a_{min}`$ is set by the Roche-lobe overflow limit, $`a_{min}=2R_{}f_R^1f_M^{1/3}`$, the probability for lensing scales as

$$p_{lens}=R_{}/a_{min}=0.5f_Rf_M^{1/3}[ln(a_{max}/a_{min})]^1.$$
(30)
Note that since $`a_{min}R_{}`$, the probability is independent of the size of the star! For a population of solar-mass stars orbiting 10 $`M_{}`$ black holes with, say, $`a_{max}=10^5a_{min}`$ and $`f_R=0.5`$, the probability of lensing is about 1%. The detected stars will be heavily concentrated toward $`a_{min}=9R_{}`$ for this choice of parameters, so $`\delta \mu _e10^3`$.

The total number of non-accreting BH-MS or NS-MS binaries in our galaxy is quite uncertain. Tutukov & Yungelson (2002) estimate that $`(36)\times 10^4`$ BH-MS binaries and $`(13)\times 10^5`$ NS-MS binaries reside in the Milky way, while Lipunov, Postnov and Prokhorov (1996) estimate $`10^{2.8}`$ BH-MS and $`10^4`$ NS-MS binaries. Given the above assumptions for $`f_R`$ and $`f_M`$, there exists a total of $`10^{23.5}`$ NS-MS transiting binaries and $`10^{0.62.6}`$ BH-MS transiting binaries for our assumed parameters, a formidably small population to find, especially considering that the signal is at a level of $`10^310^4`$ which may be swamped by stellar variability. Space-based transit searches can easily reach this precision over the transit timescale; however, they will monitor only a small fraction of the number of stars required to find BH-MS or NS-MS binaries. Globular clusters may be the best locations to search for these binaries due to the concentration of stars and the large over-representation of accreting BH-MS and NS-MS binaries relative to the rest of the Galaxy; however, the precision and time coverage of current surveys is not sufficient for a detection.

Several microlensing searches are taking place with the goal of finding deviations in the lightcurve induced by planets (Albrow, 1995; Bennett and Rhie, 1996). In the case of a terrestrial planet companion of a star lensing a source at the bulge, the Einstein radius is $`R_E=2.7\times 10^{11}\mathrm{cm}(M/M_{})^{1/2}(D/4\mathrm{k}\mathrm{p}\mathrm{c})^{1/2}`$. In the case of a Jupiter-mass planet in a globular cluster lensing a source in the same cluster, the Einstein radius is $`R_E=4.8\times 10^{10}\mathrm{cm}(M/M_J)^{1/2}(D/1\mathrm{p}\mathrm{c})^{1/2}`$ (however, the optical depth for such events is negligible). These are somewhat smaller than the size of a red-giant star, which comprise a large fraction of microlensing targets given their bright luminosities. The region of the microlensing lightcurve near the planet can be treated as a Chang-Refsdal lens which means that lensing of a red-giant star by a star-planet system will create fluctuations of the magnitude given by equation (22). This generalizes a result from Gould and Gaucherel (1997), who computed the amplification of a large source lensed by a point mass with shear, to include limb-darkening and convergence. Gould (1996) has argued for the monitoring of giant stars to look for planets given their brightness, lack of confusion, and large cross section. If such events are found, the limb-darkened lightcurves can be approximated by equation (22).

## 8 Conclusions

We have shown that lensing of a large source adds an area which is approximately constant when the source dimensions are much larger than the Einstein radius, which leads to an accurate approximation for the lensing lightcurve in the small amplification limit, typically $`\delta \mu <0.01`$ for this approximation to be relevant. We have extended this result to include limb-darkening, shear, and convergence. Since many large-area accurate photometric searches are being carried out or planned, it is likely that this formalism will be relevant for planning observations and discovering dark objects via lensing. The most relevant applications are likely to be lensing by white dwarfs in binaries and lensing perturbations by planets during Galactic microlensing events. Using the formalism of this paper we have rederived a result that the fluctuations of a distant source lensed by a random star field at low optical depth varies as the square root of the optical depth times the ratio of the Einstein radius to the source size, first derived by Refsdal and Stabell (1991).

Support for E.A. was provided by the National Aeronautics and Space Administration through Chandra Postdoctoral Fellowship Award PF0-10013 issued by the Chandra X-ray Observatory Center, which is operated by the Smithsonian Astrophysical Observatory for and on behalf of the National Aeronautics Space Administration under contract NAS 8-39073.


format: text

text:
# Comment on “Experimental Nonlocality Proof of Quantum Teleportation and Entanglement Swapping”

In a recent paper in which they demonstrate quantum nonlocality for photons that never interacted , Jennewein, Weihs, Pan, and Zeilinger mention a seemingly paradoxical situation that arises when Alice’s Bell-state analysis is delayed long after Bob’s measurements. According to the authors, “This seems paradoxical because Alice’s measurement projects photons 0 and 3 into an entangled state after they have been measured,” in other words, “This means that the physical interpretation of his (Bob’s) results depends on Alice’s later decision.” Still quoting from their paper: “Therefore, this result indicates that the time ordering of the detection events has no influence on the results and strengthens the argument of A. Peres : this paradox does not arise if the correctness of quantum mechanics is firmly believed.”

The authors are correct about the time order of the detection events. Actually, Einstein, Podolsky, and Rosen (EPR) correlations, teleportation, and entanglement swapping are different aspects of the very same phenomenon when viewed from different Lorentz frames that move relative to each other . On the other hand, they are incorrect when they conclude that the physical interpretation of Bob’s results depends on Alice’s later decision, and Peres’ statement, to which they subscribe, does little to explain why there is no paradox. The seeming paradox arises because, apparently, a measurement in the future projects the photons in an entangled state in the past  . If this were true, this would imply that not only local realism, as shown by Bell , but also nonlocal realism cannot mimic quantum mechanics, since there can manifestly be no causal connection between the events observed by Bob (even accepting the possibility of superluminal interaction). In particular, Bohm’s interpretation  would be condemned. Therefore, it is important to show why there is no paradox.

Naturally, firm belief is not the most appropriate way to try to solve a paradox. Actually, the reason why such a paradox does not occur is quite simple and has already been explained . The point is that whenever Bob performs his measurement first, photons 1 and 2 are projected into well-defined polarization states. As a consequence, Alice will perform an experiment similar to the one performed by Hong, Ou, and Mandel , but in which the photons that arrive together at the beam splitter can have different polarizations. Since the probability of coincident detection depends on the relative polarization of the photons, by selecting the events in which the detectors click in coincidence, Alice obtains a subset that behaves as if it consisted of entangled pairs of distant particles. Naturally, there is nothing paradoxical in this result. The point is that for Alice the probability of coincident detection will depend on Bob’s experimental outcome. Therefore, no influence of future actions on past events needs to be invoked. Very probably, the authors were not aware of this simple explanation, or else they would realize that there is no reason to perform the experiment discussed by Peres, since it adds nothing to what is already known. The photon-photon entanglement in type-II spontaneous parametric down-conversion and the behavior of photons that arrive together at a beam splitter have already been experimentally verified. Naturally, the authors have performed an important and interesting experiment, but their extension to the situation in which Alice’s measurement occurs after Bob’s is unnecessary and their interpretation misleading.

ACKNOWLEDGMENT

Financial support was provided by Brazilian agency CNPq-Institutos do Milênio- Informação Quântica


text:
# CAN ONE MAKE A POWDER FORGET ITS HISTORY?

## Abstract

It is shown that computer simulations can qualitatively reproduce experiments, where a powder of cohesive, round, hard particles is periodically deformed at constant volume. Two types of initial configurations are considered: Uniaxially precompacted ballistic deposits and biaxially precompacted DLA-clusters. Both initial configurations had the same volume fraction, but due to the different precompaction procedure completely different principal stresses. After a transient which lasts only less than a period, the stresses follow the same periodic function, i.e. the powder forgot its history.

1. INTRODUCTION

One of the fundamental open questions in the physics of powders is how to characterize their state uniquely, the particle configuration in general depending on the deformation history of the powder. This is in marked contrast to thermal equilibrium systems, which statistically explore all available configurations, so that temporal correlations normally die out quickly. A powder, however, is “jammed” or frozen in, if it is not driven by external forces. Even if it flows, it does not explore its alternative configurations on a time scale which is short compared to the inverse strain rate.

This has macroscopic consequences. For example the stress in a compacted cohesive powder not only depends on the solid fraction but also, by what kind of compaction the solid fraction was reached. It is believed, that at least the fabric tensor is needed in order to determine the stress state. Any microscopic theory of the dynamical evolution of the state and hence of powder flow requires an answer to the question, how to characterize such a state. Therefore it is very important also for macroscopic predictions needed e.g. in process engineering.

In this paper we address only a partial aspect of this question. It is known, that the stress components in a powder, which is quasi-statically deformed at fixed volume, approach values independent of the precompaction history . This stress state depends on the volume fraction, but also on the total deformation at constant volume: As one continues to shear, the stress components increase. The point we want to make in this paper is that there are also other ways to make the powder forget its precompaction history. As an example we investigate the periodic kneading of the powder at constant volume. In contrast to steady state flow the principal axes switch back and forth periodically in this case. This has the advantage that the sample returns to its initial shape periodically. In every period steady state flow properties may be reached or not, depending on how long the period lasts.

Periodic kneading was studied in a different context previously . There it was shown for limestone that the stresses are periodic functions with decreasing amplitude. After a few periods the amplitude becomes constant. In the present paper we compare experiments on carbonyl-iron powder and computer simulations. In contrast to limestone the particles are spherical and very hard. As we are going to show, the amplitude of the periodic stresses becomes constant already within one cycle in this case. Simultaneously all stress tensor differences due to different sample preparations vanish within this short transient, i.e., the powder forgets its history.

2. METHODS

2.1 Experimental Setup
The true biaxial shear tester allows the deformation of a brick shaped bulk solid specimen. The deformation mechanism is composed of a bottom and a top plate which are fixed at a distance of 36.5 mm from each other and four side plates which can be moved independently from a maximum distance of 130 mm to a minimum distance of 70 mm, always staying perpendicular to their adjacent ones. Thus the bulk solid specimen remains brick shaped and principal strains equal the normal strains (Fig. 1a):

$$ϵ_{ij}=0,\text{ where }ij;ϵ_{ii}=ϵ_i.$$
(1)
The side plate’s velocities can be set in steps of 0.03 $`\mu `$m/s up to a maximum velocity of 65 $`\mu `$m/s corresponding to maximum shear rates of 0.1% per second. On the specimen’s border, the complete stress state is measured with three five component load cells , which are installed in the bottom and in two perpendicular load plates, LC 1, 2, and 3 in Fig. 1b.

For the tests with the true biaxial shear tester, silicon grease is spread on bottom, lateral and top plates. Then they are covered with highly flexible rubber membranes in such a way that the membranes are evenly deformed with the borders of the powder specimen. The friction between specimen and load plates thus is minimized, hence shear stresses are negligible. Therefore, the normal stresses on the specimen borders are principal stresses. The principal stress axes coincide with the deformation axes (Fig. 1a):

$$\sigma _{ij}=0,\text{ where }ij;\sigma _{ii}=\sigma _i.$$
(2)
In this paper we use the engineering sign convention for the stress, i.e. a compressive stress is defined as positive. Deformations and stresses are continuously measured using a process computer, which also controls the deformation process . As experiments with finite strains are carried out, we define the strains using the natural strain increment  and obtain

$$ϵ_i=\mathrm{ln}\frac{L_i}{L_{i0}}$$
(3)
for the strains in $`x`$-, $`y`$\- and $`z`$-direction and for the volume strain respectively:

$$ϵ_V=\mathrm{ln}\frac{V}{V_0}=\mathrm{ln}\frac{L_xL_yL_z}{L_{x0}L_{y0}L_{z0}}=ϵ_x+ϵ_y+ϵ_z$$
(4)
In order to perform tests at constant volumes, one can easily control the condition $`ϵ_x+ϵ_y=`$const. ($`ϵ_z=0`$ during the complete measurement due to the fixed bottom and top plates). For test preparation, the powder is first sieved in $`z`$-direction into the tester. Then, it is precompacted in the same direction in order to homogenize the sample. After this, the sample is sheared with $`\dot{ϵ}_x=\dot{ϵ}_y`$, i.e., the sample’s volume remaining constant. The directions of shear are changed periodically. All experimental results presented in this paper have been obtained with carbonyl-iron powder with a median particle size of about 2 microns. This powder was chosen because its primary particles are perfectly round, in analogy to the particles used in the computer simulations.

2.2 Simulation Method
For the computer simulations we use a discrete element method, the Contact dynamics method . In contrast to soft particle molecular dynamics, the volume exclusion of perfectly rigid particles and Coulomb’s friction law are implemented exactly . The cohesion between the particles is included in a simplified model as an attractive force acting within a range leading to a force and an energy barrier to break a contact between two particles . Furthermore, rolling friction between the round particles is implemented, also leading to the stabilization of large pores . We use two dimensional simulation because the influence of the torsional degree of freedom for a contact present in three dimensions is not well understood up to now.

All simulation results are given in the following units: As unit of length we take the particle radius $`r`$, as unit of mass $`\rho r^3`$, where $`\rho `$ is the mass density, and as unit of time the inverse strain rate $`|\dot{ϵ}_x|^1=|\dot{ϵ}_y|^1`$ (which was constant during the periodic kneading). In these units the simulation parameters were: The cohesion force between two particles, $`F_\mathrm{c}=10^{22}\rho r^4|\dot{ϵ}_x|^2`$, the range of the cohesion force, $`d=10^4r`$, the friction coefficient, $`\mu =0.3`$, and the rolling friction coefficient, $`\mu _\mathrm{R}=0.2r`$ (which is the torque resisting the rolling motion divided by the normal force $`F_\mathrm{n}+F_\mathrm{c}`$ at the contact). We like to point out that for the typical values of the carbonyl-iron powder experiment \[$`r1`$ $`\mu `$m, $`\rho 10^4`$ kg/m3 and $`|\dot{ϵ}_x|10^4`$ s-1\] our simulation parameter $`F_\mathrm{c}10^6`$ N is a reasonable value for the cohesion force between the particles.

The walls were cohesive, but frictionless, so that no tangential stresses could occur, as in the experiment. This guarantees, that the principal axes of the stress and the strain tensors were collinear.

3. RESULTS

For the two-dimensional simulations, we prepared two different kinds of initial configurations in the following way (cf. Fig. 2): One is a ballistic deposit , which was precompacted in the direction of the deposition  up to a volume fraction of $`\nu =0.42`$. The other is an off-lattice DLA cluster  centered in a square box, and then precompacted biaxially to the same volume fraction, shrinking the square box symmetrically. Note that the first preparation procedure tends to lead to anisotropic, but homogeneous configurations, while the second one results in inhomogeneous, but isotropic configurations (the initial DLA cluster is essentially round and has a higher density at the center). This was done deliberately to demonstrate that even in this case both configurations reach the same (periodic) stress/strain relationship.

The precompacted systems are now used as initial configurations for the periodic shearing at constant volume. As in the experiment, the shearing was performed with constant logarithmic strain rate. The maximal strain was $`36\%`$ (with respect to the square shape). Figure 3 shows the time dependence of the two principal stresses, averaged over $`20`$ independent runs (with different, but similarly prepared initial configurations), including the precompaction. Moreover the data were smoothed by averaging them over a running time window of 5% of the period.

Clearly, the stresses are qualitatively different during the precompaction: For configuration type 1 (ballistic deposit) and uniaxial precompaction one gets anisotropic stresses, the larger stress belonging to the direction of compaction. For configuration type 2 (DLA-cluster) and biaxial precompaction both principal stresses increase in the same way, but when the square is deformed into a rectangular shape at constant volume (first 1/4 cycle), the principal stress in the expanding direction drops to a small positive value, while the other component stays large. This behavior agrees qualitatively with the experiment with the biaxial precompaction, Fig. 4.

The time dependent stresses within the first half cycle ($`t/t_{\mathrm{cycle}}=0`$ to $`1/2`$) in Fig. 3a are qualitatively different from the ones in all subsequent half cycles. Similarly the first 1/4 cycle ($`t/t_{\mathrm{cycle}}=1/4`$ to $`0`$) in Fig. 3b is different. However, after these transients, it seems that the same stress pattern is repeated period after period, also in the experiment Fig. 4. In order to compare the stress patterns for the two cases, Fig. 3a and b, we averaged them over 9 repetitions.

The results are shown in Fig.5. Within the error bars the stress pattern for the samples obtained from uniaxially precompacted ballistic deposits is the same as for those obtained from biaxially precompacted DLA-clusters. Moreover $`\sigma _y`$ simply coincides with $`\sigma _x`$ shifted by half a period, as it should. This shows, that after the transients the stresses do not show any dependence on the preparation history. This is also reflected in Fig.2e and f, which show that after one cycle the particle configurations are qualitatively indistinguishable.

Figure 5b shows the stress pattern for a higher solid fraction. Here the maximum strain compared to the square shape was 16%. Again one gets the same pattern in spite of the different types of original configurations and precompaction histories. However, the shape and amplitude of this pattern is completely different. This shows that the periodic stress pattern depends on the solid fraction.

4. CONCLUSIONS

We investigated the strain-dependent stress in a cohesive powder, which was periodically deformed at constant volume and strain rate. It was shown that contact dynamics simulations agree qualitatively with experiments on carbonyl-iron powder.

For the simulations two types of initial configurations were prepared. One had anisotropic stress and fabric tensors (uniaxially compacted ballistic deposits), the other was isotropic (biaxially compacted DLA-clusters), but both types of configurations had the same volume fraction. These initial configurations were then periodically deformed at constant volume. After a transient which lasts only less than a period, one obtains for both types of initial configurations the same periodic stress pattern. This leads us to the conclusion that periodic shear like steady state flow makes the powder forget its history quickly. It was also shown, that the periodic stress pattern depends on the solid fraction.

5. REFERENCES

6. NOMENCLATURE

$`ϵ`$ – strain tensor, dimensionless
$`\sigma `$ – stress tensor, Pa
$`L`$ – size of the sample chamber, m
$`V`$ – volume of the sample chamber, m3
$`r`$ – particle radius, m
$`\rho `$ – mass density, kg m-3
$`d`$ – range of the cohesion force, m
$`F_c`$ – cohesion force, N
$`F_n`$ – normal force, N
$`\mu `$ – coefficient of tangential friction, dimensionless
$`\mu _\mathrm{R}`$ – coefficient of rolling friction, m
$`\nu `$ – volume fraction, dimensionless

7. ACKNOWLEDGMENT

We thank Guido Bartels and Henning Knudsen for fruitful discussions. This research was supported by DFG within grants WO 577/3-1, SCHW 233/29-1 and SFB 445, as well as by Federal Mogul GmbH and BASF AG Ludwigshafen.


text:
# Cross-Correlation and Maximum Likelihood Analysis: A New Approach to Combine Cross-Correlation Functions

## 1 Introduction

Ever since the seminal works of Simkin (1974) and Tonry & Davis (1979), the cross-correlation technique to measure astronomical Doppler shifts has become extremely popular. The advent of digitized spectra and computers made it the preferred method. It has been applied in all the astronomical fields that require the measurement of radial velocities from observed spectra, ranging from binary and multiple stellar systems to cosmology. In recent years, improvements in the precision of radial velocities measured through cross-correlation led to the detection of many extrasolar planets (e.g., Mayor & Queloz, 1995).

The cross-correlation technique is conceptually simple and can be presented in an intuitive manner. Nevertheless, its properties have been studied extensively, in an effort to improve its precision and overcome its few limitations. Thus, various methods have been suggested to estimate its precision (e.g., Tonry & Davis, 1979; Connes, 1985; Murdoch & Hearnshaw, 1991). TODCOR \- a Two-Dimensional Correlation technique, was introduced as a generalization of cross-correlation meant to measure the Doppler shifts of two blended spectra (Zucker & Mazeh, 1994).

Due to the progress in detector technology, the modern spectrographs pose a new challenge to the technique, by producing multi-order spectra. In order to maximize the precision of the measured velocities, we need to combine the spectral information in all the relevant orders, potentially reducing the signal-to-noise ($`S/N`$) required to achieve a specified precision. One approach was suggested by Connes (1985) and Bouchy et al. (2001), who calculated the cross-correlation for each order separately, and then calculated a weighted average of the resulting cross-correlation functions. The weights used by Connes and Bouchy et al. reflect their estimate of the $`S/N`$ in the corresponding orders.

This paper offers another approach to study the properties of cross-correlation. It is shown that under certain assumptions, measuring the radial velocity through cross-correlation is equivalent to using maximum-likelihood analysis for the measurement. This approach has been already used in communication theory in the context of signal detection (e.g., Proakis, 1995). This approach leads in a natural way to an error estimate and to a new scheme for combining cross-correlation functions. This new scheme is shown, theoretically and through simulations, to be superior over the other schemes proposed, both in terms of precision and in terms of signal detection capability.

Section 2 shows how cross-correlation can be derived from maximum-likelihood theory. It derives an error-estimate and shows simulations to test this error estimate. The new combination scheme is derived from maximum-likelihood principles in Section 3. Several simulations to demonstrate the power of the technique are also presented in this Section. The paper is concluded by a few remarks in Section 4.

## 2 From maximum likelihood to cross-correlation

### 2.1 The basic assumptions

Let $`f(n)`$ denote the observed spectrum, whose Doppler shift is to be found by correlating it against $`g(n)`$ – the ‘template’ of zero shift. Both the stellar spectrum and the template are assumed to be described as functions of the bin number – $`n`$, where $`n=A\mathrm{ln}\lambda +B.`$ Thus, the Doppler shift results in a uniform linear shift of the spectrum (Tonry & Davis, 1979).

Now let us introduce the statistical model. Under this model we assume that the observed spectrum was produced by multiplying the template by a scaling constant ($`a_0`$), shifting it ($`s_0`$ bins) and adding a random white gaussian noise with a fixed standard deviation ($`\sigma _0`$), i.e.:

$`f(n)`$ $`=`$ $`a_0g(ns_0)+d_n,`$
$`d_n`$ $``$ $`N(0,\sigma _0^2)`$
Obviously, $`a_0`$ and $`s_0`$ are not known in advance. In the literature, $`\sigma _0`$ is usually assumed to be known, but in the following derivation we will assume no prior knowledge of $`\sigma _0`$.

A few simplifying assumptions are usually applied, which are also used in this work. Thus, the spectra involved are assumed to be “continuum-subtracted”, i.e. a best-fit low-order polynomial has been subtracted from the spectra. As a result, the spectra will have a zero mean:

$`{\displaystyle \underset{n}{}}f(n)`$ $`=`$ $`0`$
$`{\displaystyle \underset{n}{}}g(n)`$ $`=`$ $`0`$
The number of bins is typically very large (usually in the order of thousands) and thus, for small shifts (even tens of bins), one can neglect edge effects and assume that for each shift, the overlap length is $`N`$ – the total number of bins.

### 2.2 Maximum-likelihood estimation

The likelihood is defined by the probability of the observed results, under the assumed model, as a function of the model parameters. In our case:

$`L`$ $`=`$ $`{\displaystyle \underset{n}{}}\left({\displaystyle \frac{1}{\sqrt{2\pi \sigma ^2}}}\right)\mathrm{exp}\left\{{\displaystyle \frac{[f(n)ag(ns)]^2}{2\sigma ^2}}\right\}`$
$`=`$ $`\left({\displaystyle \frac{1}{\sqrt{2\pi \sigma ^2}}}\right)^N\mathrm{exp}\left\{{\displaystyle \underset{n}{}}{\displaystyle \frac{[f(n)ag(ns)]^2}{2\sigma ^2}}\right\}`$
The usual practice in maximum-likelihood estimation theory is to find the values of the parameters which maximize the natural logarithm of the likelihood fuction. In our case this logarithm is:

$$\mathrm{log}L=N\mathrm{log}\sigma \frac{1}{2\sigma ^2}\underset{n}{}[f(n)ag(ns)]^2+\mathrm{𝑐𝑜𝑛𝑠𝑡}.,$$

where $`\mathrm{𝑐𝑜𝑛𝑠𝑡}.`$ is independent of the parameters. As it turns out, the values of $`\sigma `$, $`a`$ and $`s`$ which maximize this function are the solutions of the equations:

$`\widehat{\sigma }^2`$ $`=`$ $`{\displaystyle \frac{1}{N}}{\displaystyle \underset{n}{}}[f(n)\widehat{a}g(n\widehat{s})]^2`$
$`\widehat{a}`$ $`=`$ $`{\displaystyle \frac{_nf(n)g(n\widehat{s})}{_ng^2(n)}}`$
and $`\widehat{s}`$ is the value which maximizes the cross-covariance function, which is defined as:

$$R(s)=\frac{1}{N}\underset{n}{}f(n)g(ns).$$

Let us introduce some notation which will facilitate the derivation from here on. Thus, let $`s_f`$ and $`s_g`$ denote the standard deviation of the analysed spectrum and the template:

$`s_f^2`$ $`=`$ $`{\displaystyle \frac{1}{N}}{\displaystyle \underset{n}{}}f^2(n)`$
$`s_g^2`$ $`=`$ $`{\displaystyle \frac{1}{N}}{\displaystyle \underset{n}{}}g^2(n).`$
Define $`C(s)`$ as the cross-correlation function of $`f`$ and $`g`$:

$$C(s)=\frac{R(s)}{s_fs_g},$$

where the normalization forces $`|C(s)|`$ to be smaller then unity. Since $`C(s)`$ is proportional to $`R(s)`$, $`\widehat{s}`$ can be defined as the shift which maximizes the cross-correlation function.

Using these notations, the expressions for $`\widehat{\sigma }`$ and $`\widehat{a}`$ reduce to:

$`\widehat{\sigma }^2`$ $`=`$ $`s_f^2[1C^2(\widehat{s})]`$
$`\widehat{a}`$ $`=`$ $`{\displaystyle \frac{R(\widehat{s})}{s_g^2}}.`$
Substituting these expressions again into the expression for $`\mathrm{log}L`$, we get:

$$\mathrm{log}L=\frac{N}{2}\mathrm{log}[1C^2(\widehat{s})]+\mathrm{𝑐𝑜𝑛𝑠𝑡}.$$
(1)
Thus, we see that the likelihood is an increasing monotonic function of the squared cross-correlation. The dependence on the square of the cross-correlation, instead of the cross-correlation itself, means that the likelihood rises for negative values of the correlation. This behaviour is related to the formal possibility that the scaling factor – $`a`$ – be negative. This has no meaning in the astrophysical context, since it means all the expected asborption lines in one component appear as emission lines.

### 2.3 Error estimate

Maximum-likelihood theory also includes a way to estimate the errors, or the confidence intervals, for the parameters. The covariance matrix of the parameters has been shown, at non-pathological cases, to be the negative inverse of the Hessian matrix of the log-likelihood function at its maximum (e.g., Kendall & Stuart, 1967). Let us calculate this matrix and substitute the values at the maximum, treating the shift, $`s`$, as a continuous variable:

$`{\displaystyle \frac{^2\mathrm{log}L}{a^2}}|_{maximum}`$ $`=`$ $`N{\displaystyle \frac{s_g^2}{\widehat{\sigma }^2}}`$
$`{\displaystyle \frac{^2\mathrm{log}L}{\sigma ^2}}|_{maximum}`$ $`=`$ $`{\displaystyle \frac{N}{\widehat{\sigma }^2}}{\displaystyle \frac{3}{\widehat{\sigma }^4}}{\displaystyle \underset{n}{}}[f(n)\widehat{a}g(n\widehat{s})]^2={\displaystyle \frac{2N}{\widehat{\sigma }^2}}`$
$`{\displaystyle \frac{^2\mathrm{log}L}{s^2}}|_{maximum}`$ $`=`$ $`{\displaystyle \frac{\widehat{a}NR^{\prime \prime }(\widehat{s})}{\widehat{\sigma }^2}}`$
$`{\displaystyle \frac{^2\mathrm{log}L}{a\sigma }}|_{maximum}`$ $`=`$ $`{\displaystyle \frac{2N}{\widehat{\sigma }^3}}[R(\widehat{s})\widehat{a}s_g^2]=0`$
$`{\displaystyle \frac{^2\mathrm{log}L}{as}}|_{maximum}`$ $`=`$ $`{\displaystyle \frac{NR^{}(\widehat{s})}{\widehat{\sigma }^2}}=0`$
$`{\displaystyle \frac{^2\mathrm{log}L}{\sigma s}}|_{maximum}`$ $`=`$ $`{\displaystyle \frac{2\widehat{a}NR^{}(\widehat{s})}{\widehat{\sigma }^2}}=0,`$
where $`R^{}(s)`$ and $`R^{\prime \prime }(s)`$ are the first and second derivatives of the cross-covariance function. The last two equations follow from the definition of $`\widehat{s}`$ as the shift where the maximum correlation is achieved.

Thus, we see that the Hessian matrix is, fortunately, diagonal and the inversion can be accomplished easily. The resulting squared error is:

$$\sigma _s^2=\frac{\widehat{\sigma }^2}{\widehat{a}NR^{\prime \prime }(\widehat{s})}=\left[N\frac{C^{\prime \prime }(\widehat{s})}{C(\widehat{s})}\frac{C^2(\widehat{s})}{1C^2(\widehat{s})}\right]^1,$$

where $`C^{\prime \prime }(s)`$ is the second derivative of the cross-correlation function.

It is instructive to examine the above expression. It is factored into three separate factors: $`N`$, the number of bins, which has an inverse relation to the error, as is expected intuitively; $`\frac{C^{\prime \prime }(\widehat{s})}{C(\widehat{s})}`$, which is a normalized measure of the sharpness of the cross-correlation peak, and it depends very weakly on the $`S/N`$. A larger absolute value of $`\frac{C^{\prime \prime }(\widehat{s})}{C(\widehat{s})}`$ implies a sharper peak and therefore a smaller error, and vice versa. The third factor, $`\frac{C^2(\widehat{s})}{1C^2(\widehat{s})}`$ is a measure of the “line” $`S/N`$ ratio of the spectrum, i.e., the ratio between the signal standard-deviation and that of the noise. This last statement can be easily seen by the relation of $`1C^2(\widehat{s})`$ to the noise standard deviation.

A different error estimate would result if the template is also taken to be a noisy spectrum, as is sometimes assumed in the literature (e.g., Verschueren & David, 1999) and thus the procedure would have to estimate the ’true’ template using the spectrum and the template. The formal calculations are much more complicated in this case and the resulting error estimate is approximately $`\sqrt{2}`$ times larger than the one quoted above. However, in the usual practice, where the same template is used for many observed spectra, the model presented here is more suitable, since the scatter of the estimated velocities would be influenced only by the noise in the spectra, whereas the noise in the template would contribute a systematic error.

An important feature of maximum-likelihood estimation is that asymptotically (for large $`N`$), it produces a “minimum-variance estimator” (Kendall & Stuart, 1967). This means we cannot expect to achieve an estimate with smaller errors. Since we have shown that the cross-correlation estimate of the shift is equivalent to a maximum-likelihood estimate, this shows the optimality of the cross-correlation method. It is important to bear in mind, though, that the statistical model we used was very simplistic, and other techniques may prove better for removing effects not included in this model, such as spectral mismatch or instrumental long-term trends (e.g., Chelli, 2000).

### 2.4 Simulations

The simulations presented here aim to test the above error estimate under controlled conditions. All the simulations use a spectrum of the G-dwarf HD 38858, obtained by CORALIE, as part of a program to derive precise abundances of planet-hosting and non-planet-hosting stars (Santos et al., 2000, 2001). To test the error estimate a single spectral order was used, ranging from $`\mathrm{5\hspace{0.17em}621}`$ to $`\mathrm{5\hspace{0.17em}683}`$ Å, with $`\mathrm{1\hspace{0.17em}920}`$ bins. Fig. 1 presents this template spectrum. In order to simulate the above specified statistical model, care was taken to eliminate other effects which are not tested in this work. Thus, the assumed doppler shift was always $`0`$, to avoid edge-effects, and the spectrum was rectified to a constant continuum level of unity.

The noise added to the spectra was a simple Gaussian white noise, with standard deviations of $`0.05`$, $`0.1`$ or $`0.2`$. Sample noisy spectra are shown in Fig. 2. These spectra were cross-correlated against the known template, resulting in the cross-correlation functions presented in Fig. 3. The exact peak locations and the relevant derivatives were estimated by interpolating with a parabola around the peak. Each simulation was repeated $`\mathrm{10\hspace{0.17em}000}`$ times. Figure 4 shows the distribution of the ratio between the actual offset of the cross-correlation peak ($`\mathrm{\Delta }v`$) and the error estimate ($`\sigma _v`$). The figure shows that this ratio is distributed as a zero-mean unit-variance Gaussian for all the noise-levels tested, as we would expect for a valid error estimate. The apparent Gaussian nature of the distribution, and especially its unimodality, show the regularity of the measurement error and its estimate. In later sections we will use the same procedure to test the error estimate of the proposed combination scheme and expect to obtain similar distributions.

## 3 Combining cross-correlation functions

### 3.1 The combination scheme

Suppose we now have $`M`$ separate spectra, each one corresponding to a separate template, and all are assumed to have the same Doppler shift. These separate spectra can either be separate exposures of the same object, or separate orders of the same exposure. The following analysis assumes they all possess the same number of bins – $`N`$, although this assumption can be very easily relaxed. We do not assume the same $`S/N`$ ratios for the various spectra nor the same scaling constants.

The most simple scheme to combine cross-correlation functions into one single function is the straightforward, non-weighted average ($`\mathrm{SA}`$):

$$\mathrm{SA}(s)=\frac{1}{M}\underset{i=1}{\overset{M}{}}C_i(s),$$

where $`M`$ is the number of cross-correlation functions to be combined, while $`C_i(s)`$ is the $`i`$-th cross-correlation, calculated between the $`i`$-th spectrum and the $`i`$-th template.

Another simple combination scheme is based on the so-called “determination coefficient”, i.e., the squared correlation coefficient – $`C^2(s)`$ in our case. This quantity is commonly used in regression analysis to describe the degree to which the variability in one variable is explained by the other (e.g., Hays, 1988), and is supposed to be additive. Therefore, we may try to combine these determination coefficients ($`\mathrm{DC}`$) to get:

$$\mathrm{DC}^2(s)=\frac{1}{M}\underset{i}{}C_i^2(s).$$

This value can be interpreted as a weighted sum of the correlations. Each term is weighted by itself, since its value - the correlation - is also a measure of its reliability.

The likelihood approach can provide another scheme to combine cross-correlations. First, let us see what is the overall likelihood of the observations:

$$L=\underset{i}{}\left(\frac{1}{\sqrt{2\pi \sigma _i^2}}\right)^N\mathrm{exp}\left\{\frac{1}{2\sigma _i^2}\underset{n}{}\left[f_i(n)a_ig_i(ns)\right]^2\right\}$$

Converting, as usual to the logarithm, we get:

$$\mathrm{log}L=\underset{i}{}\left\{N\mathrm{log}\sigma _i\frac{1}{2\sigma _i^2}\underset{n}{}\left[f(n)a_ig_i(ns)\right]^2\right\}+\mathrm{𝑐𝑜𝑛𝑠𝑡}.$$

Substituting the optimal values of $`\sigma _i`$ and $`a_i`$ leads eventually to:

$$\mathrm{log}L=\frac{N}{2}\underset{i}{}\mathrm{log}\left[1C_i^2(s)\right]+\mathrm{𝑐𝑜𝑛𝑠𝑡}.$$

which is very similar to eq. 1. Therefore we can equate these two expressions to get an “effective” correlation value – $`\mathrm{ML}`$:

$$NM\mathrm{log}\left[1\mathrm{ML}^2(s)\right]=\underset{i}{}N\mathrm{log}\left[1C_i^2(s)\right]$$

or:

$$\mathrm{ML}^2(s)=1\left\{\underset{i}{}\left[1C_i^2(s)\right]\right\}^{1/M}$$

As it turns out, for sufficiently small $`C_i(s)`$ the above expression reduces to the expression of the “determination coefficient” mean – $`\mathrm{DC}(s)`$. On the other extreme, if one $`C_i(s)`$ approaches unity (meaning it has a very high $`S/N`$), $`\mathrm{DC}(s)`$ will approach unity. This can be understood intuitively, since the other cross-correlations probably represent much poorer $`S/N`$ and therefore can be neglected.

In order to obtain an error estimate for the shift, the Hessian of the likelihood has to be calculated and inverted. This time, the Hessian is not diagonal, but a laborious calculation yields the following expression for the error:

$$\sigma _s^2=\left[MN\frac{\mathrm{ML}^{\prime \prime }(\widehat{s})}{\mathrm{ML}(\widehat{s})}\frac{\mathrm{ML}^2(\widehat{s})}{1\mathrm{ML}^2(\widehat{s})}\right]^1,$$

which is remarkably similar to the expression for the error in the single spectrum case, except for using the total number of bins – $`MN`$ instead of $`N`$!

### 3.2 Simulations

In order to demonstrate the effect of combining cross-correlations, $`10`$ spectral orders were used, taken from the same spectrum used in 2.4, covering the spectral range from $`\mathrm{5\hspace{0.17em}189}`$ to $`\mathrm{5\hspace{0.17em}683}`$ Å. As a first test, noise was added to these $`10`$ orders, with a very high standard deviation of $`3.5`$. Figure 5 presents the results of this test. The top panel shows a cross-correlation function corresponding to one of those orders. No prominent peak is present in the cross-correlation. The middle panel shows a simple average of the ten cross-correlation functions ($`\mathrm{SA}`$). The correct peak is evidently present. It is similarly present in the maximum-likelihood combined correlation – $`\mathrm{ML}`$, as the bottom panel shows. Since the cross-correlation values are very small, $`\mathrm{DC}`$ and $`\mathrm{ML}`$ are virtually identical. This simulation demonstrates the power of combining cross-correlation functions in general. In this case no velocity could have been measured for any individual order, but the cumulative effect of the $`10`$ orders produced a detectable peak.

The next test shows the optimality of the maximum-likelihood scheme ($`\mathrm{ML}`$) compared to the two averages $`\mathrm{SA}`$ and $`\mathrm{DC}`$. This is evident mainly in cases where the $`S/N`$ ratio is strongly varying among the analysed spectra. To show that, noise was added to the $`10`$ orders with standard deviations of $`0.1i`$, where $`i`$ is the order number. Figure 6 shows the distribution of the peak offset for each of the combination schemes, after repeating the simulation for $`1000`$ times. In addition, the Figure also includes a histogram of the velocities derived by a weighted average of the velocities obtained from each order separately ($`\mathrm{VA}`$). The weighting was based on the individual error estimates. Obviously, $`\mathrm{VA}`$ can be calculated only if the correct peaks can be identified in the cross-correlation functions of all orders. The histogram corresponding to $`\mathrm{ML}`$ is evidently less dispersed than all others. Numerically, the four standard deviations corresponding to $`\mathrm{SA}`$, $`\mathrm{DC}`$, $`\mathrm{VA}`$ and $`\mathrm{ML}`$ were $`0.128`$, $`0.078`$, $`0.099`$ and $`0.068`$ $`\mathrm{km}\mathrm{s}^1`$.

An additional criterion to distinguish among combination schemes is their ability to accentuate the correct peak relative to the spurious peaks. This ability is important for detection of faint objects, with low $`S/N`$. To test this ability, the typical height of the spurious peaks was estimated by the standard deviation of the function in a region that is safely distant from the correct peak. Specifically, this standard deviation was estimated on two segments of $`500`$ bins, on both sides of the correct peak, at a distance of at least $`100`$ bins from it. The mean ratios of the peak height to this standard deviation were $`36`$, $`58`$ and $`65`$ for $`\mathrm{SA}`$, $`\mathrm{DC}`$ and $`\mathrm{ML}`$, correspondingly. Once again, we see that $`\mathrm{ML}`$ provides the most prominent correct peak.

The above tests were repeated in numerous other conditions, in all of them $`\mathrm{ML}`$ was found to be the optimal scheme to combine cross-correlations.

After establishing the optimality of the maximum-likelihood combination scheme, its error estimate had to be validated, in a similar fashion to the validation done for the single spectrum case. Figure 7 shows the $`\mathrm{\Delta }v/\sigma _v`$ calculated for the $`\mathrm{ML}`$ functions in the $`1000`$ simulations of $`10`$ orders. Once again, we see a very good agreement with a zero-mean unit-variance Gaussian.

## 4 Concluding remarks

This paper presents a novel approach for combining cross-correlation functions. This approach is based on the relation of cross-correlation to maximum-likelihood estimation. Unlike previously suggested approaches, this approach is not based on a linear combination of the individual functions. Previous works looked for the optimal linear combination, without considering other options (e.g., Connes, 1985). In cases where the scaling constants and the signal-to-noise ratios are known, maximum-likelihood analysis will indeed lead to a linear combination. However, exact estimate of the signal-to-noise is difficult in real cases. The suggested combination does not use any such external estimate.

Although the white gaussian noise model is assumed in most analyses of the cross-correlation, it is obvious that reality is much more complex. Thus, the noise level may depend on the intensity, like in Poissonian noise, or vary across the spectrum because of illumination effects. Additional complexities mainly include long-term trends and spectral mismatch of the template to the analysed spectrum. Attempts to optimize radial-velocity measurement to deal with these problems can be based also on maximum-likelihood theory.

The suggested approach was already applied successfully to the case of HD 41004, where $`30`$ CORALIE orders were combined in order to detect a very faint companion of a K star (Zucker et al., 2003). The combined functions were not the conventional cross-correlation functions, but two-dimensional correlation (TODCOR) functions (Zucker & Mazeh, 1994). The successful detection of the faint companion allowed a more detailed solution of the system and confirmed earlier hypotheses of Santos et al. (2002). This case demonstrates the power of the technique to enhance the capabilities of traditional cross-correlation and improve its ability to detect and measure very faint objects, with very low $`S/N`$.

## acknowledgements

I am grateful to Tsevi Mazeh and Morris Podolak for useful discussions and critical reading of the manuscript. I thank the Geneva planet search group for the CORALIE spectrum of HD 38858. This research was supported by the Israeli Science Foundation (grant no. 40/00) and the Jacob and Riva Damm Foundation.