text:
# Detection of compact objects by means of gravitational lensing in binary systems

## 1 Introduction

One of the most important manifestations of gravitational lensing is the visible variability of the astrophysical objects whose emission is affected by the gravitation field of the lens. This effect becomes observable over a reasonable time period, when the velocities of relative motions of the observer, the lens and the source are great enough, i.e. their mutual location changes rapidly. Many cases of this kind of variability have been examined: from light variation of distant quasars to outburst of stars caused by the influence of planets (see, for instance, Zakharov 1997). The brightness variation of the binary system components caused by gravitational lensing was been first considered by Ingel (1972, 1974) and Maeder (1973). They noted that repetition and small characteristic times of the effect make it one of the most accessible to detection and detailed study. It has also been shown (Maeder 1973) that the large amplitudes of brightness variations ($`0\stackrel{m}{.}10\stackrel{m}{.}5`$) can be expected in binary systems consisting of compact objects – white dwarfs, neutron stars, black holes. In fact, only in these cases does the Einstein-Hvol’son radius turn out to be smaller than the lens star radius, and the image of the radiation source is practically not overlapped by it. The latter circumstance makes the gravitational lensing influence in binary systems an effective means for the detection and detailed investigation of compact objects with the aid of photometric methods alone. Their permeability is by $`2^\mathrm{m}3^\mathrm{m}`$ higher (as compared to spectroscopy), and, therefore, the number of stars amenable to study is an order of magnitude larger. Here the light curve analysis determines all the parameters of the binary system in full analogy to this task for eclipsing binary systems (with allowance made for its singularity) (see, for example, Goncharsky et al. 1985).

A complete set of data for compact object binaries and their statistical analysis gives a possibility to test binary evolution theories and to investigate the last stages of single star evolution as well. We note that only 10 double white dwarfs were detected during the last 10 years (Maxted et al. 2000) and about 40 pairs of a white dwarf and a neutron star – during 25 years (Thorsett & Chakrabarty 1999). We further note that by means of gravitational lensing it is possible to detect compact object binaries at the same rate (see Section 5). Moreover, this may help to discover binaries containing neutron stars with a low magnetic field (not pulsars). Study of these systems in combination with binary star evolution theories gives a possibility to test the detailed cooling models for white dwarfs and neutron stars and the equation of state of the latter (Nelemans et al. 2001; Yakovlev et al. 1999).

Studying compact object pairs is very important to solve a number of astrophysical problems. We present below several such examples.

The ”Standard candle” of modern cosmology, type Ia supernovae, apparently arises from merging of double CO white dwarfs (Webbink 1984; Iben & Tutukov 1984), which have not been detected at this time.

Close double compact objects have to contribute a significant part of the gravitational wave signal at low frequencies. Thus, white dwarf pairs may be a source of the unresolved noise (Evans et al. 1987; Grischuk et al. 2001). Statistical properties of the closest double compact objects in principle determine parameters of the gravitational wave signal.

Binary systems consisting of a white dwarf and a neutron star (a pulsar) are unique laboratories for high-precision tests of general relativity. To date, post-Keplerian general relativistic parameters have been measured in four such pairs (Thorsett & Chakrabarty 1999). This problem may be solved very easily for binaries with gravitational self-lensing since they are observed nearly edge-on.

It seems that the detection of black holes forming pairs with white dwarfs might be an extraordinarily important result of the search for gravitational lensing in compact object binaries. Despite the opinion of most researchers, in a certain sense black holes have not been discovered so far. Only observational data on the behaviour of matter close to the event horizon showing its presence may testify that a black hole is identified (Damour 2000). There is indirect evidence of the presence of black holes in X-ray binaries and galactic cores based on the ”mass-size” relation expected for them (Cherepashchuk 2001). The horizon neighbourhood is seen neither in X-ray binaries nor in the cores of active galaxies because they are screened out by the accreted gas (the accretion rates are very high). This means that only black holes accreting usual interstellar plasma at low rates of $`10^{14}10^{15}M_{}`$/year can be recognized as objects without a surface and with an event horizon (Shvartsman 1971). The black hole companions of white dwarfs in the binaries detected by means of gravitational lensing could be the best objects for horizon study and tests of general relativity in the strong field limit (Damour 2000). It is very easy to measure the mass and size of such a black hole due to binary edge-on orientation and investigate the radiation of gas near the horizon. This matter is accreted from interstellar medium only because its transfer from the white dwarf is absent.

The gravitational lensing effects in binary systems consisting of compact components were studied previously in several papers (Maeder 1973; Gould 1995; Qin et al. 1997). However, the authors did not go further than estimation of the probability of detecting the effect of orientation of the binary system (often underestimating it by a factor of 2-3), which turned out to be very low. Their general conclusion is that the effect cannot be actually recorded. Nevertheless, the data accumulated by now on the evolution of binary systems and their parameters make it possible to define with higher accuracy the probability of detection of the brightness enhancement in binary systems, to find the expected number of such flares and to propose in the final analysis the strategy for their search based on present-day facilities. Our paper is devoted to performing these tasks.

We examine in Section 2 the distinguishing features of light magnification in the systems consisting of two white dwarfs, a white dwarf and a neutron star, and also a white dwarf and a black hole. In Section 3 we derive probabilities of recording the effect for its different amplitudes in all three cases. In Section 4 the numbers of systems of different types which may be detected by means of gravitational lensing are estimated. In Section 5 we discuss the possibilities of the quest for brightness magnification in such systems with the aid of the telescope and equipment used in the survey SDSS (York et al. 2000).

## 2 Light curves in gravitational lensing

We will consider three types of binary systems consisting of
a) two white dwarfs with masses of $`0.7M_{}`$ (sometimes $`0.5M_{})`$,
b) a white dwarf ($`M_{WD}=0.7M_{}`$) and a neutron star ($`M_{NS}=1.4M_{}`$),
c) a white dwarf ($`M_{WD}=0.7M_{}`$) and a black hole ($`M_{BH}=10M_{}`$).
It is clear that both components of a pair play the part of a gravitational lens alternately. However, we will mention at once that the source of radiation will always be the white dwarf to which the approximation of a uniformly luminous disk can be applied. The amplitude, the shape and the duration of flares of the lensed object are defined by the relationships between the parameters of the binary system, the major semiaxis $`a`$, the masses and radii of the components and the orbital plane orientation.

The luminosity magnification of a uniformly radiating disk in gravitational lensing was analysed in detail in the papers by Refsdal (1964), Liebes (1964), Byalko (1969), Ingel (1972, 1974), Maeder (1973) and Agol (2002). On account of the conservation of the surface brightness, the task is reduced in the final analysis to a study of variations of the source image area with allowance made for its occultation by an opaque lens, as the source and lens are moving with respect to the observer (in our case – the orbital motion around the mass centre of the binary system).

Let $`r_l`$ be the distance from the observer to the gravitational lens, $`a`$ is the binary system orbit semiaxis (for the sake of simplicity we consider only circular motions), $`M_l`$ is the lens mass, $`r_g=2GM_l/c^2`$ is its gravitational radius, $`R_l`$ is the lens radius, $`R_s`$ is the source radius, $`M_s`$ is the source mass.

In our case for the light path we can use with high accuracy an approximation of a broken line, the angle between two sections of which is inversely proportional to the impact parameter.

The gravitational lens equation for a point-like source (see, for instance, Zakharov 1997) is

$$\overline{\eta }=\frac{r_l+a}{r_l}\overline{\xi }a\overline{\theta }(\overline{\xi }),$$
(1)
where $`\overline{\eta }`$ and $`\overline{\xi }`$ are the vectors that characterize the location of the points of intersection of the light ray with the planes of the source and the lens, respectively, $`\overline{\theta }(\overline{\xi })`$ is the angle of light deflection (see Fig. 1). In our designations $`\overline{\theta }(\overline{\xi })=\frac{2r_g}{\xi }\frac{\overline{\xi }}{\xi }`$. The main parameter of the task is the Einstein-Hvol’son radius $`r_e=\sqrt{2r_g\frac{ar_l}{r_l+a}}`$, the radius of the ring - the image of the point source, lying on the ”observer – lens” axis in the lens plane. In our case $`r_la`$ and

$$r_e=\sqrt{2r_ga}$$
(2)
or, in solar units,

$$r_e=310^3M^{1/2}a^{1/2}.$$

Equation (1) has two solutions and hence the source images are two points in the lens plane with coordinates:

$$\xi _{1,2}=\frac{d}{2}\left(1\pm \sqrt{1+\frac{4r_e^2}{d^2}}\right),$$
(3)
where $`d=|\overline{\eta }|`$ is the ”lens–source” vector projection onto the lens plane, i.e. displacement of the source off the ”observer–lens” axis.

For the magnification coefficient $`K`$, i.e. the ratio of the sum of brightness of the two images to the point-like source brightness, we have (Zakharov 1997)

$$K=\left|\frac{\xi d\xi }{\eta d\eta }\right|_1+\left|\frac{\xi d\xi }{\eta d\eta }\right|_2=\frac{1}{2}\left(\sqrt{1+\frac{4r_\mathrm{e}^2}{d^2}}+\frac{1}{\sqrt{1+\frac{4r_e^2}{d^2}}}\right).$$
(4)
To evaluate the possibility of using particular approximations, it is necessary to estimate the values of the task parameters. As has already been mentioned, we consider only the white dwarf as a source of radiation. The role of the lens is played by a white dwarf (WD), a neutron star (NS) or a black hole (BH). At typical velocities of the orbital motion, 100–300 km/s, the accretion luminosity of the BH is $`10^310^5`$ as low as the thermal luminosity of the WD (Ipser & Price 1982). Assuming $`M_{WD}=0.7M_{},M_{NS}=1.4M_{},M_{BH}=10M_{},R_{NS}=15,R_{BH}=r_g=2GM_{BH}/c^2`$, we have estimated the rest of the parameters of binary systems, which are presented in Table 1 in solar units.

For the WD radius we used the relationship (Nauenberg 1972):

$$R=0.01125R_{}\left(\left(\frac{M}{1.454M_{}}\right)^{\frac{2}{3}}\left(\frac{M}{1.454M_{}}\right)^{\frac{2}{3}}\right)^{\frac{1}{2}}.$$
(5)
One can see that in a system of two WDs, the Einstein radius, the sizes of the lens and the source are close. This implies that when considering the light magnification one has to take into account the source extension and occultation of its image by the companion – the lens, at least, at $`a<100R_{}`$, which, as we will see, is our primary interest. In binary systems, containing a NS and a BH, eclipsing is absent; the source, however, can be regarded as a point one only at $`a>(1020)R_{}`$, and one can apply expression (4) to derive the brightness increase.

It should be noted that when the binary system is viewed edge-on, at the moment of conjunction, a configuration in which the observer, the lens and the source are on one line $`(d=0)`$ is realized. It was shown (Liebes 1964) that in this case $`K=\sqrt{1+\frac{4r_e^2}{R_s^2}}`$, as a result of averaging of expression (4) over the surface of the uniformly radiating disk. At $`d0`$ the task gets essentially more complicated and has no analytical solution. The accurate relationships for the brightness magnification by a point lens, where elliptical Legendre integrals of the first and second kinds are used, were derived by Refsdal (1964), Liebes (1964), Byalko (1969), Ingel (1972, 1974) and Maeder (1973). There is, however, one more case in which a simple expression for $`K`$ exists – the contact between the disk edge and the lens centre:

$$\begin{array}{cc}K\hfill & =\frac{2}{\pi }[1+4(\frac{r_e}{R_s})^2]\mathrm{arcsin}[1+4(\frac{r_e}{R_s})^2]^{\frac{1}{2}}\hfill \\ & \\ & \frac{2}{\pi }\sqrt{1+4\beta ^2},\hfill \end{array}$$
(6)
where $`\beta =\frac{r_e}{R_s}`$.

The light curve of the component-source is the result of variation of magnification with time variations of $`d`$ during the system rotation period. In the approximation of a point source, according to (4)

$$K(t)=\frac{1}{2}\left[\sqrt{1+\frac{4r_e^2}{d(t)^2}}+\frac{1}{\sqrt{1+\frac{4r_e^2}{d(t)^2}}}\right].$$
(7)
Near the conjunction $`d^2(t)=d_0^2+(vt)^2`$, where $`d_0=a\mathrm{sin}i`$ is the minimum value of $`d`$, which depends on the orbit inclination $`i`$, $`v`$ is the ralative velocity of the system components. At the moment $`t=0`$, the brightness variation reaches a maximum and

$$\begin{array}{cc}K_m\hfill & =K(0)=\frac{1}{2}\left[\sqrt{1+\frac{4r_e^2}{d_0^2}}+\frac{1}{\sqrt{1+\frac{4r_e^2}{d_0^2}}}\right]\hfill \\ & \\ & =\frac{1}{2}\left[\sqrt{1+4\beta _0^2}+\frac{1}{\sqrt{1+4\beta _0^2}}\right],\hfill \end{array}$$
(8)
here $`\beta _0=\frac{r_e}{d_0}`$.

It is customary to characterize the duration of the flare (brightness increase) (see Zakharov & Sazhin 1998) by the time $`\tau `$ that is taken by the moving source projection onto the lens plane to cross the Eintstein ring. It can be easily found from the relation $`d^2(\frac{\tau }{2})=d_0^2+(\frac{v\tau }{2})^2=r_e^2`$, that:

$$\tau =\frac{2r_e}{v}\sqrt{1\frac{1}{\beta _0^2}}.$$
(9)
In a point approximation $`K(\frac{\tau }{2})=\frac{1}{2}[\sqrt{5}+\frac{1}{\sqrt{5}}]1.34`$, and $`\tau `$ defines the flare duration at a level $`\mathrm{\Delta }m_0=0\stackrel{m}{.}32`$ of the source light at a quiescent phase.

Note that the brightness magnification of one of the components by a factor of $`K`$ corresponds to magnification of the system brightness (which is the subject of detection in observations) $`K_0`$ times, and

$$K=K_0+\frac{L_l}{L_s}(K_01),$$
(10)
where $`\frac{L_l}{L_s}`$ is the lens–source brightness ratio. For a system containing two equal white dwarfs $`\frac{L_l}{L_s}=1`$, $`K_0=\frac{1}{2}(K+1)`$, and $`\mathrm{\Delta }m_0=0\stackrel{m}{.}17`$ with respect to the total system brightness off the flare. It is self-evident that if the part of the lens is played by a neutron star or a black hole, the $`K`$ is practically equal to $`K_0`$.

Present-day telescopes detect rather small light variations during short exposures (see Section 5). We are guided by these possibilities when determining the minimum amplitude of detectable flares. Say, if the amplitude of a flare $`\mathrm{\Delta }m_0=0\stackrel{m}{.}2`$ ($`K=1.2`$), then the trajectory of the source projection onto the lens plane does not cross the Eintstein ring at all, and such a flare is excluded from the analysis. Therefore, it is necessary to introduce $`\tau _{0.1}`$ at which $`K=1.1`$. Using (4) we derive

$$(\frac{r_e}{d})^2=\beta ^2=\frac{K^21+K\sqrt{K^21}}{2},$$
(11)
hence at $`K=1.1,\beta _{0.1}=0.60`$, $`d_{0.1}=1.67r_e`$. Thus

$$\tau _{0.1}=1.67\frac{2r_e}{v}\sqrt{1\frac{1}{(1.67\beta _0)^2}}.$$
(12)
In the general case, when deriving the light curves, equation (1) was solved for each point of the source profile, its image was constructed in the lens plane, and its area was found numerically. Since in gravitational lensing, the surface brightness remains constant, the light curve is determined as variation of the ratio of the image areas and the source itself, depending on the position of the binary system components with respect to the observer, which varies during the period of rotation. The contribution to the system radiation of the lens–white dwarf and occultation by it of a part of the source image was taken into account. Fig. 2 shows the brightness variations of the binary system consisting of two white dwarfs with a mass of $`0.7M_{}`$, $`a=14.2R_{}`$ and the orbit inclination $`i=0`$ (i.e. $`d_0=0`$). The similarity of the sizes of the source, the lens and the Einstein ring results in formation of a rather complex light curve. Similar events were analysed by Marsh (2001) and Agol (2002).

Comparison of the light curves in Fig. 3, plotted from the use of the results of direct calculations and with the use of point approximation (7), shows that they are virtually coincident for $`a>30R_{}`$. In this case the light curves can be constructed by means of a point approximation. The exception is the region near the maximum at $`\frac{d(t)}{r_e}<0.25`$.

Consider the effect of gravitational lensing in a binary system comprising two equal white dwarfs with $`M_{WD}=0.5M_{}`$ and $`a=100R_{}`$. According to Kepler’s third law (in solar units):

$$T=2.79a^{3/2}(M_l+M_s)^{\frac{1}{2}}\mathrm{hours},$$
(13)
its period $`T=82.4`$ days. Flares which amplitude depends on $`i`$ will be observed twice during this time. For $`d_0ai>R_s`$, point approximation (8) can be employed to estimate the amplitude. Consequently, at $`i>\frac{R_s}{a}=\frac{10^2}{10^2}=10^4`$, it follows from (2) and (8)

$$K_m=\frac{1}{2}\left[\sqrt{1+\frac{4}{i^2}\frac{2r_g}{a}}+\frac{1}{\sqrt{1+\frac{4}{i^2}\frac{2r_g}{a}}}\right]2.27.$$

Thus the amplitude of these flares $`K_0=\frac{1}{2}(K_m+1)1.64`$ or $`\mathrm{\Delta }m_00.54`$.

From (2) and (9), considering that $`v=\sqrt{\frac{G(M_l+M_s)}{a}}`$, obtain under the same condition (point approximation):

$$\tau =\frac{4a}{c}\sqrt{\frac{M_l}{M_l+M_s}}\sqrt{1\frac{i^2a}{910^6M_l}}.$$
(14)
And for the given system $`\tau =\frac{4a}{c\sqrt{2}}\sqrt{1\frac{i^2a}{4.510^6}}=662\sqrt{12.210^7i^2}`$ sec.

In Fig. 4 the light curves of this system for different $`i`$, plotted from the data of numerical computations, are displayed. The light variations are expressed in stellar magnitudes. The orbit inclination $`i`$ varies from zero to $`2.510^4`$ with a step of $`2.510^5`$ from the upper to the lower curve, which is displaced for convenience. The flare shape variations are due to these changes.

One more step brings us closer to the real condition of searching for flares. Having fixed three levels of brightness increase, $`0\stackrel{m}{.}2`$, $`0\stackrel{m}{.}5`$ and $`1\stackrel{m}{.}0`$, we determined by two methods the minimum sizes of the systems consisting of two white dwarfs and also a white dwarf and a neutron star whose light rises to these levels (Table 2). From the relation $`K=\sqrt{1+\frac{4r_e^2}{R_s^2}}`$, $`a_{min}`$ has been found analytically, while $`a_c`$ has been derived by direct calculations. The light curves of flares in these edge-on-viewed pairs, which were constructed based on numerical computations, are presented in Fig. 5a, b. It can be seen from Table 2 and this figure that the analytical estimates of the characteristic durations of flares $`\tau _{0.1}`$ are close to the results of direct calculations. The flares of a white dwarf coupled with a black hole constructed in a similar manner are exhibited in Fig. 5c (Table 2 contains only $`a_{min}`$ for this binary system), here the separation of the components $`a=3R_{}`$, while the amplitudes $`0\stackrel{m}{.}2,0\stackrel{m}{.}5,1\stackrel{m}{.}0`$ correspond to the orbital plane inclinations 0.023, 0.015 and 0.01. It is precisely this value of $`a`$ that was used in the computation since it is close to the minimum size of the pair WD–BH, which ”survived” in the process of emission of gravitation waves during the lifetime of the Galaxy. The number of closer pairs drops sharply and the probability of their detection is very low (see Section 3).

Fig. 5a, b, c shows that the characteristic durations of flares for pairs of white dwarfs are 50–200 s; for binaries containing NS and BH, the light variations are shorter, 30–130 s and 20–40 s, respectively. Singnificant variations of the flare shape are noticeable in all the cases.

Without dwelling on details, note that the set of parameters of these binary systems obtained in observations – light curves, periods, amplitudes, duration and shape of flares, can make it possible to determine physical characteristics of the components – their masses, sizes, temperatures (e.g. Cherepashchuk & Bogdanov 1995a, b). This problem can be resolved completely for two white dwarfs. At the same time, the determination of, at least, the mass of the lens component and, therefore, its identification as a NS or a BH seems to be self-important for the study of properties of relativistic objects and their extreme gravitation fields.

## 3 Probability of detection of flares and their average characteristics

The number of systems that can be found as a result of brightness increases in gravitation lensing is determined, first of all, by the probability $`F(K_0,t)`$ of recording a flare with an amplitude larger than the level $`K_0`$, during the observing time $`t`$ of a sample of systems distributed in some manner over parameters (orbital separation, masses and component brightness). This means that $`F(K_0,t)`$ is the probability of a flare detection from a certain system averaged over such a distribution

$$F(K_0,t)=\underset{\mathrm{\Omega }}{}P(K_0,t,\omega )f(\omega )𝑑\omega ,$$
(15)
where $`\omega `$ is the set of parameters, $`f(\omega )`$ is the function of their distribution, $`\mathrm{\Omega }`$ is the region in which they are specified, $`P(K_0,t,\omega )`$ is the probability of recording of the flare in a system with the parameters $`\omega `$. For the sake of simplicity, with allowance made for the low accuracy of any estimates associated with statistical characteristics of binary systems with compact companions (Masevich & Tutukov 1988; Lipunov et al. 1996; Nelemans et al. 2001; Fryer et al. 1999, where the errors of different models are defined), we restrict ourselves to an analysis of the distribution of the systems only in $`a`$. Then $`P(K_0,t,a)`$ is the probability of a joint onset of two events, favourable orientation of the orbital plane with the semiaxis $`a`$, when the level $`K_0`$ in the flare is exceeded, and the recording of a result if only one such a flare occurs during the time $`t`$, $`f(\omega )d\omega =f(a)da`$.

The first event is consistent with the geometric probability $`p_g`$ which can be easily found. It follows from (4), (7), (8) and (11) that the values $`d_0^{}<d_0=\frac{r_e}{\beta (K)}`$, where $`K=K_0+\frac{L_l}{L_s}(K_01)`$, correspond to the flares of the source with amplitudes exceeding $`K`$. The locus of the ends of the normals of the orbital planes satisfying this condition is a globe belt of width $`d_0`$ and radius $`a/2`$, its area is $`S_d=\pi ad_0`$. The locus of the ends of the normals of any orbit is the surface of a sphere of radius $`\frac{a}{2}`$, its area is $`S_a=\pi a^2`$. Hence $`p_g=\frac{S_d}{S_a}=\frac{d_0}{a}=\frac{r_e}{a\beta (K)}`$. Finally, from (2) $`p_g=310^3a^{\frac{1}{2}}M_l^{\frac{1}{2}}[\beta (K)]^1`$.

Note that in determining $`\beta (K)`$, both analytical expression (11) and the numerical relationship can be applied. Fig. 6 shows the relationships between geometric probabilities and the size of a system consisting of two white dwarfs which have been derived by these two methods (analytically in a point approximation – the dashed line, and numerically – the solid one). The contribution of the dwarf lens is taken into account here. Minimum sizes of systems capable of giving rise to a flare of amplitude $`K`$ correspond to the values $`a_c`$ in Table 2. The numerical calculations in systems WD–NS and WD–BH are virtually the same as the analytical estimates, but $`p_g`$ for WD-NS reduces to zero at $`a=a_c`$.

The probability $`p_t`$ of recording a flare during the time $`t`$ is essentially dependent on the relationship of the period of recurrence of flares $`T`$ (given by (13)), the duration of the flare $`\tau `$ and the net continuous time of observing $`t`$. Let us lay down certain conditions which must be satisfied to make the search for flares efficient enough. First of all, it is evident that the duration of an elementary continuous exposure $`t_e`$ must be close to that of flare $`\tau `$. On the other hand, the interval between the exposures must be considerably shorter than $`\tau `$. Implying by $`t=\underset{e}{}t_e`$ the duration of the total exposure when observing a specified sky region, we obtain for the probability $`p_t`$:

$$p_t=\{\begin{array}{ccc}\frac{t}{T}& \mathrm{at}& t<T\\ 1& \mathrm{at}& t>T,\end{array}$$
(16)
which corresponds to

$$p_t=\{\begin{array}{ccc}0.36\frac{t(M_l+M_s)^{1/2}}{a^{3/2}}& \mathrm{at}& a>a_t\hfill \\ 1& \mathrm{at}& a<a_t,(16^{})\hfill \end{array}$$

where $`a_t=[0.36t(M_l+M_s)^{1/2}]^{2/3}`$ is the semiaxis of binary with period $`T=t`$ and we consider $`a_ta_c`$ only since $`p_g=0`$ as $`a<a_c`$.

Hence

$$\begin{array}{c}P(K_0,t,a)=p_gp_t\hfill \\ =\{\begin{array}{ccc}0,\hfill & \mathrm{at}\hfill & a<a_c,\hfill \\ 310^3a^{\frac{1}{2}}M_l^{1/2}[\beta (K)]^1\hfill & \mathrm{at}\hfill & a_c<a<a_t,\hfill \\ & & \\ 1.0810^3ta^2\frac{[M_l(M_l+M_s)]^{1/2}}{\beta (K)}\hfill & \mathrm{at}\hfill & a_c<a;\hfill \\ & & a_t<a.\hfill \end{array}\hfill \end{array}$$
(17)
Due to the serious anisotropy of pair distributions on the sky, the optimum time of observation of a specified sky field will be shown to be greater than one night length (at least for WD-WD systems). In this case a simple expression (16) for $`p_t`$ no longer holds since observations on individual nights become statisticaly independent. The probability of detecting at least one flare in $`n`$ nights is given by the Bernoulli distribution:

$$p_n=1(1p_1)^n=\underset{1}{\overset{n}{}}(1)^{i+1}C_n^ip_1^i.$$

where $`p_1=p_t|_{t=9^h}=3.24(M_l+M_s)^{1/2}a^{3/2}`$ is the detection probability in $`9^h`$ (the observational night average length).

Thus, expression for $`P(K_0,n,a)`$ takes the following form

$$\begin{array}{c}P(K_0,n,a)=p_gp_n\hfill \\ =\{\begin{array}{ccc}0,\hfill & \mathrm{at}\hfill & a<a_c\hfill \\ 310^3a^{\frac{1}{2}}M_l^{1/2}[\beta (K)]^1,\hfill & \mathrm{at}\hfill & a_c<a<a_1\hfill \\ 310^3M_l^{1/2}[\beta (K)]^1\hfill & & \\ \times \underset{1}{\overset{n}{}}(1)^{i+1}(3.24)^iC_n^i\hfill & & \\ \times (M_l+M_s)^{i/2}a^{\frac{1}{2}(3i+1)}\hfill & \mathrm{at}\hfill & a_1<a,(17^{})\hfill \end{array}\hfill \end{array}$$

where $`a_1=a_t|_{t=9^h}=2.2(M_l+M_s)^{1/3}`$.

For randomization corresponding to (15) we must convolve $`P(K_0,t,a)`$ or $`P(K_0,n,a)`$ with the pair distribution function in $`a`$.

Now we turn to the analysis of this distribution. We must clearly see the difference in distribution for binary systems, progenitors of pairs of compact objects, and the finite distribution of these pairs themselves. In the first case the results of observations yield a distribution uniform in logarithm of semiaxis (Abt 1983): $`f(a)a^1(1\mathrm{log}(a/R_{})6)`$. In the second case the situation is not so clear because of the uncertainty in quantitative characteristics determined in the course of different versions of population synthesis and because of the qualitative distinctions of these versions (Yungelson et al. 1994; Lipunov et al. 1996; Saffer et al. 1998; Portegies Zwart & Yungelson 1998; Fryer et al. 1999; Nelemans et al. 2001). For instance, in different papers the initial mass functions with different indices are used; due to the great uncertainty in the loss of matter, the supposed masses of BH progenitors vary from $`25M_{}`$ to $`80M_{}`$. The finite numbers of, say, couples WD–BH contain an uncertainty of 5 orders of magnitude (see Fryer et al. 1999)! Nevertheless, since our analysis is qualitative, one can be content with the main characteristics of samples of objects and rough estimates.

First of all, it will be noted that couples lose energy and angular momentum radiating gravitation waves (Landau & Livshits 1983; Grishchuk et al. 2001). Since the rate of these losses depends dramatically on $`a(a^5)`$, the stars in close pairs stick together and leave the sample. Using the expression for gravitational luminosity of a binary system in a circular orbit $`\frac{dE}{dt}=\frac{32}{5}\frac{G^4}{c^5}\frac{M_l^2M_s^2(M_l+M_s)}{a^5}`$ (Landau & Livshits 1983) and taking into account that the total energy of the system $`E=\frac{GM_lM_s}{2a}`$, one can readily derive the value of $`a_r`$ of the initial semiaxis of a system that will merge during $`10^{10}`$ years, which is the life-time of the Galaxy:

$$a_r=2R_{}[M_lM_s(M_l+M_s)]^{1/4}.$$

We tabulate in Table 3 the values of $`a_r`$ and corresponding $`T_r`$ for the systems under study.

The results of population synthesis, which are partially confirmed by observations, suggest that the density destribution of pairs of compact objects rises with decreasing $`a`$ approximately as $`\frac{1}{a}`$, reaches a maximum at $`a=a_0a_r`$, after which it decreases rapidly to $`a1R_{}`$ (Lipunov et al. 1996; Hernanz et al. 1997; Saffer et al. 1998; Portegies Zwart & Yungelson 1998; Wellstein & Langer 1999).

It follows from the same papers that another important parameter, the maximum size of a binary system $`a_{max}`$, is known to be less than $`50R_{}`$ at small eccentricities.

With allowance made for the said above, we have used the function of distribution in $`a`$ in a simple form:

$$\begin{array}{c}f(a)=\gamma \mathrm{lg}eg(a)\hfill \\ =\{\begin{array}{ccc}\frac{\gamma \mathrm{lg}e}{a_0(a_01)}(a1),\hfill & \mathrm{at}& 1aa_0,\hfill \\ & & \\ \frac{\gamma \mathrm{lg}e}{a},\hfill & \mathrm{at}& a_0aa_{max}.\hfill \end{array}\hfill \end{array}$$
(18)
From the condition of normalizing $`\underset{1}{\overset{a_{max}}{}}f(a)𝑑a=1`$ one can easily find $`\gamma `$:

$$0.43\gamma =\frac{1}{\frac{a_01}{2a_0}+\mathrm{ln}\frac{a_{max}}{a_0}}.$$
(19)
Note at once that the share of systems with $`a<a_0`$ ranges from 10 % to 30 % at $`3<a_0<10`$ and $`a_{max}=50,30`$, which is close to the results of population synthesis for pairs of different type compact objects (Hernanz et al. 1997; Saffer et al. 1998; Portegiez Zwart & Yungelson 1998; Nelemans et al. 2001).

In our further estimations we will use distribution (18), varying the parameter $`a_0`$ in the range $`3R_{}÷10R_{}`$, at two values of the parameter $`a_{max}30R_{}`$ and $`50R_{}`$.

Having fixed the distribution $`f(a)`$, now we have to determine the lower limits of integrating in (15) for different types of binary systems. Referring to Table 2, Table 3 and Fig. 6 and allowing for (17) and (18), it can be readily seen that the lower limits for systems consisting of two white dwarfs and a white dwarf in pair with a neutron star coincide with $`a_c`$ and are 14 and 2.5, respectively. However, for pairs consisting of a white dwarf and a black hole integration is done over the whole range of possible values of $`a`$, i.e. from 1 to $`a_{max}`$. The values of binary periods corresponding to $`a_c`$ are 123, 7.6 and 0.85 hours for WD-WD, WD-NS and WD-BH systems respectively. This leads to the following expressions for flare detection probability obtained from (16), ($`16^{}`$), (17), ($`17^{}`$) and (18).

For WD-WD pairs
when $`t<9^h<T_c=123^h`$ since $`a|_{t=9^h}<a_c`$ and $`a_0<a_c`$,

$$\begin{array}{c}F(K_0,t)=\underset{a_c}{\overset{a_{max}}{}}f(a)1.0810^3[\beta (K)]^1ta^2\hfill \\ \times [M_l(M_l+M_s)]^{1/2}da=1.0810^30.43\gamma \hfill \\ \times [\beta (K)]^1t[M_l(M_l+M_s)]^{1/2}\underset{a_c}{\overset{a_{max}}{}}a^3𝑑a\hfill \\ =5.410^4[\beta (K)]^1t[M_l(M_l+M_s)]^{1/2}\hfill \\ \times \left[\frac{a_01}{2a_0}+\mathrm{ln}\frac{a_{max}}{a_0}\right]^1\left(\frac{1}{a_c^2}\frac{1}{a_{max}^2}\right);\hfill \end{array}$$
(20)
when the overall observation time of a specified sky region is $`n`$ nights,

$$\begin{array}{c}F(K_0,n)=\underset{a_c}{\overset{a_{max}}{}}f(a)310^3M_l^{1/2}[\beta (K)]^1\hfill \\ \times \underset{i=1}{\overset{n}{}}(1)^{i+1}(3.24)^iC_n^i(M_l+M_s)^{\frac{i}{2}}a^{\frac{3i+1}{2}}da\hfill \\ =310^30.43\gamma M_l^{1/2}[\beta (K)]^1\hfill \\ \times \underset{i=1}{\overset{n}{}}(1)^{i+1}(3.24)^iC_n^i(M_l+M_s)^{\frac{i}{2}}\underset{a_c}{\overset{a_{max}}{}}a^{\frac{3(i+1)}{2}}da\hfill \\ =310^3\left[\beta (K)\left(\frac{a_01}{2a_0}+\mathrm{ln}\frac{a_{max}}{a_0}\right)\right]^1\hfill \\ \times M_l^{1/2}\underset{i=1}{\overset{n}{}}(1)^{i+1}(3.24)^iC_n^i(M_l+M_s)^{\frac{i}{2}}\frac{2}{3i+1}\hfill \\ \times [a_c^{\frac{3i+1}{2}}a_{max}^{\frac{3i+1}{2}}].(20^{})\hfill \end{array}$$

For WD-NS pairs
when $`t<T_c=7.^h6`$, since $`a_t<a_c<a_0`$

$$\begin{array}{c}F(K_0,t)=1.0810^3[\beta (K)]^1t[M_l(M_l+M_s)]^{1/2}\hfill \\ \times \underset{a_c}{\overset{a_{max}}{}}f(a)a^2da=1.0810^3\hfill \\ \times \left[\beta (K)\left(\frac{a_01}{2a_0}+\mathrm{ln}\frac{a_{max}}{a_0}\right)\right]^1t[M_l(M_l+M_s)]^{1/2}\hfill \\ \times \left[\frac{1}{a_0(a_01)}\left[\mathrm{ln}\frac{a_0}{a_c}\left(\frac{1}{a_c}\frac{1}{a_0}\right)\right]+\frac{1}{2}\left(\frac{1}{a_c^2}\frac{1}{a_{max}^2}\right)\right];\hfill \end{array}$$
(21)
when $`7.^h6<t<9^h`$, since $`a_c<a_t<a_0`$

$$\begin{array}{c}F(K_0,t)=310^3[\beta (K)]^1M_l^{1/2}[\underset{a_c}{\overset{a_t}{}}a^{1/2}f(a)da\hfill \\ +0.36(M_l+M_s)^{1/2}t\underset{a_t}{\overset{a_{max}}{}}a^2f(a)da]\hfill \\ =310^3\left[\beta (K)\left(\frac{a_01}{2a_0}+\mathrm{ln}\frac{a_{max}}{a_0}\right)\right]^1M_l^{1/2}\hfill \\ \times \{\frac{\frac{2}{3}[0.36t(M_l+M_s)^{1/2}a_c^{3/2}]2[0.71t^{1/3}(M_l+M_s)^{1/6}a_c^{1/2}]}{a_0(a_01)}\hfill \\ +0.36t(M_l+M_s)^{1/2}[\frac{1}{a_0(a_01)}[\mathrm{ln}a_0\frac{2}{3}\mathrm{ln}(0.36t\hfill \\ \times (M_l+M_s)^{1/2})\left(0.36t(M_l+M_s)^{1/2}\right)^{\frac{2}{3}}+\frac{1}{a_0}]\hfill \\ +\frac{1}{2}(\frac{1}{a_0^2}\frac{1}{a_{max}^2})]\};(21a)\hfill \end{array}$$

when the overall observation time of a specified sky region is $`n`$ nights, since $`a_c=2.5<a_t|_{t=9^h}=2.8`$ and $`a_t<a_0`$

$$\begin{array}{c}F(K_0,n)=310^3[\beta (K)]^1M_l^{1/2}[\underset{a_c}{\overset{a_t}{}}f(a)a^{1/2}da\hfill \\ +\underset{i=1}{\overset{n}{}}(1)^{i+1}(3.24)^iC_n^i(M_l+M_s)^{\frac{i}{2}}\underset{a_t}{\overset{a_{max}}{}}f(a)a^{\frac{3i+1}{2}}da]\hfill \\ =310^3[\beta (K)(\frac{a_01}{2a_0}+\mathrm{ln}\frac{a_{max}}{a_0})\hfill \\ \times (a_0(a_01))\left]^1M_l^{1/2}\right[\frac{2}{3}[3.24(M_l+M_s)^{1/2}a_c^{3/2}]\hfill \\ 2[1.48(M_l+M_s)^{1/6}a_c^{1/2}]+3.24n(M_l+M_s)^{1/2}\hfill \\ \times \mathrm{ln}\frac{a_0}{2.2(M_l+M_s)^{1/3}}+\underset{i=2}{\overset{n}{}}(1)^{i+1}(3.24)^iC_n^i(M_l+M_s)^{\frac{i}{2}}\hfill \\ \times \frac{2}{3(i1)}\left[(3.24(M_l+M_s)^{1/2})^{(i1)}a_0^{\frac{3(i1)}{2}}\right]\hfill \\ \underset{i=1}{\overset{n}{}}(1)^{i+1}(3.24)^iC_n^i(M_l+M_s)^{\frac{i}{2}}\hfill \\ \times [\frac{2}{3i1}[(3.24(M_l+M_s)^{1/2})^{(i\frac{1}{3})}a_0^{\frac{3i1}{2}}]\hfill \\ \frac{2}{3i+1}(a_0^{\frac{3i+1}{2}}a_{max}^{\frac{3i+1}{2}})]].(21^{})\hfill \end{array}$$

For WD-BH pairs
when $`t<T_c=0.^h85`$, since $`a_t<a_c<a`$ and $`a_c<a_0`$ where $`a_c=1`$ (left distribution edge)

$$\begin{array}{c}F(K_0,t)=1.0810^3[\beta (K)]^1t[M_l(M_l+M_s)]^{1/2}\hfill \\ \times \underset{1}{\overset{a_{max}}{}}f(a)a^2da=1.0810^3[\beta (K)(\frac{a_01}{2a_0}\hfill \\ +\mathrm{ln}\frac{a_{max}}{a_0})]^1t[M_l(M_l+M_s)]^{1/2}\hfill \\ \times \left[\frac{1}{a_0(a_01)}\left(\mathrm{ln}a_01+\frac{1}{a_0}\right)+\frac{1}{2}\left(\frac{1}{a_0^2}\frac{1}{a_{max}^2}\right)\right];\hfill \end{array}$$
(22)
when $`0.^h85<t<9^h`$, since $`1<a_t<a_t|_{t=9^h}=4.8`$

$$\begin{array}{cc}F(K_0,t)=\hfill & 310^3[\beta (K)]^1M_l^{1/2}[\underset{1}{\overset{a_t}{}}a^{1/2}f(a)da\hfill \\ & +0.36(M_l+M_s)^{1/2}t\underset{a_t}{\overset{a_{max}}{}}a^2f(a)da],\hfill \end{array}$$

and for $`a_t<a_0`$

$$\begin{array}{c}F(K_0,t)=310^3\left[\beta (K)\left(\frac{a_01}{2a_0}+\mathrm{ln}\frac{a_{max}}{a_0}\right)\right]^1\hfill \\ \times M_l^{1/2}[\frac{0.24t(M_l+M_s)^{1/2}1.42t^{1/3}(M_l+M_s)^{1/6}+\frac{4}{3}}{a_0(a_01)}\hfill \\ +0.36t(M_l+M_s)^{1/2}[\frac{1}{a_0(a_01)}[\mathrm{ln}a_0\hfill \\ \frac{2}{3}\mathrm{ln}[0.36t(M_l+M_s)^{1/2}]2.8t^{2/3}(M_l+M_s)^{1/3}]\hfill \\ +\frac{1}{2}(a_0^2a_{max}^2)]],(22a)\hfill \end{array}$$

or for $`a_0<a_t`$

$$\begin{array}{c}F(K_0,t)=310^3\left[\beta (K)\left(\frac{a_01}{2a_0}+\mathrm{ln}\frac{a_{max}}{a_0}\right)\right]^1\hfill \\ \times M_l^{1/2}\{\frac{1}{a_0(a_01)}(\frac{2}{3}a_0^{3/2}2a_0^{1/2}+\frac{4}{3})+2[a_0^{1/2}\hfill \\ 1.41t^{\frac{1}{3}}(M_l+M_s)^{\frac{1}{6}}]+0.18(M_l+M_s)^{1/2}t\hfill \\ \times [3.8t^{\frac{4}{3}}(M_l+M_s)^{\frac{2}{3}}a_{max}^2]\};(22b)\hfill \end{array}$$

when the overall observation time of a specified sky region is $`n`$ nights, for $`a_1=2.2(M_l+M_s)^{1/3}=4.8<a_0`$, $`F(K_0,n)`$ coincides with ($`21^{}`$) where $`a_c`$ is replaced by 1, and for $`a_04.8`$

$$\begin{array}{c}F(K_0,n)=310^3M_l^{1/2}[\beta (K)]^1(\underset{1}{\overset{a_1}{}}f(a)a^{1/2}da\hfill \\ +\underset{i=1}{\overset{n}{}}(1)^{i+1}(3.24)^i(M_l+M_s)^{\frac{i}{2}}C_n^i\underset{a_1}{\overset{a_{max}}{}}a^{\frac{3i+1}{2}}f(a)da)\hfill \\ =310^3\left[\beta (K)\left(\frac{a_01}{2a_0}+\mathrm{ln}\frac{a_{max}}{a_0}\right)\right]^1\hfill \\ \times M_l^{1/2}\{\frac{1}{a_0(a_01)}(\frac{2}{3}a_0^{3/2}2a_0^{1/2}+\frac{4}{3})\hfill \\ +2a_0^{1/2}1.35(M_l+M_s)^{1/6}\hfill \\ +\underset{i=1}{\overset{n}{}}(1)^{i+1}(3.24)^i(M_l+M_s)^{\frac{i}{2}}C_n^i\frac{2}{3i+1}\hfill \\ \times [(2.2)^{\frac{3i+1}{2}}(M_l+M_s)^{\frac{3i+1}{6}}a_{max}^{\frac{3i+1}{2}}]\}.(22^{})\hfill \end{array}$$

We use the value of $`0\stackrel{m}{.}2`$ for $`K_0`$ for all systems since it will be shown in Section 5 that flares with such amplitudes are detectable with the present-day observing facilities.

The values of $`F(K_0,t)`$ for $`K_0=1.2(\mathrm{\Delta }m=0\stackrel{m}{.}2)`$ and $`t=9^h`$ are presented in Table 4 for a range of $`a_0`$ and $`a_{max}`$ values. These rough model estimates may change by a factor of 1.5–2 depending on the real pattern of the distribution of pairs of compact objects along the semiaxes.

Fig.7 represents the behaviour of $`F(t)F(K_0,t)`$ (for $`\mathrm{\Delta }m=0\stackrel{m}{.}2)`$ as a function of $`t`$ and $`n`$ for the values of $`a_0`$ and $`a_{max}`$ which give approximately the average values among those presented in Table 4 for different system types.

We have also estimated the average durations and periods of reccurence of the events considered and their dispersions using (20) – ($`22^{}`$). To do this correctly one has to convolve $`\tau `$ and $`T`$ not with the $`f(a)`$ but with the $`P(K_0,t,a)f(a)/F(K_0,t)`$ according to the share each parameter region brings to the probability of detecting a flare. These values, formally, depend on the values of $`t`$ or $`n`$. However, it is clear that their variations cannot influence much the real average characteristic. Actually in these calculations we use the expression for $`P(K_0,t,a)`$ and $`F(K_0,t)`$ in the region where they are linear in $`t`$, which as we will see in Section 5 corresponds to the optimal strategy for searching for the flares.

The results are presented in Table 5 for the value of amplitude $`\mathrm{\Delta }m0\stackrel{m}{.}2`$. The parameters averaged over the range of the used $`a_0`$ (from $`3R_{}`$ to $`10R_{}`$) are in the first and third columns. Variations of flare duration with varying $`a_0`$ do not exceed a factor of 2. In columns 2 and 4 the limits of the intervals which contain parameters of 90 % of flares are given. The minimum values of the parameters correspond to the closest systems of all types which are capable of giving rise to a flare with the amplitude $`\mathrm{\Delta }m0\stackrel{m}{.}2`$. In Table 5 we give parameters of events that one can hope to discover in observations.

## 4 The number of pairs being investigated in the Galaxy

As we have already noted, the present-day numbers of different pairs of compact objects are determined in the framework of different versions of population synthesis (Lipunov et al. 1996; Wellstein & Langer 1999; Fryer et al. 1999; Nelemans et al. 2001). However, the main objective of this research is the study of systems which are of interest within the scope of gravitational-wave astronomy and of the problem of origin of gamma-ray bursts. These are pairs of neutron stars, neutron stars with black holes and essentially rarely white dwarfs with black holes. By now extensive work has been done on the modeling of the origin and evolution of pairs of white dwarfs. It goes without saying that these binary systems are best understood; comparison with observations can be made for them (Hernanz et al. 1997; Nelemans et al. 2001). There are at least two relatively reliable estimates of the total number of binary white dwarfs in the Galaxy, $`510^8`$ (Lipunov et al. 1996) and $`2.510^8`$ (Nelemans et al. 2001). In principle, these data are sufficient to determine the number of gravitational-lens flares that can be detected in this sample. Nevertheless, we have estimated the number of binary white dwarfs and also white dwarfs in pairs with neutron stars and black holes in the Galaxy in the framework of assumptions of the progenitors of binary systems consisting of compact objects. We have followed Bethe and Brown (1998, 1999) who have determined the number of massive binaries containing neutron stars and black holes.

We will base our discussion on the following assumptions (as previously, solar units are used).

1. The mass distribution of progenitor stars (initial mass function) obeys Salpeter’s (1955) law in the interval $`0.1<M<120`$, i.e.

$$dW(M)=w(M)dM=1.35(\frac{M}{0.1})^{2.35}dM,$$
(23)
where $`dW(M)`$ is the probability that the mass of a star will be in the interval $`M÷M+dM`$, $`w(M)`$ is the probability density.

At the observed star formation rate $`\frac{dN}{dMdt}=0.9(\frac{M}{0.1})^{2.35}\mathrm{year}^1`$ and the Galaxy age of $`10^{10}`$ years, the number of stars is $`N1.510^{11}`$.

2. About 100 % of stars are members of binary systems. This evaluation is used in the population synthesis as well as 50 % (see for instance, Nelemans et al. 2001).

3. The distribution density in the mass ratio, $`qM_s/M_l1`$, is constant, i.e. $`\frac{dN}{dq}\mathrm{const}`$ ( e.g. Portegies Zwart & Yungelson 1998).

4. The initial masses of progenitor stars must lie within the following intervals to give:
white dwarf – $`1<M<10`$, neutron star – $`10<M<25`$, black hole – $`25<M`$.
In the last case the uncertainty is rather great; the values $`40M_{}`$ (Van den Heuvel & Habets 1984), and $`60M_{}`$ (Woosly et al. 1995) are also used. As for us, we follow Portegies Zwart et al. (1997).

5. The neutron star formed in a supernova explosion in a binary system receives a kick. As Cordes & Chernoff (1997) have shown, its distribution is well approximated by the sum of two Gaussians with standard 175 km/s and 700 km/s, the share of the former being 80 % while that of the latter being 20 %. With these parameters 43 % of binary systems survive after the explosion (Bethe & Brown 1998). As a black hole is formed, one can believe that the kick will be lower, at least, inversely proportional to its mass (Fryer et al. 1999). For this case, using the relationships between the probability of decay and the additional velocity derived by Bethe & Brown (1998), we have found that 87 % of binaries survive after the formation of a black hole.

Now we turn to estimation of the number of pairs of different types.

For the probability $`d\varphi `$ to obtain a binary system with the mass of the more massive companion (in our case the object is a lens) in the interval $`M_l÷M_l+dM_l`$ and the mass ratio in the interval $`q÷q+dq`$, considering (23) and the homogeneity of the distribution in $`q`$ ($`M_l`$ is measured in $`0.1M_{}`$), we have

$$d\varphi =dqdW=1.35(M_l)^{2.35}dqdM_l.$$
(24)
The number of pairs of type $`i`$ is $`N_i=\varphi _i\gamma _iN_{bin}`$, $`i=d`$ for a pair of white dwarfs, $`i=n`$ and $`i=b`$ for binaries containing a white dwarf with a neutron star and a black hole, respectively; $`\gamma _i`$ is the pair survival probability after the supernova explosion.

Clearly

$$\varphi _i=\underset{Q_i}{}\underset{M_i}{}𝑑\varphi =1.35\underset{Q_i}{}\underset{M_i}{}(M_l)^{2.35}𝑑q𝑑M_l,$$
(25)
where $`Q_i`$ is the interval where the values of $`q`$ lie, and $`M`$ is the interval of $`M_l`$ values. Determine $`Q_i`$ for each type of pairs, taking into account that $`10<M_s=qM_l<100`$ (assumption 5). Then $`\frac{10}{M_l}<q<1`$ for a pair WD–WD, $`\frac{10}{M_l}<q<\frac{100}{M_l}`$ for a pair WD–NS and WD–BH. Thus, proceeding from (25), we obtain for a pair of white dwarfs

$$\begin{array}{cc}\varphi _d\hfill & =1.35\underset{10}{\overset{100}{}}(M_l)^{2.35}𝑑M_l\underset{\frac{10}{M_l}}{\overset{1}{}}𝑑q\hfill \\ & =1.35\underset{10}{\overset{100}{}}(M_l)^{2.35}\left(1\frac{10}{M_l}\right)𝑑M_l.\hfill \end{array}$$
(26)
Hence

$$\begin{array}{cc}\varphi _d\hfill & =1.35\underset{10}{\overset{100}{}}M_l^{2.35}𝑑M_l13.35\underset{10}{\overset{100}{}}M_l^{3.35}𝑑M_l\hfill \\ & =0.0430.025=0.018,\hfill \end{array}$$

and allowing for the fact that the number of pairs in the Galaxy is $`N_{bin}=0.5N`$, obtain for the number of binary white dwarfs $`N_d=1.3510^9`$ pairs.

For white dwarfs in pair with neutron stars

$$\begin{array}{cc}\varphi _n\hfill & =1.35\underset{100}{\overset{250}{}}M_l^{2.35}𝑑M_l\underset{\frac{10}{M_l}}{\overset{\frac{100}{M_l}}{}}𝑑q\hfill \\ & =1.35\underset{100}{\overset{250}{}}M_l^{2.35}\frac{90}{M_l}𝑑M_l=901.35\underset{100}{\overset{250}{}}M_l^{3.35}𝑑M_l,\hfill \end{array}$$

and finally

$$\varphi _n=\frac{901.35}{2.35}\left[M_l^{2.35}|_{100}^{250}\right]=51.71.7610^5=9.110^4.$$

Because 57 % of binary systems decay after the supernova explosion which results in formation of a neutron star, $`N_n9.110^40.430.5N2.910^7`$ pairs.

The same is for pairs containing black holes:

$$\varphi _b=51.7\left[M_l^{2.35}|_{250}^{1200}\right]=51.72.2610^6=1.210^4.$$

Since the probability of their survival is 87%, then $`N_b1.210^40.870.5N7.810^6`$ pairs.

It should be emphasized that the above-given estimates are defined in many respects by the choice of the power index $`k`$ in the initial function of stellar masses. Clearly they decrease with increasing $`k`$. It is likely that the version with $`k=2.35`$ of Salpeter (1955) is most consistent with current observational data (see, for instance, Grishchuk et al. 2001; Raguzova 2001), but, nevertheless, $`k=2.5`$ and also $`k=2.7`$ are used in the population synthesis as well (Bethe & Brown 1998, 1999; Nelemans et al. 2001).

To make the picture complete, we present in Table 6 the estimates of the numbers of pairs in the Galaxy for the mentioned $`k`$.

On the other hand, in order to determine the accuracy of our estimates, one can compare them (for binary white dwarfs) with the results of population synthesis. In particular, the number of pairs of this type is $`4.510^8`$ (Lipunov et al. 1996) and $`2.510^8`$ (Nelemans et al. 2001). Note that the last number was obtained for $`k=2.7`$ and it is 2 times as low as our rough value from Table 6. The 2-4–time differences are reasonable, taking into account the uncertainty in the results of modeling of evolution and the qualitative character of our estimates.

## 5 Strategy of search and the number of detectable objects. Conclusions

To determine the number of objects that may be detected using a certain telescope one has to know not only the overall number but also their distribution in the Galaxy, their luminosity function and the distribution of absorbing matter as well.

We used the star distribution in the Galaxy of the following form:

$$\rho =\rho _0e^{\frac{R}{H}}e^{\frac{|z|}{h}},$$
(27)
where $`R`$ and $`z`$ are the point cylindrical coordinates in the Galaxy, $`H`$ and $`h`$ are radial and vertical distance scales, respectively, for which we used the following values: $`H=8`$ kpc and $`h=250`$ pc (see Dehnen & Binney 1998; Nelemans et al. 2001). Adopting the Sun’s coordinates $`R_{}=8.5`$ kpc and $`z_{}=30`$ pc we have:

$$\rho _0=\rho _{loc}e^{\frac{R_{}}{H}}e^{\frac{z_{}}{h}}.$$

The overall number of objects in the Galaxy is expressed by the integral

$$N=\underset{\mathrm{}}{\overset{\mathrm{}}{}}𝑑z\underset{0}{\overset{\mathrm{}}{}}2\pi R\rho _0e^{\frac{R}{H}}e^{\frac{|z|}{h}}𝑑R,$$

which yields

$$\rho _{loc}=\frac{Ne^{\frac{R_{}}{H}}e^{\frac{|z_{}|}{h}}}{4\pi H^2h}=\frac{N}{6.5610^{11}}\mathrm{pc}^3$$
(28)
and

$$\rho (R,z)=\frac{N}{6.5610^{11}}e^{\frac{RR_{}}{H}}e^{\frac{|z|z_{}}{h}}\mathrm{pc}^3.$$
(29)
We used distributions for WD-NS and WD-BH in the same form with the only difference in normalization coming from their fewer overall number in the Milky Way. The local densities corresponding to the overall numbers given in Table 6 are presented in Table 7.

Note that the part of WD-WD pairs among the complete number of WDs is (7-25) % (Iben et al. 1997; Fryer et al. 1999) and the WD density in the Solar neighbourhood is $`(420)10^3`$ pc<sup>-3</sup> (see Nelemans et al. 2001 and references therein). This gives us a range of local density of $`310^4510^3`$ pc<sup>-3</sup>. It is very close to our estimate shown in Table 7 and the most optimistic values of density may be $`22.5`$ times more than those.

Unfortunately, there is no data to compare with for binaries other than WD-WD. We suggest only that our estimation range may be $`23`$ times narrower. Therefore, the real number of objects may be a few times higher.

The insterstellar absorption is proportional to the interstellar medium density which is distributed in the Galaxy by the same (27) law with $`h=80`$ pc (Dehnen & Binney 1998). The normalization was done according to the local average value of absorption $`A_g^{}=1\stackrel{m}{.}16\mathrm{k}\mathrm{p}\mathrm{c}^1`$ in the $`g^{}`$ band of the SDSS equipment (Schlegel et al. 1998).

The WD luminosity function is well known and is explained fairly well in the frame of their cooling theory. We used the luminosity function of Oswalt et al. (1996). It is readily seen from this that more than half of the WDs are brighter than $`13\stackrel{m}{.}5`$ and only less than 20 percent are fainter than $`15^\mathrm{m}`$, where the function reaches its maximum. Thus, the magnitudes of the majority of WDs do not suffer much from the bolometrical corrections which is a few tenth for the stars of these spectral classes (see, e.g. Allen 1973).

In other binary systems, the situation is more complicated. It is only clear that the optical emission from both a NS and a BH can be associated with interstellar gas accretion. In the former case it is likely to be negligibly small because of the high orbital velocity of the NS ($`>100`$ km/s) (Shapiro & Teukolsky 1983). At the same time the orbital velocity of a BH in a pair with a WD is about 50 km/s, and the luminosity of accretion plasma can reach $`10^{30}`$ erg/s by moderate-optimistic estimates (Shvartsman 1971; Beskin & Karpov 2002) and will be by an order of magnitude less by a pessimistic one (Ipser & Price 1982). Apparently, its contribution to the total luminosity of the system is insufficient to correct the depth of the space of detection. However, this level of optical emission makes possible searching for its fast variations with a microsecond time resolution for investigation of accretion processes and strong gravitational fields near the horizon of events (Shvartsman 1971; Beskin et al. 1997; Beskin & Karpov 2002).

Now we have everything we need to calculate the on-sky density of objects in a given direction ($`l,b`$):

$$\begin{array}{cc}\sigma (l,b)\hfill & =\frac{d^2N}{\mathrm{cos}bdldb}\hfill \\ & =r^2\rho (R,z)I(M)𝑑r,\hfill \end{array}$$
(30)
where $`I(M)`$ is a cumulative luminosity function from Oswald et al. (1996), $`M=m_{lim}5\mathrm{lg}(r/10)A(r,l,b`$), $`m_{lim}`$ is an equipment limiting apparent stellar magnitude, $`A(r,l,b)`$ is a total extinction in the $`(l,b)`$ direction up to the distance $`r`$, and

$$\begin{array}{c}z=z_{}+r\mathrm{sin}b,\hfill \\ R=\sqrt{R_{}^2+r^2\mathrm{cos}^2l2R_{}r\mathrm{cos}b\mathrm{cos}l}\hfill \end{array}$$

are the point cylindrical coordinates expressed in terms of its galactic coordinates $`(l,b)`$ and distance $`r`$.

Unfortunately, this integration cannot be done analytically and so we will turn at once to a particular version of the equipment being discussed. In our view, an ideal tool for carrying out the proposed programme is currently the 2.5-metre telescope at APO (New Mexico). It has been used to accomplish one of the most promising projects – the Sloan Digital Sky Survey (SDSS). The telescope has a field of $`3^{}`$, in which 30 CCD chips of $`2048\times 2048`$ pixels are mounted (York et al. 2000). Sky scanning at a speed of its movement along a strip $`2\stackrel{}{.}5`$ wide is performed within the framework of the project. During an exposure of 55 s, the limiting stellar magnitude in the $`g^{}`$ filter (it is close to the $`B`$ filter of the Johnson system) is $`23\stackrel{m}{.}2`$ at a signal-to-noise ratio of $`5`$ (Oke & Gunn 1983; York et al. 2000). We have constructed a relationship between amplitude of flares recorded at a $`5\sigma `$ level and stellar magnitude of a quiet object in the $`g^{}`$ band (Fig. 8). In so doing, we used the parameters corresponding to the SDSS project: the seeing – 0.8 arcsec, the quantum efficiency of the CCD matrix – 80 %.

By comparing the data of Fig. 8 with the mean parameters of flares (Table 5) and their light curves for different binary systems (Fig. 5), it is clear that with a particular exposure of time $`t_e1030`$ s, the flares with $`\mathrm{\Delta }m0\stackrel{m}{.}2`$ of any pairs that we have examined can be recorded at a significance level of $`510^310^5`$ (duration of any flare is a few $`t_e`$). In fact, the level of the flare detection will be better due to registration in the SDSS of the field in several colour bands ($`u^{},g^{},r^{},z^{}`$) and the limiting $`r^{}`$ stellar magnitude is also $`23\stackrel{m}{.}`$ Moreover, the use and comparison of the data in the $`g^{}`$ and $`r^{}`$ filters give the possibility of easy separation of cosmic ray traces.

On the basis of (30) we have built a map of the expected on-sky densities of WD-WD pairs up to $`m_{lim}=23^\mathrm{m}`$ which is presented in Fig. 9. The value of the local density corresponding to the Salpeter mass-spectrum index of $`2.35`$ is used. One can see from the map that the maximal density regions have galactic latitudes of $`510`$ degrees and a question might arise if it is possible to observe in such star-rich regions of the sky. However, even in these latitudes the average distance between stars of up to $`23^\mathrm{m}`$ is about 3-4 arcseconds (see e.g., Zombeck, 1982, p.34) and thus, they are resolved with the SDSS seeing and this crowding therefore cannot have much negative impact on observations with this equipment.

To derive the expected number of events of the discussed type $`N_{exp}`$ one multiplies the probability of detecting it from a randomly chosen pair $`F(t)`$ for $`\mathrm{\Delta }m=0\stackrel{m}{.}2`$ (see (15), (20)-(22) and Fig.7) by the number $`D(\mathrm{\Omega })`$ of such pairs observed during a certain observational programme that covered the celestial sphere part $`\mathrm{\Omega }`$:

$$N_{exp}=F(t)\times D(\mathrm{\Omega }),$$
(31)
where

$$D(\mathrm{\Omega })=\underset{\mathrm{\Omega }}{}\sigma (l,b)\mathrm{cos}bdldb.$$
(32)
It is clear, however, that when applying to reality one deals with a telescope of limited field of view and the observational time provided for a programme is limited as well, so increasing the total exposure of a field to increase $`F(t)`$ leads to decreasing the number of such fields in the programme. It is also evident that one has to monitor at first the best, i.e. maximal density $`\sigma (l,b)`$ fields first. Fig. 10 represents the number of pairs of each type $`D(\mathrm{\Omega })`$, numerically calculated from (32), in $`\mathrm{\Omega }`$ fields with a size of $`2\stackrel{}{.}5\times 2\stackrel{}{.}5`$ (field of view of SDSS) after sorting them in decending order of pair density.

Now let the programme last for 5 years. Considering that the matter in question is detection of rather faint objects, observations have to be made on dark nights only at moon phases within $`\pm 5^d`$ from the new moon, which gives approximately $`36\%`$ of night time. Assuming the average observational night of $`9^h`$ the 5-year resource of time then will make $`n_{obs}=657`$ nights or $`T_{obs}=5913`$ hours. Thus, as every field is observed for $`t`$ hours or $`n`$ nights (and $`\mathrm{\Omega }=T_{obs}/t`$, $`n_{obs}/n`$) the expected number of detected pairs is

$$N_{exp}=D\left(\frac{T_{obs}}{t}\right)F(t),$$
(33)
or

$$N_{exp}=D\left(\frac{n_{obs}}{n}\right)F(n).(33^{})$$

These quantities derived from Figs. 7 and 10 are presented in Fig. 11. One can see from this figure that the optimal total exposure per field for WD-WD pairs is $`67`$ nights, and it is about 1 night when searching for WD-NS or WD-BH systems.

The expected numbers of detections with optimal strategies along with the optimal total exposures per field are given at Table 8.

A new project entitled ”The Dark Matter Telescope” (http://www.dmtelescope.org/index.htm) being developed involved 8.4 (effectively 6.9) metre telescope with a field of view slightly wider than that of the SDSS. The use of it in a manner similar to that described above increases the number of detectable objects by a factor of 2.2. Other wide field telescopes of 4 metres in diameter being built now, LAMOST (http://www.greenwich-observatory.co.uk/lamost.
html) and VISTA (http://www.vista.ac.uk), could detect 1.5 times more binaries than the SDSS.

Thus, for searching for the flares caused by gravitation lensing in binary systems with compact companions, the facilities of the SDSS can be used, having reduced the time of an individual exposure to 10-20 s. It is obvious that after recording a flare in a certain object, another telescope of similar class must be involved in observations. This object should be monitored to prove the effect and for detailed investigation of the detected system.

Thus, within the frame of the proposed programme, it will be possible to find several dozen binary systems comprised of compact objects. It will be recalled that in numerous surveys of microlensing, the number of reliably detected events amounts to about 350 (Alcock 1999), only in two cases the matter in question being of massive lenses (black holes?) (Bennett et al. 2001). When applying the technique suggested, there is a possibility of detecting with absolute assurance more than 10 black holes with open event horizons. An investigation of these objects with high time resolution within the MANIA experiment (Beskin et al. 1997) may permit at last the detection of observational evidence of extreme gravitational fields.

###### Acknowledgements.

This investigation was supported by the Russian Ministry of Science, Russian Foundation of Basic Researches (grant 01-02-17857), Federal Programs ”Astronomy” and ”Integration” and Science-Education Center ”Cosmion”. The authors are very grateful to the anonymous referee for very useful observations and E. Agol for important notes. AVT would like to thank D.A.Smirnov and S.V. Karpov for valuable discussions and his especially grateful to I.V. Chilingaryan for his help with various questions. We thank T.I. Tupolova for manuscript preparation.


format: text

text:
# Computational Methods for UV-Suppressed Fermions

## 1. Introduction

There has been recent interest in the so-called ‘Ultraviolet slowing down’ of fermionic simulations in lattice QCD \[Irving et al, 1998, Duncan, Eichten and Thacker, 1999, de Forcrand, 1999, Peardon, Lattice2001, A. Hasenfratz & Knechtli, 2001\]. These studies try to address algortmically large fluctuations of the high end modes of the fermion determinant. The goal is to increase the signal-to-noise ratio of the infrared modes and to accelerate fermion simulations as well.

In fact, all the computational effort needed to treat UV-modes by above algorithms can be reduced to zero by suppressing them in the first place \[Boriçi, 2002\]. The lattice Dirac operator of this fermion theory is given by:

(1.1)
$$D=\frac{\mu }{a}\mathrm{\Gamma }_5\mathrm{tanh}\frac{a\mathrm{\Gamma }_5D_{W/S/O}}{\mu }$$

where $`D_{W/s/o}`$ is the input lattice Dirac operator, $`a`$ the lattice spacing and $`\mu >0`$ is a dimensionless parameter. For Wilson (W) and overlap (o) fermions as the input theory one has $`\mathrm{\Gamma }_5=\gamma _5`$. For staggered fermions $`\mathrm{\Gamma }_5`$ is a diagonal matrix with entries $`+1/1`$ on even/odd lattice sites. The theory converges to the input theory in the contimuum limit and is local and unitary as shown in detail in \[Boriçi, 2002\]. The input theory is also recovered in the limit $`\mu \mathrm{}`$. For $`\mu 0`$ one has $`D\mu `$, i.e. a quenched theory.

Perturbative calculations with this theory are straightforward. To fix the idea I assume in the following Wilson fermions to be input fermions. The inverse fermion propagator is given by

(1.2)
$$\stackrel{~}{D}(p)=\frac{\mu }{a}\gamma _5\mathrm{tanh}\frac{a\stackrel{~}{H}_W(p)}{\mu }$$

with $`p=\{p_\nu ,\nu =1,\mathrm{},4\}`$ being the four-momentum vector. As usual, gauge fields are parametrized by $`su(3)`$ elements:

(1.3)
$$U(x)_\nu =e^{iagA(x)_\nu },A(x)_\nu su(3)$$

and the Wilson operator is written as a sum of the free and interaction terms:

$$D_W=D_W^0+D_W^I$$

The splitting of the lattice Dirac operator is written in the same form:

$$D=D^0+D^I,D^0=\frac{\mu }{a}\gamma _5\mathrm{tanh}\frac{aH_W^0}{\mu }$$

where the interaction term has to be determined. This can be done by expanding $`D`$ in terms of $`a/\mu `$:

(1.4)
$$D=D_W[\text{1l}+c_1(\frac{aH_W}{\mu })^2+c_2(\frac{aH_W}{\mu })^4+\mathrm{}]$$

where $`c_1,c_2,\mathrm{}`$ are real expansion coefficients. Calculation of $`D^I`$ is an easy task if one stays with a finite number of terms in the right hand side of (1.4). Also, the number of terms can be minimized using a Chebyshev approximation for the hyperbolic tangent. 111I would like to thank Joachim Hein for discussions related to lattice perturbation theory.

In this paper I describe computational methods needed to evaluate the effective action of the theory defined above. In particular, the complexity of the proposed Lanczos method does not depend on the input sparse matrix that describes a fermion theory on the lattice.

In the following section I derive a class of Lanczos based methods for computations with the proposed theory and then in section 3 conclusions follow.

## 2. Lanczos based methods for computations with fermions

The effective action of the theory defined above can be written as:

(2.1)
$$S_{\text{eff}}=\text{tr}f(A)$$

where $`A^{N\times N}`$ and $`f(s)`$ is a real and smooth function of $`s^+`$. The matrix $`A`$ is assumed to be Hermitian and positive definite. Since the trace is difficult to obtain one can use the stochastic method of \[Bai et al, 1996\]. The method is based on evaluations of many bilinear forms of the type:

(2.2)
$$(b,A)=b^Tf(A)b$$

where $`b^N`$ is a random vector. The trace is estimated as an average over many bilinear forms. A confidence interval can be computed as described in detail in \[Bai et al, 1996\].

The method described here is similar to the method of \[Bai et al, 1996\]. Its viability for lattice QCD computations has been demonstrated in the recent work of \[Cahill et al, 1999\]. \[Bai et al, 1996\] derive their method using quadrature rules and Lanczos polynomials. Here, I give an alternative derivation which uses familiar tools (in lattice simulations) such as sparse matrix invertions and Padé approximations. The Lanczos method enters the derivation as an algorithm for solving linear systems of the form:

(2.3)
$$Ax=b,x^N$$

### 2.1. Lanczos algorithm

I follow standard standard texts as \[Golub & Van Loan, 1989\] and notations and arguments of \[Boriçi, 1999b, Boriçi, 2000a, Boriçi, 2000b\]. $`n`$ steps of the Lanczos algorithm \[Lanczos, 1952\] on the pair $`(A,b)`$ are given by Algorithm 1.

The Lanczos vectors $`q_1,\mathrm{},q_n^N`$ can be compactly denoted by the matrix $`Q_n=[q_1,\mathrm{},q_n]`$. They are a basis of the Krylov subspace $`𝒦_n=\text{span}\{b,Ab,\mathrm{},A^{n1}b\}`$. It can be shown that the following identity holds:

(2.4)
$$AQ_n=Q_nT_n+\beta _nq_{n+1}e_n^T,q_1=b/b_2$$

$`e_n`$ is the last column of the identity matrix $`\text{1l}_n^{n\times n}`$ and $`T_n`$ is the tridiagonal and symmetric Lanczos matrix (2.5) given by:

(2.5)
$$T_n=\left(\begin{array}{cccc}\alpha _1& \beta _1& & \\ \beta _1& \alpha _2& \mathrm{}& \\ & \mathrm{}& \mathrm{}& \beta _{n1}\\ & & \beta _{n1}& \alpha _n\end{array}\right)$$

The matrix (2.5) is usually referred to as the Lanczos matrix. Its eigenvalues, the so called Ritz values, tend to approximate the extreme eigenvalues of the original matrix $`A`$.

To solve the linear system (2.3) I seek an approximate solution $`x_n𝒦_n`$ as a linear combination of the Lanczos vectors:

(2.6)
$$x_n=Q_ny_n,y_n^n$$

and project the linear system (2.3) on to the Krylov subspace $`𝒦_n`$:

$$Q_n^{}AQ_ny_n=Q_n^{}b=Q_n^{}q_1b_2$$

Using (2.4) and the orthonormality of Lanczos vectors, I obtain:

$$T_ny_n=e_1b_2$$

where $`e_1`$ is the first column of the identity matrix $`\text{1l}_n`$. By substituting $`y_n`$ into (2.6) one obtains the approximate solution:

(2.7)
$$x_n=Q_nT_n^1e_1b_2$$

### 2.2. Algorithms for the bilinear form (2.2)

The theoretical framework of the algorithm of \[Bai et al, 1996\] can be based on the Padé approximant of the smooth and bounded function $`f(.)`$ in an interval. The Padé approximation can be expressed as a partial fraction expansion. Therefore, one can write:

(2.8)
$$f(s)\underset{k=1}{\overset{m}{}}\frac{c_k}{s+d_k}$$

with $`c_k,d_k0,k=1,\mathrm{},m`$. It is assumed that the right hand side converges to the left hand side when the number of partial fractions becomes large enough. For the bilinear form I obtain:

(2.9)
$$(b,A)\underset{k=1}{\overset{m}{}}b^T\frac{c_k}{A+d_k\text{1l}}b$$

A first algorithm can already be written down at this point. Having computed the partial fraction coefficients one can use a multi-shift iterative solver of \[Freund, 1993\] to evaluate the right hand side (2.9). To see how this works, I solve the shifted linear system:

$$(A+d_k\text{1l})x^k=b$$

using the same Krylov subspace $`𝒦_n`$. A closer inspection of the Lanczos algorithm, Algorithm 1 suggests that in the presence of the shift $`d_k`$ I get:

$$\alpha _i^k=\alpha _i+d_k$$

while the rest of the algorithm remains the same. This is the so called shift-invariance of the Lanczos algorithm. From this property and by repeating the same arguments which led to (2.7) I get:

(2.10)
$$x_n^k=Q_n\frac{1}{T_n+d_k\text{1l}_n}e_1b_2$$

Using the shift-invariance of the Lanczos algorithm I obtain Algorithm 2.

Note that the residual errors $`r_i^k,i=1,\mathrm{},n,k=1,\mathrm{},m`$ are given by:

$$r_i^k=bAx_i^kd_kx_i^k$$

In exact arithmetic their norm is given by:

(2.11)
$$1/\rho _i^k=bAx_i^kd_kx_i^k_2$$

By applying Algorithm 2 one can solve the shifted linear systems on the right hand side of (2.9). The algorithm stops if the linear system with the smallest shift is solved to the desired accuracy $`ϵ`$. This is a well-known technique \[Freund, 1993\] which is used also in lattice QCD \[Frommer et al, 1995\]. However, the problem with this method is that one needs to store a large number of vectors that is proportional to $`m`$. This could be prohibitive if $`m`$ is say larger than $`10`$.

In fact, the right hand side of (2.9) can be written in terms of solutions $`x_n^k,k=1,\mathrm{},m`$ as a sum of scalars:

(2.12)
$$(b,A)\underset{k=1}{\overset{m}{}}c_kw^k,w^k=b^Tx^k$$

Therefore, it is easy to replace the vector recurrences by scalar recurrences of the form:

(2.13)
$$\stackrel{~}{w}_{i+1}^k=(\stackrel{~}{w}_i^k\alpha _i+\stackrel{~}{w}_{i1}^k\beta _{i1})/\beta _i$$

In this way one obtains the Algorithm 3. It is clear that by applying Algorithm 3 one gains substantial storage savings compared to Algorithm 2. If one has a good Padé approximant for the function $`f(.)`$ one can apply Algorithm 3.

Note that one another way to save storage is using the multi-shift Gonjugate Gradient variant of \[B.Jegerlehner, 1998\].

### 2.3. An exact method

If a Padé approximation is not sufficient or difficult to obtain, the Lanczos method is the only alternative to evaluate exactly the bilinear forms of type (2.2).

To see how this is realized I assume that the linear system (2.3) is solved to the desired accuracy using the Lanczos algorithm, Algorithm 1 and (2.7). In the application considered here one can show that:

(2.14)
$$\underset{k=1}{\overset{m}{}}b^T\frac{c_k}{A+d_k\text{1l}}b=b^2\underset{k=1}{\overset{m}{}}e_1^T\frac{c_k}{T_n+d_k\text{1l}_n}e_1$$

For the result (2.14) to hold, it is sufficient to show that:

$$b^T\frac{c_k}{A+d_k\text{1l}}b=b^2e_1^T\frac{c_k}{T_n+d_k\text{1l}_n}e_1$$

which can be shown using the orthonormality property of the Lanczos vectors and (2.10). Note however that in presence of roundoff errors the orthogonality of the Lanczos vectors is lost but the result (2.14) is still valid. The interested reader may consult the work of \[Cahill et al, 1999, Golub & Strakos, 1994\].

From this result and the convergence of the partial fractions to the matrix function $`f(.)`$, it is clear that:

(2.15)
$$(b,A)\widehat{}_n(b,A)=b^2e_1^Tf(T_n)e_1$$

Note that the evaluation of the right hand side is a much easier task than the evaluation of the right hand side of (2.2). A straightforward method is the spectral decomposition of the symmetric and tridiagonal matrix $`T_n`$:

(2.16)
$$T_n=Z_n\mathrm{\Omega }_nZ_n^T$$

where $`\mathrm{\Omega }_n^{n\times n}`$ is a diagonal matrix of eigenvalues $`\omega _1,\mathrm{},\omega _n`$ of $`T_n`$ and $`Z_n^{n\times n}`$ is the corrsponding matrix of eigenvectors, i.e. $`Z_n=[z_1,\mathrm{},z_n]`$. From (2.15) and (2.16) it is easy to show that (see for example \[Golub & Van Loan, 1989\]):

(2.17)
$$\widehat{}_n(b,A)=b^2e_1^TZ_nf(\mathrm{\Omega }_n)Z_n^Te_1$$

where the function $`f(.)`$ is now evaluated at individual eigenvalues of the tridiagonal matrix $`T_n`$.

The eigenvalues and eigenvectors of a symmetric and tridiagonal matrix can be computed by the QL method with implicit shifts \[Press et al, 1993\]. The method has an $`O(n^3)`$ complexity. Fortunately, one can compute (2.17) with only an $`O(n^2)`$ complexity. Closer inspection of eq. (2.17) shows that besides the eigenvalues, only the first elements of the eigenvectors are needed:

(2.18)
$$\widehat{}_n(b,A)=b^2\underset{i=1}{\overset{n}{}}z_{1i}^2f(\omega _i)$$

It is easy to see that the QL method delivers the eigenvalues and first elements of the eigenvectors with $`O(n^2)`$ complexity. 222I thank Alan Irving for the related comment on the QL implementation in \[Press et al, 1993\].

A similar formula (2.18) is suggested by \[Bai et al, 1996\]) based on quadrature rules and Lanczos polynomials. The Algorithm 4 is thus another way to compute the bilinear forms of the type (2.2).

Clearly, the Lanczos algorithm and Algorithm 3 has an $`O(nN)`$ complexity, whereas Algorithm 4 has a greater complexity: $`O(nN)+O(n^2)`$. However, Algorithm 4 delivers an exact evaluation of (2.2). For typical applications in lattice QCD the $`O(n/N)`$ overhead is small and therefore Algorithm 4 is the recommended algorithm among all three algorithms presented in this section.

A remark on stopping criteria is also desirable. The method of \[Bai et al, 1996\] computes the relative differences of (2.18) between two successive Lanczos steps and stops if they don’t decrease below a given accuracy. In order to perform the test their algorithm needs to compute the eigenvalues of $`T_i`$ at each Lanczos step $`i`$. This may be a large computational overhead. On the other hand the test proposed here is theoretically safe. This is clarified by the remark at the end of the proof of the lemma (2.14). However, this test may be too prudent since the prime interest here is the computation of the bilinear form (2.2).

To illustrate this situation I give an example from a lattice calculation (The lattice size and parameters are given in section 5.1). I compute the bilinear form (2.2) for:

(2.19)
$$f(s)=\mathrm{log}\mathrm{tanh}\sqrt{s},s^+$$

and $`A=H_W^2`$, $`b^N`$. The real and imaginary parts of the $`b`$elements are chosen randomly from the set $`\{+1,1\}`$.

In Fig. 1 are shown the normalized recursive residuals $`\rho _0/\rho _i=bAx_i_2/b_2,i=1,\mathrm{},n`$ and relative differences of (2.18) between two successive Lanczos steps. The figure illustrates clearly the different regimes of convergence for the linear system and the bilinear form. The relative differences of the bilinear form converge faster than the computed recursive residual. This example indicates that a stopping criterion based on the solution of the linear system may indeed be strong. Therefore, the recommended stopping criteria would be to monitor the relative differences of the bilinear form but less frequently than proposed by \[Bai et al, 1996\]. More investigations are needed to settle this issue. Note also the roundoff effects (see Fig. 1) in the convergence of the bilinear form.

## 3. Conclusion

In this paper I have described computational methods needed to evaluate the effective action of the theory with suppressed cutoff modes.

All methods described in this paper have a complexity that does not depend on the input sparse matrix that describes the fermion theory on the lattice. In this way, it may be concluded that simulation algorithms of lattice QCD which are based on the estimation of the effective fermion action have the same complexity. The UV-suppressed fermions are such an example.

## Acknowledgements

I would like to thank Philippe de Forcrand and Alan Irving for discussions and useful suggestions at different stages of this work.