EDP Sciences logo
Subscriber Authentication Point
Search
Menu
Home All issues Volume 515 (June 2010) A&A, 515 (2010) A52 Full HTML
Free Access
Issue
A&A
Volume 515, June 2010
Article Number A52
Number of page(s) 7
Section Planets and planetary systems
DOI https://doi.org/10.1051/0004-6361/200913755
Published online 08 June 2010
A&A 515, A52 (2010)

Bayesian analysis of caustic-crossing microlensing events

A. Cassan1,2 - K. Horne3 - N. Kains3 - Y. Tsapras4,5 - P. Browne3

1 - Institut d'Astrophysique de Paris, UMR 7095 CNRS, 98 bis boulevard Arago, 75014 Paris, France
2 - Université Pierre & Marie Curie, Paris, France
3 - Scottish Universities Physics Alliance, School of Physics & Astronomy, University of St Andrews, North Haugh, KY169SS, UK
4 - Las Cumbres Observatory, 6740B Cortona Dr, suite 102, Goleta, CA 93117, USA
5 - School of Mathematical Sciences, Queen Mary University of London, Mile End Road, London E1 4NS, UK

Received 27 November 2009 / Accepted 23 February 2010

Abstract
Aims. Caustic-crossing binary-lens microlensing events are important anomalous events because they are capable of detecting an extrasolar planet companion orbiting the lens star. Fast and robust modelling methods are thus of prime interest in helping to decide whether a planet is detected by an event. Cassan introduced a new set of parameters to model binary-lens events, which are closely related to properties of the light curve. In this work, we explain how Bayesian priors can be added to this framework, and investigate on interesting options.
Methods. We develop a mathematical formulation that allows us to compute analytically the priors on the new parameters, given some previous knowledge about other physical quantities. We explicitly compute the priors for a number of interesting cases, and show how this can be implemented in a fully Bayesian, Markov chain Monte Carlo algorithm.
Results. Using Bayesian priors can accelerate microlens fitting codes by reducing the time spent considering physically implausible models, and helps us to discriminate between alternative models based on the physical plausibility of their parameters.

Key words: gravitational lensing: micro - methods: analytical - methods: statistical - planetary systems

1 Introduction

Mao & Paczynski (1991) first suggested that observations of Galactic gravitational microlensing events could lead to the discovery of extrasolar planets. Microlensing involves the time-dependent brightening and then dimming of a background source star as an intervening massive object (the lens) crosses the observer line-of-sight. Light rays from the source bend in the vicinity of the lens, focusing them toward the observer. Since 1994, survey teams such as OGLE[*] (OGLE III, Udalski 2003) and MOA[*] (Bond et al. 2001) have reported more than four thousand microlensing events toward the Galactic bulge to date. Several hundreds of these events have been carefully selected and densely sampled by follow-up networks such as PLANET[*], $\mu$FUN[*], RoboNet[*], and MiNDSTEp[*]. Although microlensing teams have so far published only nine exoplanet detections, the method itself stands out because of its high sensitivity to low-mass planets with orbits of several astronomical units. It thus probes in the planet mass-separation plane a region beyond reach of any other technique, as demonstrated by the detection of the very first cool super-Earth, OGLE 2005-BLG-390Lb (Kubas et al. 2008; Beaulieu et al. 2006).

A number of microlensing events exhibit anomalous behaviour (i.e., they cannot be adequately modelled by the standard single-lens light curve, e.g., Paczynski 1986) and some of these anomalies can be attributed to lensing by binary objects. The types of light curves produced by binary lensing form a rich tapestry but, in general, binary systems with two equal mass components tend to exhibit pronounced, long anomalies in their light curves, whereas when the secondary companion is only a small fraction of the total mass, the anomalies can be quite short and subtle. It is primarily these latter types of anomalies that may be caused by star-planet binaries (Gould & Loeb 1992; Mao & Paczynski 1991). Nevertheless, because the true nature of the anomaly cannot always be established while the microlensing event is still ongoing, every binary-lens microlensing event constitutes a prime target for planet hunting.

In binary lensing, the lens system configuration delineates regions of space on the source plane that are bound by gravitational caustics. Caustics are closed curves with concave segments that meet in outward pointing cusps, defined by the location where the Jacobian determinant of the lens mapping equation vanishes, i.e., are lines of infinite point-source magnification. There are three kinds of caustic topologies, which depend on the values of the binary lens mass ratio q and the two component projected separation d in angular Einstein ring radius $\theta_{\rm E}$ (Einstein 1936)

\begin{displaymath}\theta_{\rm E}= \sqrt{\frac{4GM}{c^2} \left(\frac{D_{\rm S}-D_{\rm L}}{D_{\rm S}~D_{\rm L}}\right)}, \end{displaymath} (1)

where $D_{\rm S}$, $D_{\rm L}$ are the observer-source and observer-lens distances and M the lens total mass. In the close separation regime (cf. Fig. 1 of Cassan 2008), there are three caustics, one central (4-cusp) and two (3-cusp) planetary caustics. In the intermediate regime, there is only one (6-cusp) caustic, and in the wide separation regime, there is one central and one planetary caustic (both with 4 cusps).

In many cases, the source trajectory happens to cross a caustic. As the source crosses the caustic curve and enters the enclosed area, a new pair of images appears, causing a sudden increase in the observed brightness. In a similar way, when the source exits the area defined by the caustics, the two images merge and disappear, causing a rapid drop in the observed brightness. These dramatic changes in magnification result in readily recognisable jumps in microlensing light curves. As emphasised by Cassan (2008), the ingress and egress times $t_{\rm in}$ and $t_{\rm out}$ may be restricted to within very tight intervals when caustic crossing features have been identified in the light curve, and thus advantageously used as alternative modelling parameters.

The new set of binary-lens modelling parameters introduced by Cassan (2008) have the advantage that two of these parameters are very closely related to features that can be directly identified in the light curve. Using this new formulation to analyse the data of OGLE 2007-BLG-472 in its most straightforward implementation as a maximum likelihood analysis (``minimising $\chi^2$''), Kains et al. (2009) unveiled a subtle aspect of binary-lens modelling: relatively improbable physical models with very large values of ${t_{\rm E}}$ were found with $\chi^2$ values lower than other more plausible models. To avoid finding parameter combinations that are physically unlikely, dramatic progress can be achieved by switching to a Bayesian analysis. This is desirable as the Bayesian approach makes use of prior information on the underlying physical parameters, while $\chi^2$ says nothing about parameter plausibility.

In this article, we show how to derive Bayesian priors for the caustic-crossing binary-lens parameters defined by Cassan (2008). These are based on physical priors on quantities that can be estimated from Galactic models or calculated from already observed events (Sects. 2 and 3). In Sect. 4, we describe an implementation of this Bayesian formalism within a Markov chain Monte Carlo fitting scheme, using in particular priors on the Einstein time ${t_{\rm E}}$ (time for the source to travel an angular distance $\theta_{\rm E}$).

2 Maximum likelihood versus Bayesian fitting

Cassan (2008) introduced a new parameterisation of the binary lens microlens light curve model that is well suited to describing caustic-crossing events. In this formalism, the caustic curve in the source plane is parameterised by a curvilinear abscissa (or arc length) from 0 to 2. The trajectory of a source crossing a caustic, which is classically parameterised by its impact parameter ${u_{\rm0}}$ and position angle $\alpha$, can alternatively be defined by giving the values $s_{\rm in}$ at ingress and $s_{\rm out}$ at egress[*]. The two parameters timing the trajectory, ${t_{\rm E}}$ (time to cross one Einstein radius) and ${t_{\rm0}}$ (date at minimum impact parameter ${u_{\rm0}}$), are then replaced by the ingress and egress times $t_{\rm in}$ and $t_{\rm out}$. The caustic curve is specified in the source (i.e., caustic) plane by a complex function $\zeta(s)=\xi(s)+{\rm i}\eta(s)$ (see Sect. 3.2), and once $s_{\rm in}$ and $s_{\rm out}$ are specified, the source trajectory is fully defined. This bijective switch of parameters, $({u_{\rm0}},\alpha,{t_{\rm E}},{t_{\rm0}}) \mapsto (s_{\rm in},s_{\rm out},t_{\rm in},t_{\rm out})$, takes advantage of the relatively high precision with which $t_{\rm in}$ and $t_{\rm out}$ can be inferred from the observations (Kubas et al. 2005; Kains et al. 2009).

Using these new parameters, Kains et al. (2009) analysed the caustic crossing event OGLE 2007-BLG-472. The approach taken was a maximum likelihood procedure, quantifying the ``goodness-of-fit'' by a $\chi^2$ statistic, and minimising the $\chi^2$ to optimise the fit. A grid search in (d,q) with even spacing in $\log{d}$ and $\log{q}$ was conducted. For each (d,q) caustic configuration, a genetic algorithm was used to explore widely the remaining parameter space. While $(s_{\rm in},s_{\rm out})$ covered the full range of possibilities, $[0,2]\times[0,2]$, $t_{\rm in}$ and $t_{\rm out}$ evolved in very tight intervals based on the values inferred from the light curve features (caustic crossing magnification peaks). These first fits were refined using a Markov chain Monte Carlo (MCMC) algorithm, again holding (d,q) fixed while optimising the remaining parameters. The best-fit models in each of the identified best-fit regions were then found by allowing all parameters to vary.

As expected for binary lens events, the resulting $\chi^2(d,q)$ maps uncovered a variety of widely-separated model parameter regions where a relatively low $\chi^2$ could be achieved. The lowest $\chi^2$ models corresponded to very low q, in the planet-mass regime. But with a short duration between the caustic entry and exit, and a planetary caustic size scaling as q1/2, these models implied an extremely long Einstein time ${t_{\rm E}}\sim (t_{\rm out}-t_{\rm in})/q^{1/2} \sim 10^4$ days, which is very unlikely according to kinematics of stars motions within the Milky Way. These best-fit maximum likelihood models were therefore rejected on this physical argument. This need to reject the lowest $\chi^2$ models highlights a weakness in the maximum likelihood approach, which neglects prior distributions on the parameter space. On the other hand, Bayesian parameter estimation takes proper account of prior distributions in the parameter space (see e.g., Trotta 2008, for a review of astrophysical applications).

In a Bayesian analysis, the posterior probability distribution over the model parameters $\theta$ is a function of the data D

\begin{displaymath}P(\theta\vert D) = \frac{P(D\vert\theta) P(\theta)} {\int P(D\vert\theta) P(\theta)~\d\theta}, \end{displaymath} (2)

where $P(\theta)$ is the prior probability distribution on the parameters, and the denominator partition function ensures proper normalisation of the posterior as a probability distribution over the parameters $\theta$. The likelihood $P(D\vert\theta)$ is a function of the parameters $\theta$ and a probability distribution over the data D. For Gaussian measurement errors with N data points having individual standard deviations $\sigma_i$, the likelihood is

\begin{displaymath}L(\theta) \equiv P(D\vert\theta) = \frac{ \exp{ \left\{ - \f... ...} } } {\left(2\pi\right)^{N/2}\prod_{i=1}^{N}\sigma_i} \cdot \end{displaymath} (3)

Since maximising the likelihood corresponds to minimising

\begin{displaymath}-2\ln{L(\theta)} = \chi^2 + 2 \sum_{i=1}^{N} \ln{ \sigma_i } + \frac{N}{2} \ln{ 2\pi } ~, \end{displaymath} (4)

a maximum likelihood is equivalent to a minimum in $\chi^2$ when the error bars $\sigma_i$ are known, and is essentially a Bayesian analysis that implicitly assumes a prior that is uniform on the chosen parameters intervals. As we show in the next section, assuming more realistic priors could substantially affect the fitting process.

3 A Bayesian prior for (s$_{\rm in}$,s $_{\rm out}$)

3.1 Distribution of (s$_{\rm in}$,s $_{\rm out}$) for isotropic trajectories

A uniform prior probability distribution in the parameter square $(s_{\rm in},s_{\rm out})$ is implicit in the maximum likelihood analysis. Because of the non-linear correspondence between the two sets of parameters, it should correspond to a rather unlikely prior for the $({u_{\rm0}},\alpha)$ source trajectory parameters. A more plausible prior would for example arrange for the source trajectories to be uniformly distributed and isotropic in orientation.

In Fig. 1, the top panel shows an intermediate caustic with d=1.1 and q=0,1 (i.e., six cusps, in orange) with several crossing trajectories. It can be seen that a straight line may cross the caustic at two (black line), four (red line), or six (blue line) locations, depending on the number and orientation of the cusps. In the bottom panel, $\sim $104 of these trajectories were randomly shot and their corresponding position in the $(s_{\rm in},s_{\rm out})$ square reported, using the same colour convention. Trajectories with a single pair of ingress and egress map into unique black points, while for red and blue trajectory lines, there are respectively two and three possible pairs of ingress and egress points.

\begin{figure} \par\hspace*{5mm}\includegraphics[width=7.4cm]{13755fig1a.eps}\\ \includegraphics[width=9cm]{13755fig1b.eps} \end{figure} Figure 1:

The top panel illustrates the three kinds of possible source trajectories crossing a caustic: the black line has a single pair of ingress and egress points, while the red and blue lines have, respectively, two and three ingress/egress points. The bottom panel shows in the $(s_{\rm in},s_{\rm out})$ square the locations of a random distribution of $\sim $104 of these trajectories crossing a (d=1.1, q=0,1)-caustic, with the same colour convention. The solid vertical and horizontal black line mark the s-locations of the caustic cusps.

Open with DEXTER

We can understand some of the structures in the bottom panel of Fig. 1 as follows: vertical and horizontal lines marking the s values at the cusps divide the $(s_{\rm in},s_{\rm out})$ square into boxes. No trajectories appear in the boxes along the diagonal because the caustics curve concavely outward. It is thus impossible for a line that enters at some position between two cusps to exit at any point between those two cusps. In a similar way, other empty regions correspond to ingress/egress pairs that cannot be realised by straight lines crossing the caustic. The trajectories are seen to bunch up around ingress/egress pairs occurring close to a cusp. This happens because any trajectory entering close to a cusp is very likely (for a wide range of angles) to also exit near the same cusp.

3.2 Analytical formulation

We develop a mathematical formulation that allows us to compute analytically priors on $(s_{\rm in},s_{\rm out})$. The lens equation for a binary lens with separation d and mass ratio q defines the mapping of the position of a point-source $\zeta$ on the source plane to the positions of its three or five images at z on the lens plane

\begin{displaymath}\zeta = z - \frac{1}{1+q}~\left(\frac{1}{\overline{z}} + \frac{q}{\overline{z}+d}\right), \end{displaymath} (5)

where the more massive body is at the centre and the companion on the left-hand side. Following Witt (1990), the caustic lines $\zeta$ are parametrised by a parameter $\phi \in [0,2\pi]$

\begin{displaymath}\frac{1}{1+q}\left[\frac{1}{z^2} + \frac{q}{(z+d)^2}\right] = {\rm e}^{-{\rm i}\phi} , \end{displaymath} (6)

where z and $\zeta$ satisfy Eq. (5). For a given angle $\phi$, it is possible to solve a fourth order polynomial equation in z to obtain the corresponding caustic points $\zeta$. While the parameter $\phi$ is used here to write the useful formulae, in practice we use instead the equivalent parameter $s = s(\phi)$ (bijection) introduced by Cassan (2008), which has the advantage of sampling the caustics evenly.
\begin{figure} \par\mbox{\includegraphics[width=3.5cm]{13755fig2a.eps} \includeg... ...13755fig2d.eps} \includegraphics[width=3.5cm]{13755fig2e.eps} } \end{figure} Figure 2:

Bayesian prior $P{(s_{\rm in},s_{\rm out})}$ as a function of $s_{\rm in}$ (horizontal axis) and $s_{\rm out}$ (vertical axis), where we have assumed isotropic source trajectories and uniform distributions for ${t_{\rm E}}$ and event rate. Higher values of P appear in white (linear scale). From left to right, the caustic configurations are: a) intermediate with d=1.1 and q=0.1; b) wide+central and c) wide+secondary caustic, both for d=2 and q=0.1; d) close+central and e) close+secondary caustic, both for d=0.5 and q=0.1.

Open with DEXTER

To write more condensed formulae, we use notations that resemble two-dimensional vector operations. Given two complex numbers $\zeta_1 = \xi_1 +{\rm i} \eta_1$ and $\zeta_2 = \xi_2 +{\rm i} \eta_2$, we write $\zeta_1\wedge\zeta_2 = \xi_1 \eta_2 - \eta_1 \xi_2$ (``wedge product'') and $\zeta_1\cdot\zeta_2 = \xi_1 \xi_2 + \eta_1 \eta_2$ (``scalar product''), which are both real numbers. Moreover, a quantity related to a caustic entry (exit) is indicated by a subscript ``${\rm in}$'' (``${\rm out}$''). Using the usual convention that ${u_{\rm0}}> 0$ when the origin of the coordinate system stays on the right-hand side of the source trajectory, one can write

\begin{displaymath}{u_{\rm0}}= \frac{\zeta_{\rm out}\wedge\zeta_{\rm in}}{\vert\zeta_{\rm out}-\zeta_{\rm in}\vert} , \end{displaymath} (7)

\begin{displaymath}\alpha = \arctan\left(\frac{\eta_{\rm out}-\eta_{\rm in}}{\x... ...}}\right) + \pi ~ H\left(\xi_{\rm in}-\xi_{\rm out}\right) , \end{displaymath} (8)

\begin{displaymath}{t_{\rm E}}= \frac{t_{\rm out}-t_{\rm in}}{\vert\zeta_{\rm out}-\zeta_{\rm in}\vert} , \end{displaymath} (9)

\begin{displaymath}{t_{\rm0}}= \frac{t_{\rm out}+t_{\rm in}}{2} - (t_{\rm out}-t... ...eft\vert\zeta_{\rm out}-\zeta_{\rm in}\right\vert^2}\right], \end{displaymath} (10)

where H is the Heaviside step function, and $\alpha = \frac{\pi}{2} ~ {\rm sign}\left(\eta_{\rm out}-\eta_{\rm in}\right)$ if $\xi_{\rm out}=\xi_{\rm in}$. The transformation between the two sets of parameters is given by the Jacobian

\begin{displaymath}J = \left\vert \frac{\partial\left({u_{\rm0}},\alpha,{t_{\... ..._{\rm out},t_{\rm in},t_{\rm out}\right)} \right\vert \cdot \end{displaymath} (11)

Since the dependencies of the classical parameters with respect to the new ones are ${u_{\rm0}}(s_{\rm in},s_{\rm out})$, $\alpha(s_{\rm in},s_{\rm out})$, ${t_{\rm E}}(s_{\rm in},s_{\rm out},t_{\rm in},t_{\rm out})$, and ${t_{\rm0}}(s_{\rm in},s_{\rm out},t_{\rm in},t_{\rm out})$, J reads (using $\phi$ instead of s)

\begin{displaymath}J = \left\vert \begin{array}{cccc} \frac{\partial{u_{\rm0}... ...\partial\left(t_{\rm in},t_{\rm out}\right)}\right\vert \cdot \end{displaymath} (12)

After some algebra, we find for the components of the two latter Jacobians
                                $\displaystyle \frac{\partial{u_{\rm0}}}{\partial\phi_{\rm in}}$ = $\displaystyle \frac{\partial{u_{\rm0}}}{\partial\xi_{\rm in}}~\frac{{\rm d}\xi_... ...rm0}}}{\partial\eta_{\rm in}}~\frac{{\rm d}\eta_{\rm in}}{{\rm d}\phi_{\rm in}}$  
  = $\displaystyle \frac{\left(\zeta_{\rm out}-\zeta_{\rm in}\right)\cdot\zeta_{\rm ... ...m in}\right) \wedge \frac{{\rm d}\zeta_{\rm in}}{{\rm d}\phi_{\rm in}}\right] ,$ (13)
$\displaystyle \frac{\partial{u_{\rm0}}}{\partial\phi_{\rm out}}$ = $\displaystyle - \frac{\left(\zeta_{\rm out}-\zeta_{\rm in}\right)\cdot\zeta_{\r... ...in}\right) \wedge \frac{{\rm d}\zeta_{\rm out}}{{\rm d}\phi_{\rm out}}\right] ,$ (14)
$\displaystyle \frac{\partial\alpha}{\partial\phi_{\rm in}}$ = $\displaystyle \frac{\partial\alpha}{\partial\xi_{\rm in}}~\frac{{\rm d}\xi_{\rm... ...alpha}{\partial\eta_{\rm in}}~\frac{{\rm d}\eta_{\rm in}}{{\rm d}\phi_{\rm in}}$  
  = $\displaystyle - \frac{\left(\zeta_{\rm out}-\zeta_{\rm in}\right) \wedge \frac{... ...{\rm d}\phi_{\rm in}}}{\left\vert\zeta_{\rm out}-\zeta_{\rm in}\right\vert^2} ,$ (15)
$\displaystyle \frac{\partial\alpha}{\partial\phi_{\rm out}}$ = $\displaystyle \frac{\left(\zeta_{\rm out}-\zeta_{\rm in}\right) \wedge \frac{{\... ...\rm d}\phi_{\rm out}}}{\left\vert\zeta_{\rm out}-\zeta_{\rm in}\right\vert^2} ,$ (16)
$\displaystyle \frac{\partial{t_{\rm E}}}{\partial t_{\rm in}}$ = $\displaystyle -\frac{1}{\left\vert\zeta_{\rm out}-\zeta_{\rm in}\right\vert} ,$ (17)
$\displaystyle \frac{\partial{t_{\rm E}}}{\partial t_{\rm out}}$ = $\displaystyle \frac{1}{\left\vert\zeta_{\rm out}-\zeta_{\rm in}\right\vert} ,$ (18)
$\displaystyle \frac{\partial{t_{\rm0}}}{\partial t_{\rm in}}$ = $\displaystyle \frac{1}{2} + \left[\frac{1}{2}\frac{\zeta_{\rm out}+\zeta_{\rm i... ...zeta_{\rm in}} {\left\vert\zeta_{\rm out}-\zeta_{\rm in}\right\vert^2}\right] ,$ (19)
$\displaystyle \frac{\partial{t_{\rm0}}}{\partial t_{\rm out}}$ = $\displaystyle \frac{1}{2} - \left[\frac{1}{2}\frac{\zeta_{\rm out}+\zeta_{\rm i... ...zeta_{\rm in}} {\left\vert\zeta_{\rm out}-\zeta_{\rm in}\right\vert^2}\right] ,$ (20)

so that
                         $\displaystyle \left\vert\frac{\partial\left({u_{\rm0}},\alpha\right)} {\partial\left(\phi_{\rm in},\phi_{\rm out}\right)}\right\vert$ = $\displaystyle \frac{\left\vert\left(\zeta_{\rm out}-\zeta_{\rm in}\right) \wedg... ...rm out}}\right\vert} {\left\vert\zeta_{\rm out}-\zeta_{\rm in}\right\vert^3} ~,$ (21)
$\displaystyle \left\vert\frac{\partial\left({t_{\rm E}},{t_{\rm0}}\right)} {\partial\left(t_{\rm in},t_{\rm out}\right)}\right\vert$ = $\displaystyle \frac{\frac{\partial{t_{\rm0}}}{\partial t_{\rm in}} + \frac{\par... ...}\right\vert} = \frac{1}{\left\vert\zeta_{\rm out}-\zeta_{\rm in}\right\vert} ,$ (22)

which gives

\begin{displaymath}J = \frac{\left\vert\left(\zeta_{\rm out}-\zeta_{\rm in}\ri... ... {\left\vert\zeta_{\rm out}-\zeta_{\rm in}\right\vert^4}\cdot \end{displaymath} (23)

The derivatives ${\rm d}\zeta/{\rm d}\phi$ evaluated at the caustic entry and exit are given by (Cassan 2008)

\begin{displaymath}\frac{{\rm d}\zeta}{{\rm d}\phi} = \frac{{\rm d} z}{{\rm d}\p... ...\rm e}^{{\rm i}\phi}~\overline\frac{{\rm d} z}{{\rm d}\phi} , \end{displaymath} (24)

where

\begin{displaymath}\frac{{\rm d}{z}}{{\rm d}\phi} = \frac{{\rm i}}{2}~\frac{({z}+d)^2+q~{z}^2} {({z}+d)^3+q~{z}^3}~({z}+d)~{z}~. \end{displaymath} (25)

In the limit of cusp-crossing trajectories, i.e., $\zeta_{\rm out}-\zeta_{\rm in}\rightarrow0$, J behaves like $1/\vert\zeta_{\rm out}-\zeta_{\rm in}\vert^2$.

As expected, the Jacobian J is a function of the two parameters $(s_{\rm in},s_{\rm out})$, while the bijection between the two set of parameters was possible by involving $(t_{\rm in},t_{\rm out})$. However, J is not yet the Bayesian prior $P{(s_{\rm in},s_{\rm out})}$ we seek. We have yet to consider two aspects. Firstly, the parameters $({u_{\rm0}},\alpha,{t_{\rm E}},{t_{\rm0}})$ are themselves affected by prior probability distributions; this is discussed at the end of this section and is the topic of Sect. 4.1. Secondly, caustic crossing points are either entries or exits since the trajectory is orientated from $t_{\rm in}$ to $t_{\rm out}$ ( $t_{\rm out}\geq t_{\rm in}$), which is not accounted for in Eq. (23). To solve this second issue, we calculate the outward normal vector to the caustics at point $\zeta$,

\begin{displaymath}N_c = {\rm i}~\frac{{\rm d}\zeta}{{\rm d}\phi} \left/ \: \left\vert\frac{{\rm d}\zeta}{{\rm d}\phi}\right\vert\right. \end{displaymath} (26)

as well as the normalised and orientated trajectory vector

\begin{displaymath}N_t =\frac{\zeta_{\rm out}-\zeta_{\rm in}}{\left\vert\zeta_{\rm out}-\zeta_{\rm in}\right\vert} ~, \end{displaymath} (27)

and check whether $N_{c,{\rm in}} \cdot N_{t,{\rm in}} < 0$ (inward motion at $\zeta_{\rm in}$) and $N_{c,{\rm out}} \cdot N_{t,{\rm out}} > 0$ (outward motion at $\zeta_{\rm out}$). If these conditions are fulfilled, we write $P{(s_{\rm in},s_{\rm out})}=J$, and 0 otherwise.

Defined in this way, $P{(s_{\rm in},s_{\rm out})}$ is thus the prior on $(s_{\rm in},s_{\rm out})$ that we seek, in the special case of isotropic source trajectories (uniform distributions for ${u_{\rm0}}$ and $\alpha$), uniform microlensing events rate ( ${t_{\rm0}}$ is a random number), and uniform Einstein time ${t_{\rm E}}\geq 0$. In Fig. 2, we have plotted $P{(s_{\rm in},s_{\rm out})}$ for various (d,q) configurations as a function of $s_{\rm in}$ (horizontal axis) and $s_{\rm out}$ (vertical axis), higher values of P appearing in white (linear scale). From left to right, these configurations are: (a) intermediate with d=1.1 and q=0.1; (b) wide+central and (c) wide+secondary, both configurations for d=2 and q=0.1; (d) close+central and (e) close+secondary caustic, both for d=0.5 and q=0.1. One can compare the intermediate case plot with Fig. 1. In Sect. 4.1, we investigate how assuming different priors on the Einstein time ${t_{\rm E}}$ affect the prior on $(s_{\rm in},s_{\rm out})$.

3.3 Extended sources

When the source approaches the caustic curves (at typically less than three projected source radii), one needs to take into account extended source effects in the modelling. As for $t_{\rm in}$ and $t_{\rm out}$, it is usually possible to extract from the light curve a new parameter that can be used instead of the source radius.

It is well known that when the source crosses a straight line caustic (which is in many cases a good approximation of a real caustic), one can easily infer the duration of the crossing from the shape of the caustic crossing feature itself (Cassan et al. 2004; Schneider & Wagoner 1987; Albrow et al. 1999). Here, we define this duration as the time for the source to cross the caustic line by its full radius (i.e., from centre to limb), so that $\Delta t_{\rm cc} = \rho_{\rm *}/v_{\perp}$. In this definition, $\rho_{\rm *}$ is the source radius in Einstein ring radius units, $v_{\perp}$ is the component of the source velocity perpendicular to the caustic, and the subscript ``cc'' refers to either the caustic entry (``in'') or exit (``out''). For a given absolute velocity $1/{t_{\rm E}}$, the source will take longer to cross the caustic if the trajectory makes a tangential angle with it. More precisely, the normal velocity is proportional to the cosine of the angle between the trajectory and the caustic normal $v_{\perp} = \vert N_{c,{\rm cc}} \cdot N_{t,{\rm cc}}\vert/{t_{\rm E}}$. Inserting into this equation the expressions for ${t_{\rm E}}$, $N_{c,{\rm cc}}$, and $N_{t,{\rm cc}}$ (Eqs. (9), (26), and (27), respectively), we can compute the source radius $\rho_{\rm *}$ as a function of $\Delta t_{\rm cc}$

\begin{displaymath}\rho_{\rm *}= \frac {\left\vert(\zeta_{\rm out}-\zeta_{\rm i... ...cc}}{{\rm d}\phi_{\rm cc}}\right\vert} ~\Delta t_{\rm cc} . \end{displaymath} (28)

This expression would be exact if the crossed caustic were a perfect and infinite straight line. In reality, however, caustic curves always have a curvature, and sometimes the source partly crosses a cusp. Nevertheless, there is no arguing that $\Delta t_{\rm cc}$ is more suitable than $\rho_{\rm *}$ for parameterising the observed caustic crossing, since its rough value can be estimated from the light curve, in contrast to $\rho_{\rm *}$. In practice, we choose the caustic crossing that provides the most comprehensive data coverage and which ressembles most closely a straight line caustic crossing to extract the starting value when fitting $\Delta t_{\rm cc}$.

4 Markov Chain Monte Carlo fitting

4.1 Examples of prior probability distributions

For a given set of fitting parameters $(s_{\rm in},s_{\rm out},t_{\rm in},t_{\rm out},\Delta t_{\rm cc})$, the prior of the probed model is given by

\begin{displaymath}P({\rm model})= P{(s_{\rm in},s_{\rm out})} ~ P(t_{\rm in}, t_{\rm out}, \Delta t_{\rm cc}) . \end{displaymath} (29)

The prior $P{(s_{\rm in},s_{\rm out})}$ is computed as explained in Sect. 3.2, and may include priors that have been defined using the other parameters ${u_{\rm0}}$, $\alpha$, ${t_{\rm E}}$, ${t_{\rm0}}$, or $\rho_{\rm *}$ by properly weighting $P{(s_{\rm in},s_{\rm out})}$. Given Eqs. (9) and (10), a prior $P(t_{\rm in}, t_{\rm out})$ is equivalent to a prior $P({t_{\rm0}}, {t_{\rm E}})$ with a corresponding change in the prior $P{(s_{\rm in},s_{\rm out})}$. We now discuss different priors for the various parameters that could realistically be used in the Bayesian analysis. In Fig. 2 for example, we illustrate the case of isotropic trajectories, which corresponds to uniform priors for the parameters ${u_{\rm0}}$ and $\alpha$. This choice is justifiable, since the direction of the binary lens axis is random.

\begin{figure} \par\includegraphics[width=5cm]{13755fig3.eps} \end{figure} Figure 3:

Prior $P{(s_{\rm in},s_{\rm out})}$ for the solution configuration of caustic crossing event OGLE 2002-BLG-069 (close+central caustic), with an underlying (uninformative) uniform prior in $\log{t_{\rm E}}$. The red cross shows the location of the found caustic crossing, (1.3,0.3), which falls in a region of relative high probability.

Open with DEXTER

The first class of priors that we can use are uninformative priors. Since the prior expresses information about the values of parameters before any data has been taken, we know that parameters such as ${t_{\rm0}}$, ${t_{\rm E}}$, or $\Delta t_{\rm cc}$ must have uninformative priors, because we can only estimate their values by examining the light curve. Although it is natural to use uniform priors for ${t_{\rm0}}$, $\alpha$ or ${u_{\rm0}}$, for strictly positive parameters such as $\Delta t_{\rm cc}$ or ${t_{\rm E}}$, it is more suitable and commonly decided to use an uninformative prior that is uniform in the logarithm of the parameter.

We illustrate the use of an uninformative prior (uniform priors in $\log{t_{\rm E}}$, in ${u_{\rm0}}$, $\alpha$, and ${t_{\rm0}}$) by computing $P{(s_{\rm in},s_{\rm out})}$ for the solution configuration of the binary lens event OGLE 2002-BLG-069 (Kubas et al. 2005; Cassan et al. 2004). The configuration for that event was that of a source crossing the central caustic of a close binary lens with parameters d=0.46, q=0.58 and $t_{\rm out}-t_{\rm in}\simeq 14.5$ days. The resulting prior $P{(s_{\rm in},s_{\rm out})}$ is plotted in Fig. 3, where the red cross shows the location of the caustic crossings at $s_{\rm in}\simeq 1.3$, $s_{\rm out}\simeq 0.3$. This falls within a region of high probability, meaning that the corresponding $P{(s_{\rm in},s_{\rm out})}$ prior would have been a reasonable choice for this event.

The second class of priors are those that we can derive using information known before the event is observed. In microlensing, a convenient parameter on which such a prior can be placed is the Einstein time ${t_{\rm E}}$. This parameter depends on the relative distances between the source, the lens, and the observer, the kinematics of both the lens and the source and the lens' mass function. Combining all these data can help us to determine which ranges of values of ${t_{\rm E}}$ are more likely to be observed. For the event OGLE-2007-BLG-472 (Kains et al. 2009), no prior information was included on ${t_{\rm E}}$ (or the prior was assumed to be uninformative), which cause the best-fit models to have unrealistically long ${t_{\rm E}}$.

The method presented here can indeed be extended to include informative priors on parameters other than ${t_{\rm E}}$, such as the source flux distribution, the blending light due to the lens, the relative proper motion of the source and lens, or the source-radius caustic crossing-time, but this would require us to link the analysis to a Monte Carlo model of the Galaxy. Although our approach can be generalised to these possible extensions, they are beyond the scope of the present paper. Using ${t_{\rm E}}$ also has the advantage that its statistical distribution is fairly well-constrained by observed single-lens light curves, since this parameter is common to single- and binary-lens events.

\begin{figure} \par\includegraphics[width=9cm, bb= 33 0 760 443]{13755fig4a.eps}... ...e*{2cm}\includegraphics[width=5cm]{13755fig4b.eps}\hspace*{2cm}} \end{figure} Figure 4:

In the top panel, the histogram (blue rectangles) shows the distribution of ${t_{\rm E}}$ found after fitting 788 single-lens microlensing events from OGLE 2006-2007 seasons. The solid black line shows the model prediction of Wood & Mao (2005), which is in good agreement with the data. The bottom panel displays the prior $P{(s_{\rm in},s_{\rm out})}$ for the same intermediate caustic configuration as for Fig. 2 (d=1.1, q=0.1), but assuming an underlying prior for ${t_{\rm E}}$ given by the above distribution.

Open with DEXTER

Empirical distributions of ${t_{\rm E}}$ can be obtained by modelling a large number of observed microlensing events. The top panel of Fig. 4 shows a histogram (blue rectangles) of ${t_{\rm E}}$ values found by fitting 788 single-lens microlensing events from the 2006-2007 OGLE seasons (including blending). As expected, the distribution is far from uniform but instead appears roughly log-normal with a peak close to $\log{t_{\rm E}}\simeq 1.32$ and $\sigma_{\log{t_{\rm E}}} \simeq 0.4$. Theoretical distributions of ${t_{\rm E}}$ can also be based on predictions obtained with a Galactic model, such as the distribution advocated by Wood & Mao (2005). This is plotted as a solid black line on top of our histogram (Fig. 4, top panel) and is seen to closely match the empirical distribution. Nevertheless, the distribution of Wood & Mao (2005) lacks both extremely long events (say ${t_{\rm E}}> 300$ days) that can be interpreted as black hole lenses, and extremely short events (say ${t_{\rm E}}< 3$ days) that can be interpreted as evidence of a population of free floating planets. But selection effects cause these extreme events to be under-represented in the observed ${t_{\rm E}}$ distribution, as can be seen in Fig. 4. For these exceptional cases, special treatment would be required, for example using a prior on ${t_{\rm E}}$ that is more generous to extreme values in an attempt to compensate for selection effects. For most of binary lens events, however, a mild discrimination against black hole or loose planet lenses seems appropriate.

Using the Wood & Mao (2005) distribution as a prior, we compute and plot (Fig. 4, bottom panel) the corresponding distribution $P{(s_{\rm in},s_{\rm out})}$ by assuming $(t_{\rm out}-t_{\rm in}) = 20$ days, d=1.1, and q=0.1 (the same intermediate configuration as Fig. 2). Figure 4 (bottom panel) shows that with this prior, cusp-crossing trajectories are far less likely to happen. For a trajectory near the cusps, this is because the source has only a short distance to travel between the entry and exit, while $(t_{\rm out}-t_{\rm in})$ is constant, meaning that the source's motion has to be very slow, leading to large values of ${t_{\rm E}}$, which are now ruled out by the prior[*]. This effect can be seen directly in the plot of $P{(s_{\rm in},s_{\rm out})}$, where strong ``wing'' features at the cusps disappear, and other features appear (compare with Fig. 2).

4.2 Posterior probability distributions: MCMC fitting

In practice, these and other statistics related to the posterior parameter distribution can be evaluated efficiently using a Markov chain Monte Carlo to evaluate the probability-weighted integrals in Bayes' theorem. A random walk in the parameter space is undertaken by taking random steps drawn from a distribution of the parameters $\theta$. Each proposed step is accepted or rejected based on the probability of the new point relative to the old one exceeding some threshold, which is adjusted to maintain the acceptance rate above roughly 20-30%. The resulting chain locates and wanders around a local minimum, sampling the parameters with a weight proportional to the posterior probability.

For a maximum likelihood analysis, the relative probability used to accept or reject new steps is $\exp{\left\{-\Delta\chi^2/2\right\}}$ alone, where $\Delta\chi^2$ is the $\chi^2$ difference between the new and old points; in a full Bayesian analysis, we multiply this exponential factor by the ratio of new to old values of the prior $P({\rm model})$, following Eq. (29). The posterior probability that the parameters $\theta$ lie in a defined region $\Theta$ is then

\begin{displaymath}P(\theta \in \Theta) = \int_\Theta P(\theta\vert D)~ \d \theta . \end{displaymath} (30)

The expected value of any function of parameters, $g(\theta)$, is

\begin{displaymath}\left< g \right> \equiv \int g(\theta)~ P(\theta\vert D)~ \d \theta , \end{displaymath} (31)

and the variance about that expected value is

\begin{displaymath}{\rm Var} \left[ g(\theta) \right] \equiv \int \left( g(\th... ... - \left< g \right> \right)^2~ P(\theta\vert D)~ \d \theta . \end{displaymath} (32)

In a similar way, confidence intervals, parameter covariances, and confidence intervals can all be evaluated easily in the usual manner given the posterior probability distribution found with the MCMC algorithm, providing us with a complete statistical picture of the parameter space that we explore.


5 Conclusion

We have investigated plausible priors for Bayesian analysis of caustic-crossing microlensing light curves, based on an alternative parameterisation introduced by Cassan (2008). We have developed a mathematical formulation that allows us to compute analytically Bayesian priors for these parameters, given the knowledge we have about the physical quantities on which they depend. A number of relevant priors that may be used in a Bayesian, Markov chain Monte Carlo implementation of the given equations have been explored.

In the context of the rapid development of a new generation of networks of classical and robotic telescopes (e.g., Tsapras et al. 2009), as well as space-based observations such as with the ESA project satellite Euclid (Beaulieu et al. 2010), a current challenge facing the microlens planet search community is to fully automate the fitting of binary lens light curves in real time, after having detected an anomaly (e.g., Horne et al. 2009). This would enable anomalies that are detected in the observed light curves to be characterised as quickly as possible and for us to ascertain whether the anomalous behaviour is caused by a planet-mass companion of the lens star. Identifying parameters that could be estimated automatically by analysing the light curve (e.g., a magnification jump due to a caustic crossing) is already a step forward in accelerating the fitting codes by exploring a far more tighter parameter space. This was the motivation of Cassan (2008) in defining a new set of parameters. In this work, we have added the possibility of including Bayesian priors in the analysis, which would avoid the need to explore combinations of parameters that are unlikely to happen.


References

  1. Albrow, M. D., Beaulieu, J.-P., Caldwell, J. A. R., et al. 1999, ApJ, 522, 1022 [NASA ADS] [CrossRef] [Google Scholar]
  2. Beaulieu, J. P., Bennett, D. P., Batista, V., et al. 2010, unpublished[arXiv:1001.3349] [Google Scholar]
  3. Beaulieu, J.-P., Bennett, D. P., Fouqué, P., et al. 2006, Nature, 439, 437 [NASA ADS] [CrossRef] [PubMed] [Google Scholar]
  4. Bond, I. A., Abe, F., Dodd, R. J., et al. 2001, MNRAS, 327, 868 [NASA ADS] [CrossRef] [Google Scholar]
  5. Cassan, A. 2008, A&A, 491, 587 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  6. Cassan, A., Beaulieu, J. P., Brillant, S., et al. 2004, A&A, 419, L1 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  7. Einstein, A. 1936, Science, 84, 506 [NASA ADS] [CrossRef] [PubMed] [Google Scholar]
  8. Gould, A., & Loeb, A. 1992, ApJ, 396, 104 [NASA ADS] [CrossRef] [Google Scholar]
  9. Horne, K., Snodgrass, C., & Tsapras, Y. 2009, MNRAS, 396, 2087 [NASA ADS] [CrossRef] [Google Scholar]
  10. Kains, N., Cassan, A., Horne, K., et al. 2009, MNRAS, 395, 787 [NASA ADS] [CrossRef] [Google Scholar]
  11. Kubas, D., Cassan, A., Beaulieu, J. P., et al. 2005, A&A, 435, 941 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  12. Kubas, D., Cassan, A., Dominik, M., et al. 2008, A&A, 483, 317 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  13. Mao, S., & Paczynski, B. 1991, ApJ, 374, L37 [NASA ADS] [CrossRef] [Google Scholar]
  14. Paczynski, B. 1986, ApJ, 304, 1 [NASA ADS] [CrossRef] [Google Scholar]
  15. Schneider, P., & Wagoner, R. V. 1987, ApJ, 314, 154 [NASA ADS] [CrossRef] [Google Scholar]
  16. Trotta, R. 2008, Contemporary Physics, 49, 71 [NASA ADS] [CrossRef] [Google Scholar]
  17. Tsapras, Y., Street, R., Horne, K., et al. 2009, AN, 330, 4 [Google Scholar]
  18. Udalski, A. 2003, Acta Astron., 53, 291 [NASA ADS] [Google Scholar]
  19. Witt, H. J. 1990, A&A, 236, 311 [NASA ADS] [Google Scholar]
  20. Wood, A., & Mao, S. 2005, MNRAS, 362, 945 [NASA ADS] [CrossRef] [Google Scholar]

Footnotes

... OGLE[*]
http://www.astrouw.edu.pl/ ogle
... MOA[*]
http://www.phys.canterbury.ac.nz/moa
... PLANET[*]
http://planet.iap.fr
...$\mu$FUN[*]
http://www.astronomy.ohio-state.edu/ microfun
... RoboNet[*]
http://robonet.lcogt.net
... MiNDSTEp[*]
http://www.mindstep-science.org
... egress[*]
We use the notations ``in'' and ``out'' in place of ``entry'' and ``exit'' of Cassan (2008) to write more condensed formulae.
... prior[*]
More precisely, when ${t_{\rm E}}\rightarrow\infty$, Wood & Mao (2005) ${t_{\rm E}}$distribution behaves like $1/{t_{\rm E}}^3 \sim \vert\zeta_{\rm out}-\zeta_{\rm in}\vert^{3}$, and since $J\sim 1/\vert\zeta_{\rm out}-\zeta_{\rm in}\vert^2$, the net result is that near cusps, $J\sim\vert\zeta_{\rm out}-\zeta_{\rm in}\vert\rightarrow 0$.
Copyright ESO 2010

All Figures

  \begin{figure} \par\hspace*{5mm}\includegraphics[width=7.4cm]{13755fig1a.eps}\\ \includegraphics[width=9cm]{13755fig1b.eps} \end{figure} Figure 1:

The top panel illustrates the three kinds of possible source trajectories crossing a caustic: the black line has a single pair of ingress and egress points, while the red and blue lines have, respectively, two and three ingress/egress points. The bottom panel shows in the $(s_{\rm in},s_{\rm out})$ square the locations of a random distribution of $\sim $104 of these trajectories crossing a (d=1.1, q=0,1)-caustic, with the same colour convention. The solid vertical and horizontal black line mark the s-locations of the caustic cusps.

Open with DEXTER
In the text

  \begin{figure} \par\mbox{\includegraphics[width=3.5cm]{13755fig2a.eps} \includeg... ...13755fig2d.eps} \includegraphics[width=3.5cm]{13755fig2e.eps} } \end{figure} Figure 2:

Bayesian prior $P{(s_{\rm in},s_{\rm out})}$ as a function of $s_{\rm in}$ (horizontal axis) and $s_{\rm out}$ (vertical axis), where we have assumed isotropic source trajectories and uniform distributions for ${t_{\rm E}}$ and event rate. Higher values of P appear in white (linear scale). From left to right, the caustic configurations are: a) intermediate with d=1.1 and q=0.1; b) wide+central and c) wide+secondary caustic, both for d=2 and q=0.1; d) close+central and e) close+secondary caustic, both for d=0.5 and q=0.1.

Open with DEXTER
In the text

  \begin{figure} \par\includegraphics[width=5cm]{13755fig3.eps} \end{figure} Figure 3:

Prior $P{(s_{\rm in},s_{\rm out})}$ for the solution configuration of caustic crossing event OGLE 2002-BLG-069 (close+central caustic), with an underlying (uninformative) uniform prior in $\log{t_{\rm E}}$. The red cross shows the location of the found caustic crossing, (1.3,0.3), which falls in a region of relative high probability.

Open with DEXTER
In the text

  \begin{figure} \par\includegraphics[width=9cm, bb= 33 0 760 443]{13755fig4a.eps}... ...e*{2cm}\includegraphics[width=5cm]{13755fig4b.eps}\hspace*{2cm}} \end{figure} Figure 4:

In the top panel, the histogram (blue rectangles) shows the distribution of ${t_{\rm E}}$ found after fitting 788 single-lens microlensing events from OGLE 2006-2007 seasons. The solid black line shows the model prediction of Wood & Mao (2005), which is in good agreement with the data. The bottom panel displays the prior $P{(s_{\rm in},s_{\rm out})}$ for the same intermediate caustic configuration as for Fig. 2 (d=1.1, q=0.1), but assuming an underlying prior for ${t_{\rm E}}$ given by the above distribution.

Open with DEXTER
In the text

Copyright ESO 2010

Current usage metrics show cumulative count of Article Views (full-text article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.

Data correspond to usage on the plateform after 2015. The current usage metrics is available 48-96 hours after online publication and is updated daily on week days.

Initial download of the metrics may take a while.