Signatures of rotating binaries in micro-lensing experiments
Gravitational micro-lensing offers a powerful method to probe a variety of binary lens systems as the binarity of the lens introduces in the event light curves deviations from the typical (single lens) Paczyński behaviour. Generally, a static binary lens is considered to fit the observed light curve and, when the orbital motion is taken into account, an over-simplified model is usually employed. In this paper, we treat in a realistic way the binary lens motion and focus on simulated events well fitted by a Paczyński curve. We show that, most often, an accurate timing analysis of the residuals (calculated with respect to the best fit Paczyński model) is sufficient to infer the orbital period of the binary lens. It goes without saying that the independently estimated period may be used to further constrain the orbital parameters obtained by the best fitting procedure that often gives degenerate solutions. We also present a preliminary analysis on the event OGLE-2011-BLG-1127 / MOA-2011-BLG-322 which was recognized to be due to a binary lens. The period analysis results in a periodicity of days which confirms the oscillation of the observed data around the best fit model. The estimated periodicity is likely associated to an intrinsic variability of the source star, thus opening the possibility to use this technique to investigate either the intrinsic variability of the source and the effects induced by the binary lens orbital motion.
keywords:gravitational lensing: micro
In the last decades the gravitational micro-lensing technique, initially developed to search for Massive Compact halo objects (MACHOs) in the Galactic halo and in the Galactic disk (Paczyński, 1986; Alcock et al., 1993; Paczyński, 1996; Roulet & Mollerach, 1997, 2002; Zakharov & Sazhin, 1998), has been widely used to infer the presence of exo-planets111At the time of writing, the number of exo-planets detected via the micro-lensing method is 24 (see e.g. the extra-solar planet encyclopedia available at http://exoplanet.eu/). orbiting around the main lensing stars (see e.g. the reviews by Perryman 2000; Perryman et al. 2005; Bennett 2009). Indeed, as shown by Mao & Paczyński (1991) the presence of a planet around its hosting star forms a binary lens that can micro-lens a background source, inducing non-negligible deviations with respect to the usual symmetric Paczyński light curve (Witt & Mao 1994; Gould 1994; Alcock et al. 1997). In particular, the light-curve analysis of highly magnified events is sensitive to the presence of lens companions when the binary components are separated by a distance of the order of Einstein radius associated to the whole lens system. In general, such signatures are characterized by short duration deviations lasting from a few hours to a few days (depending on the parameters of the binary lens system). As a matter of fact, the micro-lensing technique is so sensitive that it allows to detect exo-planets in a rather large range of masses, spanning from Jupiter-like planets down to few Earth-mass objects (Bennett & Rhie, 1996; Bond et al., 2004; Udalski et al., 2005; Beaulieu et al., 2006; Gould et al., 2006; Gaudi et al., 2008; Bennett, 2009; Dominik, 2010; Gaudi, 2011).
Micro-lensing has also become a powerful tool to study several aspects of stellar astrophysics. In fact, a sufficiently amplified micro-lensing event gives the opportunity to investigate the source limb-darkening profile, i.e. the variation of intensity from the center of the disc up to its border, thus implying the possibility to test among different atmosphere models (Abe et al. 2003; Gaudi & Gould 1999). This opens the unique possibility to investigate the Galactic bulge star atmospheres, otherwise hardly to be studied due to their distance. Besides the brightness profile of the background source disc, highly magnified events with large radii sources allow to measure the lens Einstein radius once the physical radius of the source star is known from different methods. More recently, it has also been proposed (see, e.g., Ingrosso et al. 2012 and references therein) that, when finite size source effects become relevant, a characteristic polarization signal may arise in micro-lensing events. This is due to the differential magnification induced during the crossing of the source star over the lens (either single or binary). This produces a polarization signal up to 0.04 for late type stars and up to a few percent for cool giants, depending on the physical processes at the basis of the polarization and the atmosphere parameters of the source star. Such a signal may be observable with the currently available technology. It was also shown (Ingrosso et al. 2013) that polarization measurements may help in disentangling the particular degeneracy between the binary or planet lens solution which occurs in some micro-lensing events.
Note that finite source effects usually work against the appearance of exo-planet signatures in micro-lensing light curves, as these features tend to be smeared out. However, although the finite source effects are generally tiny, they cannot be ignored when modeling a light curve. This is particularly true for large planetary masses or in case of binary micro-lensing events characterized by caustic crossing, i.e. when the source passes through fold and/or cusp caustics (Schneider, Ehlers & Falco, 1992).
The main reason for considering binary events is that the best fit procedure to a micro-lensing light curve may allow, in principle, to derive the parameters (the projected lens separation and the mass ratio ) of the lens system. Typically, a static binary lens is considered and this implies the minimization of a functional depending on seven free parameters. These are, in addition to and , the time of closest approach to the lens system, the impact parameter (in units of the Einstein radius), the Einstein time of the event , the angle that the background source trajectory forms with respect to the binary lens separation, and the source radius . Considering in the fit procedure the orbital motion of the lens system (see e.g. Dominik 1997; Penny et al. 2011 a, 2011 b) is extremely time consuming since a good modeling of such motion would imply six additional parameters, i.e. the semi-major axis , the orbit eccentricity , the time of passage at periastron , the angle between the normal to lens orbital plane and the line of sight, the orientation of the orbital plane in the sky, and the orbiting versus (either clockwise or counterclockwise). As discussed in Park et al. (2013), for the determination of the lens parameters one should also include the relative lens-source parallax due to Earth motion around the Sun which involves two more parameters. As a result of the large number of parameters involved222Note that in the most general case one should also consider the baseline magnitude and the blending parameter., when a best fit procedure is attempted, some tricks are required in order to make the fit converging to reliable results.
In this paper, we will concentrate on binary systems with orbital parameters for which the resulting micro-lensing light curve is very close to a Paczyński curve (i.e. a planetary case). Indeed, one expects that the presence of planets rotating around the hosting star should cause, most often, only small perturbations (see also Bozza 1999) to the Paczyński light curve associated to an equivalent single lens event. We consider two cases: 1) the binary system orbital period is lower than the typical Einstein time of the event, and 2) is comparable or larger than . In the latter, we required that the light curve is long enough to contain at least three full cycles. Then we showed that an accurate timing analysis of the residuals (calculated with respect to the Paczyński model) is sufficient to infer .
The plan of the paper is as follows: in Section 2, we describe the model adopted for simulating a micro-lensing light curve taking into account the orbital motion of the binary system and the finite source size effects. In Section 3, for selected simulated events, we apply the timing analysis to best fit residuals and, finally, in Section 4, we discuss the advantages of our procedure and present a preliminary analysis on the OGLE-2011-BLG-1127 / MOA-2011-BLG-322 event, and address future perspectives.
2 Simulating a rotating binary micro-lensing event
2.1 The complex notation for a binary with point-like masses
It is well known (see e.g. Schneider, Ehlers & Falco 1992) that in a binary micro-lensing event involving a point-like source the angular positions of the images () and binary lens components ( and ) in the lens plane are related to the angular source position () in the source plane via the lens equation
where and are the masses of the binary lens components normalized to the total system mass with , and . All the positions are measured in units of the angular Einstein radius of the whole lens system, i.e.
being and (with ) the distance from the observer to the source and lens, respectively. Here, is the gravitational constant and the speed of light. Eq. (1) can be easily generalized in order to account for multiple lenses as in Kayser, Refsdal & Stabell (1986). By using complex notation (see Witt 1990; Witt & Petters 1993; Witt & Mao 1995), eq. (1) can be easily rewritten in terms of complex variables. Let us define the variables and , being and the dimensionless components of the vector and , respectively. With this substitution one has
where the bar over a symbol indicates the operator of complex conjugation.
To determine the image position and magnification, one has to take the complex conjugate of eq. (3) and substitute the expression for back in it. After some algebra, one gets a fifth-order polynomial in , i. e. with coefficients depending on , and , whose solutions give the image positions. Here, represents the separation of the two binary lens components in a reference frame () with origin in the middle point and axis oriented from the star to the planet and so that and , i.e. the star and its companion are always on the (real) axis. The vertical axis is oriented in order to have a right-handed reference frame. A (fixed) reference frame , with origin in the center of mass and axes and parallel to and , will be always used in the following.
The lensing phenomenon separates the source star image into several pieces () and conserves the source brightness (see, e.g., Schneider, Ehlers & Falco 1992). Hence, for a point-like source, the summation over the magnification of each of the images gives the total magnification (Witt & Mao, 1995)
where is the parity of the i-th image and the determinant of the Jacobian is
The point-source magnification diverges if any of the images appear on the critical curve of the lens, i.e. the curve along which vanishes. This curve, once mapped back in the source plane via the lens equation in eq. (1) or, equivalently, in eq. (3), gives the caustic curves. As shown by Witt (1990), the critical curves can be obtained by solving
with real phase varying from to . The previous equation corresponds to a fourth-order complex polynomial in z, the roots of which333In the present work we use the new and robust algorithm to solve polynomial equations presented in Skowron & Gould (2012) and optimized for the micro-lensing case. The code is a few times faster than similar methods designed for more general purposes. lie on the critical curves. By varying and searching for the roots of eq. (6), one finally gets the caustic curves via the lens equation444For a binary system of point-like masses, the caustic curve is constituted either by a single continuous curve or two (or three) separate paths depending on the binary parameters and (Schneider, Ehlers & Falco 1992). For a complete picture of the caustic zoo, we refer to Pejcha & Heyrovský (2009)..
2.2 The Inverse Ray Shooting (IRS) method
The complex method described above fails when the source is very close to the lens caustics. In such conditions, a more accurate (but more time consuming) approach is given by the IRS technique (see e.g. Schneider, Ehlers & Falco (1992); Kayser, Refsdal & Stabell (1986); Wambsganss (1997)). Note that the IRS method is a general technique, thus allowing one to consider all those effects which can not be easily implemented when generating a synthetic light curve. The method consists in solving numerically the lens equation in eq. (3) by shooting light rays backwards from the observer through the lens plane up to the source plane. The rays which are collected in the source plane, i.e. those for which the two sides of the lens equation differ at least by a chosen tolerance, and the density of the rays at a particular location in the source plane is proportional to the magnification at this point. Also in this case the magnification map depends on the mass ratio and the projected binary lens separation . As one can easily verify, the complex and IRS methods give comparable results beyond a distance (in Einstein radii in the source plane) from the fold and cusp caustics. Hence, we use the IRS method for any source distance less than , and the complex method otherwise. We verified that a good compromise between the calculation speed and the quality of the output simulated light curve was obtained by choosing .
2.3 Finite source effects
For real astrophysical source stars, the point-like source approximation described above breaks down where the magnification varies non linearly on the scales comparable to the source radius . In this situation, the magnification is obtained by integrating in eq. (4) - weighted with the surface brightness - over the source disk and dividing by the unamplified source flux. In particular,
where is the center of the background source star on its trajectory at a given time and the surface integral is extended over the whole surface .
We numerically solve555The numerical integration is performed by using the Cuhre method implemented in Hahn (2005). the previous integral when the source (or part of it) is close enough (or crosses) the caustic curves associated to the binary lens system, i.e. when the IRS method has to be employed ( ). Otherwise, for any source position for which the complex formalism (Witt, 1990) is a good approximation, we use the hexadecapole technique (Gould 2008). The hexadecapole method substantially consists in calculating the magnification in 13 locations of the source disk: one position coincident with the star center -the monopole approximation-, eight positions on the source limb -with radius -, and four positions on a ring of radius . As a last remark, we remind that we assumed the source surface brightness to be described by the linear limb-darkening profile (Choi et al., 2012)
where the parameter depends on wavelength, spectral type, surface gravity and metallicity of the source star. In the following, we assume as a typical value (Claret, 2000).
2.4 The lens system orbital motion
For a binary system constituted by two masses and orbiting around the common center of mass, the reduced mass moves on an ellipse, with semi-major axis and eccentricity . The orbit equation666In this work, we assume that the binary lens motion occurs with respect to a right-handed reference frame () with origin in one of the ellipse foci. The trajectory is in the plane, being the axis orthogonal to it. In general, the orbit of the reduced mass is also characterized by a semi-major axis forming a phase angle (measured counterclockwise) with the axis. Finally, the orbital plane is seen by the distant observer with an inclination angle between the axis and the line of sight. Note that three rotations and a traslation (depending on time ) are required in order to switch between the and reference frames. is (see e.g. Smart 1980)
where is the radial distance from the center of mass and the true anomaly, respectively. The true anomaly is related to the mean anomaly by the Kepler equation
being a generic instant of time during the orbital motion, the time of the pericenter passage and the angular frequency given, as usual, by
where is the keplerian orbital period. Since we are giving the semi-major axis in units of the Einstein radius, the orbital period results to be given in Einstein time.
We then solved the Kepler equation by using a series expansion of Bessel functions (see e.g. Watson 1966). In particular, it can be easily shown that
where is the Bessel function of order , and the mean anomaly contains the time dependence. Note that a few terms in the previous summation are usually required to get the solution of the Kepler equation within a few percent of accuracy. Hence, for any time , the previous equation gives the true anomaly with the reduced mass position obtained by means of eq. (9). With this formalism, the radial coordinates of the two lenses on their orbits are
Of course, the separation (which enters in the complex and IRS methods described previously) between the two lens components is now time dependent and it is given by the projection onto the lens plane of the vector .
3 The output synthetic light curve and the timing analysis on the residuals
The formalism introduced in the previous Section allows us to obtain simulated micro-lensing light curves taking into account the orbital motion of the binary lens. Here, we concentrate on binary systems characterized by orbital periods lower than the typical Einstein time of the event and by orbital parameters and giving rise to event light curves close to the Paczyński behaviour.
In Fig. 1 (left side), we give a snapshot at a generic time of the position of the source, the rotating lens system and the associated caustic curves (inset) in the source plane (upper panel). In the additional material available on-line, we present a movie of the micro-lensing event corresponding to Fig. 1. By measuring all the distances in units of the Einstein radius, we simulated a binary system rotating counterclockwise with parameters , , , and . The corresponding binary system has an orbital period , being the Einstein time of the event. The purple ellipse represents the planet trajectory as projected in the plane of the sky. The background source (with radius ) is moving at the impact parameter (as measured from the binary center of mass) at an angle with respect to the axis. With the assumed binary system and micro-lensing parameters, we expect to have only minor deviations (a planetary case) with respect to the single lens light curve. Indeed, in the middle panel, we show the magnification light curve (red line) resulting from eq. (7) to which the best fit light curve (green line) obtained with a Paczyński model is superimposed. As one can see, the two curves are almost indistinguishable in this case. As it is clear from the bottom panel of Fig. 1, the residuals between the simulated light curve and the best fit model are as low as a few parts per thousand but clearly show a periodic behaviour (here, for simplicity, we are ignoring the effect of any noise in the data strain). As a guide for the eye, we also show a series of vertical solid lines separated in time by one orbital period . In the right side of the figure, we simulated a micro-lensing event due to a binary system with the same mass and semi-major axis as before, but different values of eccentricity, inclination and phase angle (, , and ).
Before discussing the results obtained with the period analysis on the residuals, we note that, due to the fact that the simulated light curve and the best-fit Paczyński model do not have the peak at the same time (and this is a general characteristic of such events), the residuals show an asymmetry with respect to and a number of extra peaks (especially concentrated around the event maximum) that are not correlated to the real binary system period. As we will discuss later, this introduces spurious peaks in the period analysis that might make impossible to recover the simulated period. Indeed, any algorithm used for the period search looks for features in phase and fails around the micro-lensing maximum to give the right period.
We tested different methods to extract the period from the residual light curve for several simulated events: the Fast Fourier Transform -FFT- (Press et al., 2007); the Phase Dispersion Minimization -PDM- (Stellingwerf, 1978); the Lomb-Scargle periodogram (Lomb, 1976; Scargle, 1982); its generalized form given in Zechmeister & Kuerster (2009) and the epoch folding technique. All the algorithms give consistent results, but here we prefer to present those obtained via the classical Lomb-Scargle method (the generalized version did not improved the period detection) because its statistical behaviour is well known and the technique can be easily implemented.
The method requires to specify the minimum () and maximum () frequencies to be searched for in the input signal, together with the frequency step . We chose to set these parameters depending on the observational data set, i.e. , , oversampling by a factor 10, being the duration of the observation and the associated time step. Note that by using the minimum frequency we are implicitly requiring to have at least three full cycles per observational window.
The blind application of the Lomb-Scargle method to the residual light curves resulted in the periodograms shown in Fig. 1, where the dashed vertical lines (always at period ) mark the periodicity detected in the signal. Note that the Lomb-Scargle periodogram (as well as the other period search algorithms) does not give exactly the simulated period since the frequency of the periodic signal is affected by the relative source-lens motion. As discussed above, the central part of the residual light curve has spurious peaks that disturb the timing analysis since they introduce additional power at frequencies different from the fundamental one, i.e .
We have found that the best results in the period search are obtained by removing a central region in the residual curve, around the event peak. Without this cut, the Lomb-Scargle periodogram may return other peaks in addition to that corresponding to the simulated one. It is remarkable that, in all the cases we have considered, it is sufficient to remove a very small region around the maximum magnification with size of the order of a fraction of the orbital period in order to get the true periodicity (indicated as ).
The second case we discuss is a rotating binary similar to that considered in Penny et al. (2011 b) with orbital parameters , , , and rotating with a period . Since the system is face-on and the eccentricity is null, the projected distance does not change during the micro-lensing event, thus implying that the caustics simply rotate without any deformation. This particular geometry appears in the periodogram as a peak always at half of the simulated period (see the left bottom panel of Fig.2) as a consequence of the North-South symmetry in the caustic plane. In the right side of the same figure, we simulate a binary lens event with the same orbital parameters but with and an observational window of . In spite of the fact that there are only a few full cycles, the periodogram gives again half of the expected period.
We note that the geometry of the binary lens described in the two examples of Fig. 2 is rather unlikely to appear in real observations and so, in Fig. 3, we consider a more realistic binary configuration with and , being all the other parameters unchanged. As one can see, in this case the Lomb-Scargle analysis allows us to recover the right period even if the associated periodogram peak is quiet broad.
4 Result and discussions
We have considered the orbital binary lens motion in simulated micro-lensing events that have light-curves similar to the Paczyński behaviour and have shown that in most cases a timing analysis of the residuals allows one to recover the binary period . Since the simulated light curve and the best-fit Paczyński model do not have, in general, the peak at the same time, the residuals show an asymmetry with respect to and a number of extra peaks that are not correlated to the real binary system period. This introduces spurious features in the analysis that can be overcome by excluding a central region of the light curve whose size is reduced smoothly until the true period is lost.
The importance of this procedure is that of inferring the orbital period of the binary lens without performing a numerically heavy multi-dimensional fit procedure which, generally, do not give unique solutions. Here, by considering Paczyński-like events, we showed that the procedure is robust enough to recover the simulated orbital period, provided that the event duration contains at least three full cycles.
A systematic search for such signatures on real (archival) observed micro-lensing events will be presented elsewhere. Here, we show a preliminary analysis on the event OGLE-2011-BLG-1127 / MOA-2011-BLG-322 (Shvartzvald et al. 2013) possibly due to a binary lens. By using the publicly available OGLE data (Udalski et al. 2003, http://ogle.astrouw.edu.pl/) we have fitted the data with a simple Paczyński model (see Fig. 4, panel a) and analyzed the residuals (Fig. 4, panel b) by using either the generalized Lomb-Scargle and the discrete fast Fourier transform techniques (DFFT treated via the CLEAN algorithm of Roberts et al 1987). These are the best-performing period search methods for this kind of data. Both techniques return the presence of a period at about days (see Fig. 4, panels c and d)777Here we note that the region around the event peak has not been removed since already poor of data..
It seems that this periodic feature is remarkably stable since it always appears when a time resolved periodogram is performed, i.e. when the period search algorithm is applied on different parts of the residual light curve, provided that each part contains a sufficiently large number of cycles. It is therefore possible that the periodicity at days is associated to an intrinsic variability of the source star as if it is in a binary system or intrinsically variable (see e.g. Wyrzykowski et al. 2006) or contaminated by a close variable star. Note also that some period search techniques, such as the cleaned DFFT (Fig. 4, panels c) and the generalized Lomb-Scargle periodogram, give significant power also at days and days. Since we test trial periods in a given range, it is possible to evaluate the significance of each feature compared to the power at all other frequencies. Hence, following Horne & Baliunas (1986) (see also Zechmeister & Kuerster 2009), we evaluate the false alarm probability and obtained the significance levels at , , and which are shown at the bottom-right panel of Fig. 4 as solid, dotted and dashed lines, respectively. As one can note, the periods at days and days are not significant.
Note also that interpreting the days feature as a true periodicity could be problematic since it is comparable to the event Einstein time days. In addition, a time resolved period analysis does not return a stable result. As a matter of fact, we note that the light curve of the event OGLE-2011-BLG-1127 shows an oscillation around the best fit (static) model, as it is apparent in Fig. 1 of Shvartzvald et al. (2013), with a time scale of days. If the bump-like structure observed days after the event peak is really present along the whole light curve (as also mentioned in Shvartzvald et al. 2013) then our analysis is returning a periodicity related to the source or to the binary lens motion.
- Abe et al. (2003) Abe, F., et al., 2003, A&A, 411, L493.
- Alcock et al. (1993) Alcock, C., et al., 1993, Nature, 365, 621.
- Alcock et al. (1997) Alcock, C., et al., 1997, ApJ, 491, 436.
- Beaulieu et al. (2006) Beaulieu, J.-P., et al. 2006, Nature, 439, 437.
- Bennett (2009) Bennett, D. P., 2009, Exoplanets. Mason J., editor. Springer, Berlin; 2009.
- Bennett & Rhie (1996) Bennett, D. P. & Rhie, S. H., 1996, ApJ, 472, 660.
- Bond et al. (2004) Bond, I. A., et al., 2004, ApJ, 606, L155.
- Bozza (1999) Bozza, V., 1999, ApJ, 348, 311.
- Claret (2000) Claret, A., 2000, A&A, 363, 1081.
- Choi et al. (2012) Choi, J.-Y., Shin, I.-G., Park, S.-Y. et al., 2012, ApJ 751, 41.
- Dominik (1997) Dominik, M., 1997, Galactic Microlensing Beyond The Standard Model, Ph.D. Thesis at University of Dortmund.
- Dominik (2010) Dominik, M., 2010, Gen. Rel. Grav., 42, 2075.
- Gaudi & Gould (1999) Gaudi, B. S. & Gould, A., 1999, ApJ, 513 619.
- Gaudi et al. (2008) Gaudi B. S., et al., 2008, Science, 319, 927.
- Gaudi (2011) Gaudi, B. S., Exoplanets. Seager S., editor. Tucson: Univ. Arizona Press; 2011. p. 79.
- Gould (1994) Gould, A., 1994, ApJ, 421, L71.
- Gould et al. (2006) Gould, A., et al., 2006, ApJ, 644, L37.
- Gould (2008) Gould, A., 2008, ApJ, 681, 1593.
- Hahn (2005) Hahn, T., 2005, Computer Physics Communications, 168, 2, 78.
- Horne & Baliunas (1986) Horne, J.H., & Baliunas, S. L., 1986, ApJ, 302, 757.
- Kayser, Refsdal & Stabell (1986) Kayser, R., Refsdal, S. & Stabell, R., 1986, A&A, 166, 36.
- Ingrosso et al. (2012) Ingrosso, G., et al., 2012, MNRAS, 426, 1496.
- Ingrosso et al. (2013) Ingrosso, G., et al., 2013, Physica Scripta. 2013. preprint arXiv:1310.5866.
- Lomb (1976) Lomb, N.R., 1976, Ap&SS, 39, 447.
- Mao & Paczyński (1991) Mao, S., & Paczyński, B., 1991, ApJ, 374, L37.
- Paczyński (1986) Paczyński, B., 1986, ApJ, 304, 1.
- Paczyński (1996) Paczyński, B., 1996, ARA&A, 34, 419.
- Park et al. (2013) Park, H., et al., 2013, eprint arXiv:1306.3744.
- Pejcha & Heyrovský (2009) Pejcha, O., & Heyrovský, D., 2009, ApJ, 690, 1772.
- Penny et al. (2011 a) Penny, M.T., Mao, S., & Kerins, E., 2011, MNRAS, 412, 607.
- Penny et al. (2011 b) Penny, M.T., Kerins, E., & Mao, S., 2011, MNRAS, 417, 2216.
- Perryman (2000) Perryman, M., 2000, Rep. Prog. Phys., 63, 1209.
- Perryman et al. (2005) Perryman, M., et al., 2005, Report by the ESA–ESO Working group on Extra-Solar Planets, astro-ph/0506163.
- Press et al. (2007) Press, W. H.. Teukolsky, S.A., Vetterling, W.T. & Flannery, B.P., Numerical Recipes. The Art of Scientific Computing, 3rd edn. Cambridge: Cambridge University Press, 2007.
- Roberts et al (1987) Roberts, D. H., Lehar, J., & Dreher, J. W., 1987, AJ, 93, 968.
- Roulet & Mollerach (1997) Roulet, E. & Mollerach, S., 1997, Phys. Rep., 279, 2.
- Roulet & Mollerach (2002) Roulet, E. & Mollerach, S., Gravitational Lensing and Microlensing. Singapore: World Scientific, 2002.
- Scargle (1982) Scargle, J.D., 1982, ApJ, 263, 835.
- Schneider, Ehlers & Falco (1992) Schneider, P., Ehlers, J. & Falco E. E., , Gravitational Lenses. Berlin: Springer, 1992.
- Shvartzvald et al. (2013) Shvartzvald, Y., et al., 2013, preprint (arXiv:1310.0008).
- Skowron & Gould (2012) Skowron, J., & Gould, A., 2012, preprint (astro-ph/1203.1034).
- Smart (1980) Smart, W., M., Textbook on Spherical Astronomy, Cambridge: Cambridge University Press, 1980.
- Stellingwerf (1978) Stellingwerf, R.F., 1978, ApJ, 224, 953.
- Udalski et al. (2003) Udalski, A., et al., 2003, Acta Astronomica, 53, 4.
- Udalski et al. (2005) Udalski, A., et al., 2005, ApJ 628, L109.
- Wambsganss (1997) Wambsganss, J., 1997, MNRAS, 284, 172.
- Watson (1966) Watson, G.N., A Treatise on the Theory of Bessel Functions, Camdridge: Cambridge University Press, 1966.
- Witt (1990) Witt, H. J.,1990, A&A, 236, 311.
- Witt & Mao (1994) Witt, H. J., & Mao, S., 1994, ApJ, 430, 505.
- Witt & Mao (1995) Witt, H. J., & Mao, S., 1995, ApJ, 447, L105.
- Witt & Petters (1993) Witt, H. J., & Petters, A.O., 1993, Journal of Mathematical Phsyics, 34, 4093.
- Wyrzykowski et al. (2006) Wyrzykowski, L., et al., 2006, Acta Astronomica, 56, 145.
- Zakharov & Sazhin (1998) Zakharov A.F. & Sazhin, M.V., 1998, Soviet Physics Uspekhi, 41, 945.
- Zechmeister & Kuerster (2009) Zechmeister, M. & Kuerster, M., 2009, A&A, 496, 577.