BINARY MICROLENSING EVENT OGLE-2009-BLG-020 GIVES VERIFIABLE MASS, DISTANCE, AND ORBIT PREDICTIONS

On 2009 February 15, Heliocentric Julian Date (HJD) ∼ 2454878, the Optical Gravitational Lensing Experiment (OGLE) team announced the ongoing microlensing event OGLE-2009-BLG-020 detected by the Early Warning System (EWS)⁵⁵ and observed on the 1.3 m Warsaw Telescope in Las Campanas Observatory in Chile. The event was also monitored by the Microlensing Observations in Astrophysics (MOA) 1.8 m telescope at Mt. John University Observatory in New Zealand. On HJD' ∼ 4915 (HJD' = HJD − 2450000), it could be seen that the light curve was deviating from the standard Paczyński (1986) model, and follow-up observations by other telescopes began: First with the 2.0 m Faulkes South (FTS) telescope in Siding Spring, Australia, and the 2.0 m Faulkes North (FTN) telescope in Haleakala, Hawaii, operated by RoboNet⁵⁶; the 36 cm telescope at Kumeu Observatory,⁵⁷ New Zealand, and the 36 cm telescope at Bronberg Observatory, Pretoria, South Africa, as a part of Microlensing Follow-Up Network (μFUN).⁵⁸ Then, shortly before the first caustic crossing occurred (HJD' ∼ 4917.3), observations began on the 1 m telescope on Mt. Canopus belonging to University of Tasmania (part of PLANET Collaboration) and the 36 cm telescope in Farm Cove Observatory⁵⁹ (also μFUN). Essential data near highest magnification (HJD' ∼ 4917.6, 27 March 27) were gathered by the Bronberg telescope. During the caustic exit (HJD' ∼ 4917.75), there were observations started on the SMARTS 1.3 m Cerro Tololo Inter-American Observatory (CTIO) and the 40 cm telescope in Campo Catino Austral Observatory (CAO),⁶⁰ Chile. Some data were gathered by other observers but too short coverage of the event, or, in some cases, big error bars, prevented us from including these in this analysis. Intense monitoring of the event was performed until HJD' ∼ 4920.

The whole light curve consists of nine years of data with 121 days during the course of the visible magnification, of which 5 days constitute intense follow-up observations. The OGLE telescope performed observations in V and I bands, and the CTIO telescope in V, I, and H bands, which permit measurements of the color of the magnified source star. Together we have 5333 data points in the light curve with 2247 during the course of the magnification event.

We take the OGLE-III V- and I-band light curves from the projects' final data reductions (Udalski et al. 2008) and calibrate against the Galactic bulge photometric maps (Szymański et al. 2011), so that the OGLE magnitudes reported in this paper are standard V (Johnson) and I (Cousins) magnitudes.

The Bronberg data were reduced using MicroFUN's image subtraction pipelines based on Woźniak (2000). The MOA data were reduced using the survey image subtraction pipeline (Bond et al. 2001). CTIO, CAO, FCO, and Kumeu data were reduced using the DoPHOT package (Schechter et al. 1993). The UTAS, FTN, and FTS data were reduced using pySIS3 (Bramich 2008; Albrow et al. 2009).

We remove the effects of a differential extinction from the unfiltered Bronberg data by using light curves of non-variable field stars that have color similar to that of the source (for details see Dong et al. 2009a).

The CTIO telescope's camera, ANDICAM, takes simultaneous exposures in the H and I bands. Because the data quality in the H band is lower than that in the I band, we include only the latter into the microlensing fit. However, we use the H-band data to calculate the (I − H) color of the source. To do so, we cross-match the H-band data with the Two Micron All Sky Survey (2MASS) catalog⁶¹ and calibrate neighboring stars of the lens/source to 2MASS H-band magnitudes.

As the calculation of microlensing magnification is very time-consuming, we bin the data to speed up the fitting procedure. On every part of the light curve and for every individual observatory, we carefully choose the number of bins and their start and end points depending on the local slope of the light curve and specifically avoid incorporating seasonal gap (for wide bins) and daily gaps (for short bins) into any bin. Before the actual binning we fit straight lines to the points in the planned bins and check that the points look consistent with the lines. We also check that the slopes for adjacent bins are comparable. The 3σ outliers from the straight line fits are marked by the program as well as all data points with observational error bars bigger than five times the median of the error bars of nearby points. We remove outliers after visual inspection. If the structure is seen inside the planned bin that might be interpreted as a real deviation, we redesign bins in this region. Our human supervised binning procedure coupled with extensive use of helper algorithms yields 374 binned data points in total.

We rescale the error bars to have χ²_d/dof ∼ 1 for each separate data set (usually coming from different observatories). This is an iterative process wherein we use our best-fit model to evaluate χ²_{d, i} values for each of the data points (where d enumerates data sets and i is the index of the given data point in this set), then based on these values we find parameters needed to rescale the error bars. After rescaling we repeat the fitting procedure, evaluate new χ²_{d, i} values, and find the new values of rescaling parameters.

The error bar rescaling formula is given by $\sigma _{d,i}^{{\rm new}} = Y_d \ssty\sqrt{\sigma _{d,i}^2 + S_d^2}$, where σ_{d, i} and σ^new_{d, i} are the uncertainties before and after rescaling and the coefficients Y_d and S_d are the parameters we find in the following way. For each data set d, we sort all points by the magnification given by the microlensing model. Then we construct the cumulative distribution of χ²_{d, i} as a function of the sorted index. The goal is to choose S_d to make the cumulative distribution of χ²_{d, i} a straight line. Then we can choose the scale Y_d to match the requirement χ²_d/dof ∼ 1.

If we see some systematic variations in the residuals from the best-fit model that are unlikely to be caused by the microlensing phenomena, we investigate whether the deviations are supported by more than one data set or whether it is likely that observations were affected by high airmass or large seeing. If so, we remove these data as outliers.

Figure 1 presents the light curve of the event. All data points are aligned to the best-fit microlensing model. The inset shows the portion of the light curve during highest magnification—when the source is crossing the caustics.

Zoom In Zoom Out Reset image size

The initial mathematical model used to describe this event is constructed using seven parameters: the mass ratio of the lensing binary (q); the projected separation of its components (s₀) in units of the Einstein radius; the angle of the lens–source relative motion projected onto the sky plane with respect to the binary axis (α₀); the Einstein crossing time (t_E), i.e., the time required for the lens to travel a distance of one Einstein radius; the time of the closest approach of the adopted center of the lens to the source (t₀); the lens–source separation at this time in units of the Einstein radius (u₀); and the source radius in the same units (ρ). See Appendix A for conventions and full definition of all parameters. We define the primary/secondary as the heavier/lighter component, respectively.

In addition, there are two parameters for each observatory that connect the microlensing magnification with the instrumental fluxes in a given band. These are the values of the source flux and an additional flux not being magnified but observed in the same direction in the sky (a blend). If there is any light coming from the lens, it will be included in this blend.

To find a microlensing model we use the method developed by S. Dong and described in Dong et al. (2006) and its modification in Dong et al. (2009b, Section 3).

We introduce a few modifications. For example, we calculate the "caustic width" (w) used for the (w, q) grid of lens geometries when searching for initial solutions. Instead of the planetary regime (q ≪ 1), we are now in a stellar binary regime, i.e., q ∼ 1. We can calculate the shorter diameter of the central caustic with the equation

$$\begin{equation} w = \frac{4 q s^2}{ (1+q)^2 }, \end{equation} \tag{ 1 }$$

based on the work by An (2005). Unlike the formula of Chung et al. (2005, Equation (12)) used by Dong et al. (2009b), the above formula is useful for any mass ratio, although it is accurate only to about 20% for lens separations 0 < s < 0.5 and accurate to a factor of 2 for s < 0.7 and stellar mass ratios (above s ∼ 0.7, we are close to intermediate caustic geometry with resonant caustic). Usually, a much better accuracy is not required as one uses the caustic dimensions as a basic scale for further searches, although for our event's geometry and chosen parameterization a more accurate formula is beneficial.

We find that combining the ideas used in both formulae gives us much more accurate results. Explicitly, we divide Equation (12) of Chung et al. (2005) by the factor (1 + q)², making it usable for any mass ratio. This modified recipe is good to about 0.5% for binary separations 0 < s < 0.5 and to about 25% for 0.5 < s < 0.7, as we get closer to the resonant caustic regime. The final formula by which we calculate the shorter diameter of the central caustic is derived to be

$$\begin{equation} w = \frac{4 q |{\rm sin} ^3 \phi _c| }{ (1+q)^2 (s+s^{-1} - 2 \cos \phi _c)^2 }, \end{equation} \tag{ 2 }$$

with the parameter ϕ_c given by Equation (10) of Chung et al. (2005).

For wide lenses (e.g., s > 2 for stellar mass ratios), Equation (2) could be multiplied by $\sqrt{1+q}$ to produce an estimate of the short diameter of the caustic that is usable for our purposes but is a less accurate approximation.

The (1 + q)² factor can be used as well to modify Equation (11) of Chung et al. (2005) formula for calculating the longer diameter of the central caustic. Although it is usable for stellar binary mass ratios, it is accurate only to the factor of 2 for 0 < s < 0.5 (5 for 0.5 < s < 0.7).

In our method, we also modify the parameterization of the Markov chain Monte Carlo (MCMC) χ² minimization by using (u₀/w) rather than u₀ (impact parameter). When we observe a near-cusp caustic crossing or cusp approach in the light curve, we can expect that the value of (u₀/w) will correlate less with the changes in lens geometry given by q and s. This is the case in our model for which the dependence of u₀ on q and s is strong but nearly vanishes when u₀ is divided by w.

First, we fit a seven-parameter microlensing model by (after Dong et al. 2009b) creating a broad grid in three parameters (log q, log w, α) and performing minimization in the remaining four parameters, which are chosen to be t₀, (u₀/w), t_eff(≡ u₀t_E), and t_*(≡ ρt_E). The values of the fluxes (two for each observatory) are calculated analytically with the least-squares method. With this standard parameterization, we are not able to find a satisfactory solution, in a sense where all data sets can be aligned in a coherent way.

We then extend our model by taking into account the Earth's orbital motion during the course of the event. This so-called parallax effect is described with two additional parameters. We use the geocentric parallax formalism (An et al. 2002; Gould 2004), wherein we utilize $\bm {\pi }_{{\rm E}}=(\pi _{{{\rm E}},N},\, \pi _{{{\rm E}},E})$ as model parameters (see Appendix A.2 for definition). The parameters α₀, t₀, and u₀, which previously described a straight trajectory of the lens in front of the source, now describe the trajectory we would see if the Earth's velocity was constant throughout the course of the event and not subject to gravitational acceleration. For the value of this constant velocity, we take the real Earth velocity at some fiducial time t_{0, par} ≡ 4917.52, and since the shift between real source position and the straight trajectory is by definition equal to zero at this chosen time, it is convenient to fix it close to some important features of the light curve; in our case, it is a time between the caustic crossings.

By performing minimization in a similar manner as previously, but with a nine-parameter model, we find a solution capable of explaining the shape of the observed variability. The parameters of this model are gathered in Table 1. We have tested whether there are other solutions emerging from known degeneracies, and we find the corresponding solution through the binary ecliptic degeneracy only (see Equation (16)), with the positive u₀ solution, being slightly preferred. Our best model has a binary separation smaller than the Einstein radius. We find no solution with a "wide" geometry given by the s↔s⁻¹ degeneracy. This degeneracy, which is very important for planetary mass ratios, is broken in this stellar binary case (q = 0.27).

The magnitude of the parallax (π_E) together with the angular size of the Einstein radius (θ_E) yields the total mass of the lens (M) and a distance to the lens (D_l) (e.g., An et al. 2002). After Gould (2000), we have

$$\begin{equation} M_l = \frac{\theta _{{\rm E}}}{\kappa \pi _{{\rm E}}}, \quad \kappa \equiv 4 G/(c^2 \,{\rm AU}) \approx 8.1 \, {\rm mas}\,{M_\odot }^{-1}, \end{equation} \tag{ 3 }$$

$$\begin{equation} D_l = {\rm AU}/\pi _l, \quad \pi _l = \theta _{{\rm E}}\pi _{{\rm E}}+ \pi _s, \end{equation} \tag{ 4 }$$

where we derive θ_E from the source radius (ρ) in units of r_E and the measurement of the angular radius of the source in physical units (θ_*; see Section 4.1.1) with the equation: θ_E = θ_*/ρ. We assume the source to be in the Galactic Bulge at D_s = 8 kpc, and thus the source parallax is π_s = 0.125 mas.

The total mass derived from the best "parallax-only" solution is 0.84 M_☉ and the distance is 0.61 kpc. This is in rough agreement with the position of the lens (assuming that all blended light is coming from the lens) on the color–magnitude diagram (CMD) and the assumption that the binary companions belong to the main sequence.

However, there are still some structures in the residuals of the fit, which we hope will be reduced by extending our model by introducing new parameters describing the orbital motion of the lens.

3.3.1. Linear Orbital Motion Approximation

3.3.2. Full Keplerian Orbit

3.4.1. Priors on Orbital Motion

3.4.2. Priors on Lens Parameters

We also include priors on the lens parameters as expected from the simple Galactic model. Our trial fits to the light curve (with only Keplerian priors) show that the lens has substantial proper motion (of order 11–16 mas yr⁻¹) and is located close to the Earth (0.5–2.4 kpc). These suggest that most of the observed proper motion is due to lens velocity itself (Earth, at the time of the event, is traveling toward the Galactic Bulge, and the effect of the source velocity on the relative proper motion is suppressed by the distance factor). This high proper motion yields linear velocities of 25–150 km s⁻¹. Also, the allowed directions of the lens' proper motion lie between 90° and 180° relative to the Galactic North (i.e., between Galactic East direction and Galactic South direction); this, in connection with high proper motion, yields high probability of a substantial component perpendicular to the disk. This picture is more consistent with Galactic thick disk kinematics rather than thin disk. Thus, in order to impose priors on the lens parameters, we assume that the lens belongs to the thick disk population. This is a less constraining assumption than thin disk, as thick disk allows for grater velocity dispersions.

We follow the approach of Batista et al. (2011, Section 5.1), but using parameters suitable for a thick disk. In the exponential disk we use scale height =0.6 kpc and scale length =2.75 kpc. The probability density of the lens mass is assumed to be proportional to M⁻¹. We expect the velocity dispersion of the thick disk stars to be (40, 55) km s⁻¹ in the North and East Galactic directions, respectively. The expected mean velocity is (0, 200) km s⁻¹, where we account for asymmetric drift of the stars of 20 km s⁻¹ behind the Galactic rotation. The expected velocity of the Galactic Bulge sources is zero and its velocity dispersion is (100, 100) km s⁻¹.

From the disk model we assign prior probability density for certain values of lens distance, mass, and lens–source relative proper motion (D_l, M and $\bm {\mu }$). We translate these priors to microlensing parameters (t_E, θ_E, $\bm {\pi }_{{\rm E}}$) using the Jacobian derived by Batista et al. (2011):

$$\begin{equation} j_{\rm gal} =\left\Vert \frac{\partial (M,D_l,{\bm {\mu }})}{\partial (\theta _{{\rm E}},t_{{\rm E}},{\bm {\pi }}_{{\rm E}})}\right\Vert = \frac{2\pi _{\rm rel} M\mu ^2}{t_{{\rm E}}\theta _{{\rm E}}\pi _{{\rm E}}^2}\frac{D_l^2}{{\rm AU}}, \end{equation} \tag{ 5 }$$

where M = θ_E/κπ_E, D_l = AU/(π_rel + π_s), π_rel = θ_Eπ_E, and μ = θ_E/t_E.

The results from MCMC are presented in Figure 3, where likelihoods obtained from the chain are projected onto the number of 2D planes, one for each pair of parameters. Those likelihoods are weighted by the priors as calculated above. We present plots in a "typically expected" scale for each parameter to show how well each of them was measured. We also provide insets to show the detailed shape of the likelihood contours (in under diagonal panels). The diagonal panel for a given parameter shows the posterior likelihood marginalized along all other dimensions. More exact numerical values of parameters from the region of maximal likelihood can be found in Table 1.

To preserve clarity, Figure 3 does not show symmetric solutions that originate from the static binary ecliptic degeneracy, which is discussed in Section 5.1 and given by Equation (16). These can be seen in the right panel of Figure 5 projected onto the π_{E, N}–γ_⊥ plane.

From OGLE-III photometric maps for the Galactic Bulge (Udalski et al. 2008; Szymański et al. 2011), we construct a CMD of stars within 4' around the location of the event (Figure 2). The position of the centroid of the Red-Clump Giant stars on the CMD is derived to be

$$\begin{eqnarray} & ((V-I), V, I)_{{\rm RC,OGLE}} = (&1.93 \pm 0.02,17.39 \pm 0.04,\nonumber\\ && 15.46 \pm 0.04). \end{eqnarray} \tag{ 6 }$$

Zoom In Zoom Out Reset image size

We take the intrinsic centroid of the red clump (RC) to be M_{I, RC, 0} = −0.25 ± 0.05 (Bennett et al. 2010), and the color (V − I)_{RC, 0} = 1.08 ± 0.06 after Bensby et al. (2010, Section 4.5), we calculate M_{V, RC, 0} = 0.83 ± 0.08 and assume the distance to the Galactic Center to be 8.0 ± 0.3 kpc (Yelda et al. 2011). The location of the observed microlensing event is (l, b) = (1°33'29'', −3°49'21'') in Galactic coordinates. When we take into account the inclination of the Galactic Bar, we infer that the RC stars in the field are slightly closer than the Galactic Center. From Nishiyama et al. (2005) we find a distance modulus smaller by 0.05 mag. This leads to an RC distance modulus of 14.47 ± 0.08, which yields estimates of the reddening along the line of sight:

$$\begin{eqnarray} &(E(V-I), A_V, A_I) = (&0.85 \pm 0.06, 2.09 \pm 0.12, \nonumber\\ &&1.24 \pm 0.10). \end{eqnarray} \tag{ 7 }$$

This gives R_VI = A_V/E(V − I) = 2.4, which is comparable to other estimates of the extinction law toward the Bulge (e.g., Bennett et al. 2010).

We assume that the microlensing source is located behind the dust causing this reddening. Applying this average extinction to the source brightness obtained directly from the microlensing fit to the calibrated OGLE data:

$$\begin{eqnarray} &((V-I), V, I)_{{\rm OGLE}} = (&1.929 \pm 0.002,18.36 \pm 0.04, \nonumber\\ &&16.43 \pm 0.04), \end{eqnarray} \tag{ 8 }$$

we obtain the dereddened color and brightness of the source star:

$$\begin{eqnarray} &((V-I), V, I)_{0,{\rm OGLE}} = (&1.08 \pm 0.06, 16.27 \pm 0.13, \nonumber\\ && 15.19 \pm 0.11). \end{eqnarray} \tag{ 9 }$$

4.1.1. Source Radius

Using color–color relations from Bessell & Brett (1988) we infer that the source is a giant star of the spectral type K and we find (V − K)₀ = 2.50 ± 0.15. We assume solar metallicity for the source star. Then, from Houdashelt et al. (2000), we find that the temperature for a cool solar–metallicity giant star with V − I = 1.08 should be about 4650 K ± 100 K. This is in agreement with Ramírez & Meléndez (2005), who give the temperature equal to 4600 K ± 100 K.

To find the value of the surface gravity we use the Berdyugina & Savanov (1994) empirical relation for G-K giants and subgiants given by

$$\begin{eqnarray} &&\log g = 8.00 \log T + 0.31 M_{0,V} + 0.27 [{\rm Fe/H}]-27.15,\nonumber\\ &&\hbox{with an accuracy of } \pm 0.25\, \mathrm{dex,} \end{eqnarray} \tag{ 10 }$$

where we use M_{0, V} = 1.80 ± 0.10 to obtain log g = 2.71 ± 0.26.

Kervella et al. (2004) calibrated the color–brightness relations for giants. We take their relation (14) with dereddened visual magnitude of the source star V₀ = 16.27 ± 0.13 and (V − K)₀ color found above. The resulting angular radius of the source is

$$\begin{equation} \theta _{*} =(4.45 \pm 0.40)\, \mu \hbox{as.} \end{equation} \tag{ 11 }$$

4.1.2. Limb-darkening Coefficients

Knowing the temperature of the source star and its surface gravity, we find the values of the linear limb-darkening coefficients. We use Table 30 of Claret (2000) with the assumption of solar metallicity and turbulence velocity of 2 km s⁻¹. Uncertainties on these parameters have a very minor influence on the values of limb-darkening coefficients, with uncertainties in temperature being a dominant factor. For bands in which the observations have been performed, we find

$$\begin{equation} u_V = 0.795, \quad u_R = 0.716, \quad u_I = 0.618, \quad u_H = 0.429. \end{equation} \tag{ 12 }$$

For calculations of the magnification of the limb-darkened source, we use coefficients in the form of Afonso et al. (2000, Section 3.1). Their linear limb-darkening coefficient Γ is related to commonly used coefficient u with the relation u = 3Γ/(Γ + 2). Thus, we have

$$\begin{equation} \Gamma _V = 0.721, \quad \Gamma _R = 0.627, \quad \Gamma _I = 0.519, \quad \Gamma _H = 0.334. \end{equation} \tag{ 13 }$$

There is a substantial degeneracy in the microlensing model between the angular velocity of the lens in the plane of the sky (γ_⊥) and the north component of the microlensing parallax (π_{E, N}), which can be seen clearly in Figure 3 (the intersection of row 8 and column 6).

Zoom In Zoom Out Reset image size

This degeneracy can be predicted to be present in a wide range of microlensing events. It is a purely geometrical effect—both rotation of the lens axis and parallax motion of the Earth have similar effects on the apparent source trajectory in the lens plane. This modification in both cases could be described, to the first order, as a curvature of the source trajectory. Hence, the curvature needed to explain observed changes in magnification can be a linear combination of both effects.

Here we stress the importance of this degeneracy. Our model minimization shows that it could be very severe and have a significant impact on the final value of the distance to the lens and its mass. Figure 4 shows two models, both fitting the observed light curve well, with extremely different values of π_{E, N} and γ_⊥, which lead to a very similar effective trajectory of the source on the plane of the lens and finally to similar goodness of fit. Masses and distances derived from these solutions differ by a factor of 2. Note that both trajectories overlap more closely when plotted in units of source radius rather than time.

Zoom In Zoom Out Reset image size

In the general case of microlensing events observed toward the Galactic Bulge and with the parallax signal measured, it is very common that the π_{E, E} component is well measured and π_{E, N} has substantial uncertainty. Looking at Figure 3 of Park et al. (2004) one sees that, if we decompose the parallax vector ($\bm {\pi }_{{\rm E}}$) into components parallel (π_∥) and perpendicular (π_⊥) to the direction of the projected Sun's acceleration (near the peak of the event), actually it is π_∥ that is well measured and the big uncertainties are in π_⊥. Since the Earth's (and thus the Sun's) acceleration vector most of the year lies near the east–west direction, the uncertainty in π_⊥ usually translates to uncertainties in π_{E, N} for Galactic Bulge events.

For this event, the caustic crossing occurred on 2009 March 27, so very close to the vernal equinox when the apex of the Earth motion was almost directly toward the Galactic Bulge. That is why the direction of the Sun's acceleration was very close to west, so π_{E, N} (≈ π_{E, ⊥} in this case) is not well constrained.

Because of this "curvature" degeneracy, we are unable to accurately determine the distance and the mass of the lens from the microlensing fit only. We need to use other known system parameters such as observed magnitudes and colors of the lens (where we assume that all blended light is coming from the lens). We therefore restrict the set of solutions to those that are consistent with theoretical color–color and color–magnitude relations.

4.2.1. Choosing Solutions Consistent with the Theoretical Isochrones

We take a set of Y² (Yonsei-Yale) Isochrones⁶² (Demarque et al. 2004) and from the microlensing fit to the V-, I-, and H-band data, we infer the magnitudes of the blend. Each link in the chain of solutions gives slightly different values of the observed magnitudes of the blend but, for reader information, we quote the most common ones:

$$\begin{eqnarray} &((V-I), I, (I-H))_b =(&1.316 \pm 0.01,15.680 \pm 0.01,\nonumber\\ && 1.409 \pm 0.05). \end{eqnarray} \tag{ 14 }$$

We assume that all non-magnified light comes from the binary lens and that the components of the lens are main-sequence stars. Given the mass ratio q  =  0.27 from the microlensing fit, we neglect the light input from the secondary and check the consistency of the mass of the primary, its distance, and its luminosity against the theoretical isochrones.

For initial tests we take the isochrone with solar metallicity and age. For each solution (link) we calculate the distance and mass of the primary and from the interpolated isochrone we find the theoretical magnitudes in V, I, and H that it should have if seen with no reddening. Then, assuming the slope of the reddening curve A_I/A_V = 0.6 and A_H/A_V = 0.17, we find the amount of interstellar reddening in the V band (A_V) that we should subtract in order to shift the observed magnitudes closest to the theoretical values. We throw out all solutions requiring negative reddening or more reddening than the total amount seen toward the Galactic Bulge source. We also add Gaussian weights to each individual link depending on how close it could be placed to the theoretical magnitudes with respect to the error bars, and a weight corresponding to the value of the reddening assuming that it should be about 1 ± 0.5 mag kpc⁻¹ in V.

The set of solutions with the highest weights lies near dA_V/dD = 1.4 mag kpc⁻¹, distance ≈1 kpc, and mass ≈1 M_☉ for the primary. From the most probable values of parallax parameters (π_{E, N} ≈ −0.2 and π_{E, E} ≈ 0.2), we see that the likely direction of the transverse velocity of the lens relative to the line of sight is perpendicular to the Galactic plane. The value of the relative proper motion taken from the equation μ = θ_E/t_E is near 12 mas yr⁻¹, which at the distance of 1 kpc translates to a linear transverse velocity of ∼60 km s⁻¹. These further constrained kinematic parameters strengthen our prediction (from Section 3.4.2) that it is more likely that the lens belongs to a thick disk population. Thus, for final isochrone consistency checks, instead of a Solar-like isochrone, we choose one with an age of 10 Gyr, a sub-solar metallicity ([Fe/H] = −0.5) and an α-enhanced mixture ([α/Fe] = 0.3) which would be more typical for stars with this kind of kinematics.

We prepare a subset of our chain links in which every link has weight added depending on how well it matches with the theoretical isochrone. The region of parameter space that coincides with the magnitudes from the isochrone has

$$\begin{eqnarray} D_l &=& 1.1 \pm 0.1 \, {\rm kpc}, \quad M_1 = 0.84 \pm 0.03 \, M_\odot, \nonumber\\ &&\frac{{d}A_V}{{d}D} = 1.4 \pm 0.2 \, {\rm mag}\,{\rm kpc}^{-1}. \end{eqnarray} \tag{ 15 }$$

4.2.2. Choosing Solutions Consistent with the Spectroscopic Mass

Nine of the microlensing parameters are known from the fit with very high precision. However, there are two pairs of parameters for which the values are imperfectly measured (see Figure 3). The first, obvious, pair is the position and velocity along the line of sight (s_z and γ_z), which do not have direct effects on the observed magnification. They are only limited by the requirement of e < 1 and by the shape of the Keplerian orbit projected on the plane of the sky—to be precise, only by the segment of the orbit that was covered by the lens components during the high-magnification period. The other pair of parameters that was not measured well are π_{E, N} and γ_⊥, which are degenerate as described in Section 4.2. We partially break this degeneracy into discrete regions by choosing only those solutions that are in agreement with other information we have on the event (Sections 4.2.1 and 4.2.2).

Because the mass of the lens depends on the magnitude of the vector $\bm {\pi }_{{\rm E}}= (\pi _{{{\rm E}},N}, \pi _{{{\rm E}},E}),$ it does not change when the sign of the π_{E, N} component changes. This leads to two regions along the π_{E, N}–γ_⊥ degeneracy that yield the same mass.

There is one more degeneracy from which our model suffers, i.e., the orbiting binary ecliptic degeneracy (see Appendix A.4). We note that the ecliptic degeneracy happens when the direction of the Sun's acceleration is constant: the one obvious case is when the source lies on the ecliptic (hence the name of this degeneracy), but the other is when the timescale of the event is much shorter than the characteristic timescale of the changes of the direction of acceleration. This is more likely to occur if the microlensing event happens when Earth is moving toward or away from the source.

This is the case in OGLE-2009-BLG-020, for which the peak magnification happened near the vernal equinox. The Sun's acceleration at that time is directed toward West and its direction varies very slowly in time. Note also that at that time π_{E, ⊥} ≈ π_{E, N} so, following (A16), from any set of parameters we can obtain an analogous set by changing

$$\begin{equation} (u_0,\alpha _0,\pi _{{{\rm E}},N},\gamma _\perp) \rightarrow -(u_0,\alpha _0,\pi _{{{\rm E}},N},\gamma _\perp). \end{equation} \tag{ 16 }$$

In Figure 5, we see the projection of our set of solutions on the plane of π_{E, N} and γ_⊥. The two allowed regions correspond to positive u₀ (on top) and negative u₀ (bottom) solutions. The left panel shows the full chain while the right panel shows the subset of links that are consistent with the isochrones chosen in Section 4.2.1. Each of the continuous regions from the right panel is split into two regions (for positive π_{E, N} and negative π_{E, N}), giving four separate regions in total.

Zoom In Zoom Out Reset image size

Note that none of the degeneracies described in this section is exact, which leads to different likelihoods between the four regions.

Changing the sign of γ_⊥ without changing the other components of the velocity flips the orbit, which manifests itself as a change of some Euler angles: Ω_node → −Ω_node and i → π − i.

As microlensing cannot directly measure position and velocity along the line of sight, every solution has its copy mirrored by the plane of the sky, i.e., s_z → −s_z and γ_z → −γ_z. For orbit elements, it means Ω_node → π − Ω_node and ω_peri → ω_peri − π.

Without loss in generality, we perform all calculations assuming γ_z > 0 at the chosen time (t_{0, kep}). RV measurements will have information on the sign of the radial velocity—we can then simply mirror our set of solutions using the above prescription.

Every link in the MCMC chain represents a complete set of parameters of the binary lens. It yields not only the mass, distance, and separation in physical units, but also all Keplerian parameters of the orbit. Thus, we can calculate the RV curve for any period of the binary.

Assuming that the observations of the RV curve of the primary are taken, we can assign to every link the likelihood that it is consistent with those data. We derive the radial velocity for any time the data were taken and, allowing the systemic velocity of the center of mass of the binary to vary, we calculate the likelihoods for every point. In this way, we construct a new set of links weighted by both the microlensing light curve and the RV curve. This new set of solutions yields new, most probable, values of all lens parameters, which may or may not coincide with the values derived from the microlensing-only solution. This would be a test of the microlensing solution.

If the new values of the parameters lie outside the 3σ limits of the microlensing solution, we can say that our method failed. By contrast, if solutions, which are consistent with the RV curve, lie near the best-fit values we obtained from microlensing, we can not only believe our solution but we can also read off all parameters of the binary, which are not given by one or the other method alone. For example, the RV curve will give us the period and systemic radial velocity, which we cannot read from the microlensing light curve alone. However, microlensing will yield the inclination, orientation on the sky, and the 2D velocity of the binary projected on the plane of the sky.

5.3.1. Example of the Test

We illustrate this test in Figure 6. The figure shows (in the background) the likelihoods derived from the microlensing solution projected onto 2D planes of six orbital parameters, the mass, distance, and mass ratio (this is a subset of the whole chain chosen using the method described in Section 4.2.1). We take one exemplary set of binary parameters and simulate the RV curve consisting of 15 measurements taken from 2011 March to October and 15 measurements in a similar period in 2012. We then test our chain against this RV curve assuming a measurement precision of 0.5 km s⁻¹. The solutions with the highest joint likelihoods (from microlensing and RV curve fitting) are overplotted on each panel (in color on under-diagonal panels, and as 10σ contour on over-diagonal panels). The initial set of parameters we chose to generate the RV curve is marked with open circles.

Zoom In Zoom Out Reset image size

We see that all orbital parameters of the "seed" set of parameters are retrieved by the RV comparison process (as the links selected to be consistent with the RV curve lie near the circles). This shows the selective power of the RV data, which can prove whether the microlensing solution is wrong or strongly suggest that it is right.

5.3.2. Useful Information for RV Observations

For point-lens events, the almost universally accepted three basic parameters are

$$\begin{equation} (t_0,u_0,t_{{\rm E}})\quad \rm (basic\ point\ lens). \end{equation} \tag{ A1 }$$

These are the time of closest approach to the lens system "center," the lens–source projected separation at that time (in units of the Einstein radius), and the Einstein crossing time. The meaning of system "center" is obvious in the case of a point lens, but will require generalization for more complicated systems. Moreover, even t_E will require more exact specification. We retain these three notations, but defer discussion of the generalizations of their meaning until later on in this paper. A derived parameter, which is sometimes used as an independent fitting parameter in place of either u₀ or t_E, is the "effective timescale"

$$\begin{equation} t_{\rm eff} \equiv u_0 t_{{\rm E}}. \end{equation} \tag{ A2 }$$

There are three higher-order parameters that can in principle be measured for point-lens events, and these lead to a fourth derived parameter. The first two parameters are the "vector microlens parallax"

$$\begin{equation} {\bm {\pi }}_{{\rm E}}\equiv (\pi _{{{\rm E}},N},\pi _{{{\rm E}},E})\equiv (\cos \phi _\pi,\sin \phi _\pi)\pi _{{\rm E}}. \end{equation} \tag{ A3 }$$

Here, $\pi _{{\rm E}}={\rm AU}/\tilde{r}_{{\rm E}}$, where $\tilde{r}_{{\rm E}}$ is the Einstein radius projected onto the observer plane and ϕ_π is the direction of the lens motion relative to the source expressed as a counterclockwise angle north through east. The mere introduction of ${\bm {\pi }}_{{\rm E}}$ brings with it a large number of symmetries and questions of convention, which will grow yet more complicated as binary lenses come into play. We therefore carefully delineate these in their simpler form here.

First, it has become customary to adopt the "geocentric framework" (An et al. 2002; Gould 2004), in which all parameters are measured in the instantaneous frame that is at rest with respect to the Earth at a specifically adopted time:

$$\begin{equation} t_{0,\rm par}\quad \rm (parameter\ reference\ time), \end{equation} \tag{ A4 }$$

which is not a fit parameter. Note that the subscript stands for "parameter" (not "parallax"). Thus, for example, u₀ is the distance of closest approach of the source and the "lens center" (i.e., for point lenses, simply the lens) in this geocentric frame. Similarly, t₀ is the time of this closest approach, and t_E is the time it would take to cross the Einstein radius if the lens–source relative motion were the same as it is as seen from the Earth at t_{0, par}. Since the Earth velocity is constantly changing, all of these parameters depend on the choice of t_{0, par}. Finally, the direction ϕ_π (but not the magnitude π_E) of the parallax vector also depends on this choice. In practice, for point lenses, t_{0, par} is chosen quite close to t₀ so that this is hardly an issue. But the issue will become more important for binary lenses.

Second, point-lens parallaxes are subject to four related "degeneracies," which one might also call "symmetries." To properly express these, we introduce a different set of basis vectors ${\bm {\pi }}_{{\rm E}}= (\pi _{{{\rm E}},\parallel },\pi _{{{\rm E}},\perp })$, where π_{E, ∥} is the component parallel to the apparent acceleration of the Sun (projected on the sky) in the Earth frame at t_{0, par} and the pair (π_{E, ∥}, π_{E, ⊥}) is right-handed (Figure 3 of Gould 2004). We must also define a convention for the sign of u₀:

$$\begin{equation} u_0>0 \quad \Leftrightarrow\quad \rm (lens\ passes\ source\ on\ its\ right), \end{equation} \tag{ A5 }$$

as in Gould (2004, Figure 2). Then:

• 1.  
  π_{E, ⊥} Degeneracy. Typically, π_{E, ∥} is much better determined than π_{E, ⊥} (Gould et al. 1994) because the former is determined at third order in time and the latter at fourth order (Smith et al. 2003). This can lead to elongated error ellipses in the ${\bm {\pi }}_{{\rm E}}$ plane. Only for sufficiently long events (or strong parallax signal) is this degeneracy broken.
• 2.  
  u₀ Degeneracy. In the limit that the projected acceleration of the Sun can be regarded as constant, there is a perfect symmetry between ±u₀ solutions, with only minor adjustment of other parameters (Smith et al. 2003).
• 3.  
  Ecliptic Degeneracy. In the limit of constant direction of acceleration (as would be the case for a source on the ecliptic), there is a perfect degeneracy

  $$\begin{equation} (u_0,\pi _{{{\rm E}},\perp }) \rightarrow -(u_0,\pi _{{{\rm E}},\perp }) \end{equation} \tag{ A6 }$$

  (Jiang et al. 2004; Poindexter et al. 2005). Hence, sources lying near the ecliptic (i.e., all Galactic Bulge sources) may suffer an approximate degeneracy. We note that at times when the Earth is moving toward or away from the Galactic Bulge (near equinoxes), the direction of the Sun's acceleration projected on the sky varies slowly in time. This significantly extends the effect of this degeneracy.
• 4.  
  Jerk-Parallax Degeneracy. Gould (2004) discovered a fourfold degeneracy in point-lens events with parallax. This degeneracy can be broken partially or completely depending on circumstances and features of the light curve. Two solutions out of four that are more likely to survive are (u₀, π_{E, ⊥}) → −(u₀, π_{E, ⊥} + π_{j, ⊥}). For the case when the source is on the ecliptic, this is equal to the ecliptic degeneracy mentioned above.

The third higher-order point-lens parameter is ρ, the radius of the source in units of the Einstein radius. Of course, ρ is only rarely measurable in point-lens events and is very frequently measured in binary events, and so is usually called a "binary-lens" parameter, but it logically precedes the introduction of binary events, so we include it here. The angular radius of the source θ_* can almost always be measured from the instrumental color and magnitude of the source (Yoo et al. 2004). When ρ is also measurable, then one can measure the angular Einstein radius, and so the proper motion:

$$\begin{equation} \theta _{{\rm E}}= \frac{\theta _*}{\rho },\qquad \mu = \frac{\theta _{{\rm E}}}{t_{{\rm E}}}. \end{equation} \tag{ A7 }$$

If both θ_E and π_E are measured, this immediately gives both the lens mass and the lens–source relative parallax (Gould 1992)

$$\begin{equation} M = \frac{\theta _{{\rm E}}}{\kappa \pi _{{\rm E}}},\quad \pi _{\rm rel} = \theta _{{\rm E}}\pi _{{\rm E}}, \quad \kappa \equiv \frac{4G}{c^2\,\rm AU}. \end{equation} \tag{ A8 }$$

Note that the "proper motion" in Equation (A7) is geocentric (i.e., in the instantaneous frame of the Earth at t_{0, par}). To compare with heliocentric proper motions derived from astrometry at successive epochs, one must convert (Janczak et al. 2010)

$$\begin{equation} {\bm {\mu }}_{\rm helio} -{\bm {\mu }}_{\rm geo} = \frac{{\bf v}_{\oplus,\perp }}{{\rm AU}} \pi _{\rm rel} = \frac{{\bf v}_{\oplus,\perp }}{{\rm AU}} \frac{\theta _{{\rm E}}^2}{\kappa M}, \end{equation} \tag{ A9 }$$

where v_{⊕, ⊥} is the Earth velocity at t_{0, par} projected on the plane of the sky, and ${\bm {\mu }}_{\rm geo} = \mu {\bm {\pi }}_{{\rm E}}/\pi _{{\rm E}}$. Clearly, this conversion can only be carried out precisely if ${\bm {\pi }}_{{\rm E}}$ is known although the right-hand side shows that the magnitude of the difference can be strongly constrained even if there are only fairly crude limits on π_E.

Note also that future astrometric microlensing measurements may determine the (geocentric) vector proper motion ${\bm {\mu }}_{\rm geo}$ even if no parallax information is obtained (Høg et al. 1995; Miyamoto & Yoshii 1995; Walker 1995).

At this point, we also introduce

$$\begin{equation} t_*\equiv \rho t_{{\rm E}}, \end{equation} \tag{ A10 }$$

the source self-crossing time (in fact, a source radius crossing time). Often, t_* is used in place of ρ as a fitting parameter. Note that sometimes ρ is written as ρ_* (in analogy to t_*), but for ρ there is no need to subscript because there are no competing quantities with this same name. Hence, we advocate simplifying the notation and dropping the subscript "*."

A static binary requires three additional parameters:

$$\begin{equation} (q,s_0,\alpha _0) \quad \rm (static\ binary\ parameters). \end{equation} \tag{ A11 }$$

These are the mass ratio of the two components (q = m₂/m₁ = M₂/M₁), their projected separation in units of the Einstein radius, and the direction of lens–source relative motion (i.e., lens motion relative to the source) with respect to the binary axis (which points from primary toward secondary). At the beginning, one should specify which component of the binary one calls a "primary." (The angle α₀ is counterclockwise. The fractional mass of the primary, m₁, is defined as M₁/M.) There are many points to note.

First, the separation is frequently called "d " (rather than "s"). However, this is the standard symbol for the derivative operator, and so should not be used for other quantities that are likely to appear in the same expressions with this operator.

Second, heretofore, s₀ and α₀ have been written simply as s and α. And indeed, for static binaries there is no confusion in doing so, since these are time-invariant quantities. And, for this reason, it remains appropriate to drop the "0" subscript in static analysis. Nevertheless, we introduce this subscript here to maintain consistency with notation developed below.

Third, the choice of t_{0, par} is no longer obvious. For high-magnification events, it might be taken as the approximate time of closest approach to the center of magnification, which closely parallels the choice for point-lens events. But it might also be chosen to be a particularly well-defined time, like a caustic crossing. In any case, since it is not even approximately obvious, it must be specified.

Fourth, in generalizing from point lens to binary lens, there is no longer a unique system "center" by which to define u₀ and t₀. However, since this does not present substantive problems until there is orbital motion, we defer discussion of this point until the next section.

Fifth, in the absence of parallax effects, static binaries are subject to an exact degeneracy

$$\begin{equation} (u_0,\alpha _0)\rightarrow -(u_0,\alpha _0) \quad \rm (binary\ degeneracy). \end{equation} \tag{ A12 }$$

Moreover, even for orbiting binaries, it is always possible to express solutions in a form with u₀ > 0 and 0 ⩽ α₀ < 2π. Hence, in our view, negative u₀ should be reserved for solutions that include parallax.

Finally, even if parallax is incorporated into the solution, static binaries will be subject to a "static binary ecliptic degeneracy":

$$\begin{eqnarray} &&(u_0,\alpha _0,\pi _{{{\rm E}},\perp })\rightarrow -(u_0,\alpha _0,\pi _{{{\rm E}},\perp }) \nonumber\\ &&\rm (static\ binary\ ecliptic\ degeneracy). \end{eqnarray} \tag{ A13 }$$

Note that for lenses seen toward the Galactic Bulge, the directions of positive (π_{E, ∥}, π_{E, ⊥}) are typically (west, north) for austral summer and autumn or (east, south) for austral winter and spring. Hence, from a practical standpoint, it is often easy to locate degenerate π_{E, ⊥} solutions by seeding with a sign reversal of π_{E, N}.

We begin by introducing the components of orbital motion

$$\begin{equation} {\bm {\gamma }}\equiv (\gamma _\parallel,\gamma _\perp) \quad {\rm (at\ }t_{0,{\rm kep}}\rm), \end{equation} \tag{ A14 }$$

which are the instantaneous components of velocity of the secondary relative to the primary, respectively parallel and perpendicular to the primary–secondary axis at a fiducial time t_{0, kep}. These are essentially γ_∥ = (ds/dt)/s₀ and γ_⊥ = −dα/dt, and they are detected effectively from the changing shape and changing orientation of the caustic, respectively. In this form, their relation to the projected physical orbital velocity is particularly simple,

$$\begin{equation} \Delta {\bf v}= D_l\theta _{{\rm E}}s_0{\bm {\gamma }}, \end{equation} \tag{ A15 }$$

where D_l is the distance to the lens. Note that (γ_∥, γ_⊥) is a right-handed system on the plane of the sky, just like (N, E) and (π_{E, ∥}, π_{E, ⊥})

The introduction of these two parameters brings with them two degeneracies. First, allowing for orbital motion reintroduces the π_{E, ⊥} continuous degeneracy. This appeared originally for point lenses because of rotational symmetry, which is the physical reason that π_{E, ⊥} was fourth-order in time while π_{E, ∥} was third-order (see Appendix A.2). Now it reappears because of a degeneracy between π_{E, ⊥} and the rotational degree of freedom of orbital motion, i.e., as a correlation between π_{E, ⊥} and γ_⊥ (cf. Section 4.2). Second, to the degree that the ecliptic degeneracy is present (i.e., to the degree that the Sun's acceleration does not change direction during the event), it takes the form

$$\begin{eqnarray} &&(u_0,\alpha _0,\pi _{{{\rm E}},\perp },\gamma _\perp) \rightarrow -(u_0,\alpha _0,\pi _{{{\rm E}},\perp },\gamma _\perp)\nonumber\\ && \rm (orbiting\ binary\ ecliptic\ degeneracy). \end{eqnarray} \tag{ A16 }$$

Logically, t_{0, par} and t_{0, kep} can be different, and so we have defined them separately. There may be cases for which one would want them to be different, but we cannot think of any. Therefore, we suggest in general that

$$\begin{equation} t_{0,\rm par}\equiv t_{0,\rm kep}, \quad \rm (suggested). \end{equation} \tag{ A17 }$$

We will adopt this convention in what follows. (This will also help to avoid confusion involving α₀, which is defined at t_{0, par}, and s₀, which, together with other phase-space parameters of the orbit, is defined at t_{0, kep}.)

We now return to the problem of specifying a "system center." For simplicity, let us begin by assuming that the binary center of mass is chosen. Then, the problem of determining t₀ and u₀ for a given model is identical to that of a single lens (as discussed in Appendix A.2) except that t_{0, par} may then differ substantially from t₀ (since it is convenient to choose t_{0, par} near the highest magnification), but this poses no problem of principle. And, in particular, if s < 1, then the center of mass and center of magnification coincide, which greatly simplifies the choices. However, if s > 1, one might for example choose the system center to be the "center of magnification" and t_{0, par} to be approximately the time closest to this center. Then, the difference between t_{0, par} and t₀ would be small, but "center of magnification" would not, in general, be in rectilinear motion (together with the center of mass). In this example, the "system center" would be offset from the center of mass by

$$\begin{equation} \Delta (\xi,\eta)= \biggl [\frac{q}{1+q}\biggl (\frac{1}{s(t_{0,{\rm par}})} - s(t_{0,{\rm par}})\biggr),0\biggr ], \end{equation} \tag{ A18 }$$

where (ξ, η) are the coordinates on the lens plane parallel and perpendicular to the primary–secondary axis. Note that this offset properly describes the true position of the "center of magnification" only at the one, chosen time (t_{0, par}) and not at the fit time t₀. Again, for the case just specified, these two times would be very close, so there is little practical impact. However, in other cases, particularly for resonant caustics of roughly equal-mass binaries, t_{0, par} might be taken to be near a caustic crossing that is very far from t₀. Then, specification of this time becomes extremely important.

For binary lenses with higher-order parameter measurements, one almost always measures θ_E. And, because γ_⊥ is often degenerate with π_{E, ⊥}, it is usually (but not always) inappropriate to attempt to measure ${\bm {\gamma }}$ without also measuring ${\bm {\pi }}_{{\rm E}}$. This means that when ${\bm {\gamma }}$ is measured, one can usually calculate the ratio of kinetic to potential projected energy (Batista et al. 2011),

$$\begin{equation} \frac{E_{\perp,{\rm kin}}}{E_{\perp,{\rm pot}}} = \frac{\kappa M_\odot \pi _{{\rm E}}(|{\bm {\gamma }}|{\rm yr})^2 s_0^3}{ 8\pi ^2\theta _{{\rm E}}(\pi _{{\rm E}}+ \pi _s/\theta _{{\rm E}})^3}. \end{equation} \tag{ A19 }$$

Evaluation of this ratio requires specification of one additional parameter, the source parallax π_s. This is usually known, at least approximately, since most sources are in the Galactic Bulge.

This ratio must absolutely be less than unity, if the system is bound. Moreover, typically it is expected to be in the range 0.25–0.6. Hence, evaluating this ratio provides a good plausibility check on the solution. If θ_E is measured in milliarcseconds and ${\bm {\gamma }}$ is measured in yr⁻¹ (as we advocate), then the numerical coefficient in Equation (A19) is 8.14/(8π²) = 0.103.

Full Kepler solutions require two parameters in addition to those already mentioned. Since the parameters already specified are in the Cartesian system, we advocate making the last two parameters also Cartesian, namely the instantaneous position and velocity in the direction perpendicular to the plane of the sky at time t_{0, kep}. Specifically, these are s_z and γ_z, which are defined so that the physical relative 3D position and velocity of the secondary relative to the primary are given by

$$\begin{equation} \Delta {\bf r} = D_l\theta _{{\rm E}}(s_0,0,s_z), \quad \Delta {\bf v} = D_l\theta _{{\rm E}}s_0(\gamma _\parallel,\gamma _\perp,\gamma _z) \end{equation} \tag{ A20 }$$

(see Figure 7). Recall that the x-axis is defined by the position of the secondary relative to the primary at t_{0, kep}. These are right-handed triads, which means that the radial direction points toward the observer, which is opposite to the convention of RV work. This conflict is not desirable but is virtually unavoidable. The identification of the x-axis with the binary axis is extremely firmly rooted in microlensing tradition. There is then only one way to maintain right-handed 2D and 3D systems for ${\bm {\gamma }}$: the one adopted above. If we were to choose a left-handed system, it would give rise to substantial confusion in the calculation of orbits, which involve several cross products. Finally, it will be very rare that any RV data are obtained on microlensing orbiting binaries. When they are, one will just have to remember to reverse the sign.

Zoom In Zoom Out Reset image size

Note that full Keplerian solution is possible for events with π_E and θ_E measured. For other events one can always use simplified 2D description of the orbital motion (Appendix A.4), assuming that the fraction of the orbit traveled by the binary components during the event is small.

There is one final degeneracy associated with two introduced (radial) parameters. As with visual binary orbits, there is a perfect degeneracy of

$$\begin{equation} (s_z,\gamma _z) \rightarrow -(s_z,\gamma _z), \end{equation} \tag{ A21 }$$

which can only be broken with the aid of RV measurements. Hence, microlensing solutions should choose either s_z ⩾ 0 or γ_z ⩾ 0. The choice will depend on which quantity is more "isolated" from zero when both are permitted free range. However, this choice should definitely be made to avoid multiple representations of what are in fact identical solutions.

Point-lens events are specified by up to six parameters $(t_0,u_0,t_{{\rm E}},{\bm {\pi }}_{{\rm E}},\rho)$, where ${\bm {\pi }}_{{\rm E}}=(\pi _{{{\rm E}},N},\pi _{{{\rm E}},E})$ is the 2D parallax vector. This is always a geocentric system, explicitly if ${\bm {\pi }}_{{\rm E}}$ is specified and implicitly if it is not. Because ${\bm {\pi }}_{{\rm E}}$ is a geocentric quantity that depends on the Earth's velocity, the definition of all of these quantities (except ρ) requires that a specific time t_{0, par} be adopted at which the Earth's velocity and position are evaluated. The only exception is if ${\bm {\pi }}_{{\rm E}}$ is itself not specified. These requirements carry over to binary lenses as well.

If ρ is measured, then it is almost always possible to measure θ_E = θ_*/ρ because θ_* can be determined from the source flux parameters in two bands.

A total of eight parameters are required to specify a binary orbit. Two of these are the secondary/primary mass ratio q and the total mass of the system M = θ_E/κπ_E. The remaining six parameters are the six Cartesian phase–space coordinates of the secondary relative to the primary. In order to relate microlensing parameters to physical units, one must multiply by r_E ≡ D_lθ_E. The lens distance can be expressed in terms of parameters already specified, plus the source parallax: D_l = AU/(π_Eθ_E + π_s). Then, all but one of the six phase–space coordinates appear directly as microlensing parameters (s₀, 0, s_z) and ${\bm {\gamma }}=(\gamma _\parallel,\gamma _\perp,\gamma _z)$, for the positions and velocities, respectively, with all quantities specified at t_{0, kep}. The missing degree of freedom represented by "0" in the spatial vector is recovered in α₀, which specifies the direction of lens–source relative motion, measured with respect to the primary–secondary axis at time t_{0, par}. We recommend $t_{0,\rm kep}= t_{0,\rm par}$.

For binary lenses, the "system center" must be given explicitly, in order to define t₀, u₀, t_E. This center must be defined as coordinates on the lens plane as determined at t_{0, par}, not at t₀; see Appendix A.1.

All coordinate systems are right-handed, either in two dimensions [(N, E), (π_{E, ∥}, π_{E, ⊥}), (γ_∥, γ_⊥)] or three dimensions (γ_∥, γ_⊥, γ_z). This means, in particular, that positive γ_z points toward the observer.

All relative motion conventions are defined by the motion of the lens (with the source thought of as fixed). Thus, first, the sign of u₀ is positive if the lens passes the source on its right. Second, ϕ_π = atan2(π_{E, E}, π_{E, N}) is the angle of lens motion, measured counterclockwise relative to north. And, third, α₀ is the angle of the lens motion, measured counterclockwise relative to the primary–secondary axis.

There are two approximate degeneracies, one discrete and one continuous, which propagate through the analysis as more parameters are added. In their "final form" when binary orbital motion is modeled, the discrete degeneracy is a generalized "ecliptic degeneracy," which reverses the sign of the four parameters (u₀, π_{E, ⊥}, α₀, γ_⊥). The continuous degeneracy is between π_{E, ⊥} and γ_⊥ (with u₀ maintaining the same sign). These degeneracies persist in simpler form when some of these parameters are not measured. Finally, there is a perfect degeneracy that reverses the signs of (s_z, γ_z), so that the solutions must be restricted to either s_z ⩾ 0 or γ_z ⩾ 0 to avoid duplicating identical solutions.

We have two systems of coordinates, the first one is relative to the plane of the orbit and the second is related to the microlensing event. The microlensing system, as described in Appendix A, has its first axis set by the binary axis projected onto the plane of the sky, and the third axis toward the observer. For fitting, we use positions (s) and velocities ($\bm {\gamma }$) in microlensing units. However, for the sake of this appendix we will be using positions and velocities (r, v) in physical units of AU and AU yr⁻¹. A simple conversion between these units is given by r = s · D_lθ_E and ${\bf v} = r_x \bm {\gamma }$, where D_l and θ_E can be calculated at any time from the set of microlensing parameters.

In the system of coordinates relative to the orbit, we set the first axis as the direction of periastron, and the third axis as that of the angular momentum vector. The position of the body (r') in this system is given by

$$\begin{equation} {\bf r}^{\prime } = \left(\begin{array}{@{}c@{}} \cos {\mathcal {E}}- e\\ \sqrt{1-e^2}\sin {\mathcal {E}}\\ 0 \end{array} \right) a \quad r^{\prime } = a(1-e \cos {\mathcal {E}}), \end{equation} \tag{ B2 }$$

where a is the semi-major axis, e is the eccentricity, and ${\mathcal {E}}$ is the eccentric anomaly, which is defined as an implicit function of time t by the Kepler equation,

$$\begin{equation} n (t-t_{\rm peri}) = {\mathcal {E}}- e \sin {\mathcal {E}}, \qquad n\equiv \sqrt{\frac{{\rm GM}}{a^3}}. \end{equation} \tag{ B3 }$$

Here, M is the system mass, n is the mean motion, and t_peri is the time of periastron.

The velocity in the plane of the orbit (v') is given by the derivative of position:

$$\begin{eqnarray} {\bf v}^{\prime } &=& \frac{{d}{\bf r}^{\prime }}{{d}t} = \left(\begin{array}{@{}c@{}} {-} {\rm sin}\ {\mathcal {E}} \\ \sqrt{1-e^2} \cos {\mathcal {E}}\\ 0 \end{array} \right) a \frac{{d}{\mathcal {E}}}{{d}t},\nonumber\\ &&\frac{{d}{\mathcal {E}}}{{d}t} = \frac{n}{1-e\cos {\mathcal {E}}}, \end{eqnarray} \tag{ B4 }$$

where the last result is found by implicit differentiation of Equation (B3). Implicit differentiation also yields

$$\begin{equation} \frac{\partial {\mathcal {E}}}{\partial t_{\rm peri}} = \frac{-n}{\sin {\mathcal {E}}}\frac{\partial {\mathcal {E}}}{\partial e} = \frac{2 a}{3\Delta t}\frac{\partial {\mathcal {E}}}{\partial e} = -\frac{{d}{\mathcal {E}}}{{d}t}. \end{equation} \tag{ B5 }$$

The equations for r' and v' will be useful for calculating relative positions of the binary component at any time and for evaluating the Jacobian of the transformation between the phase-space coordinates used in the microlensing fit and the orbital elements (j_kep). This in turn will be used to construct the full Jacobian (B1) used for weighting the solution, together with the priors on all "physical" parameters (cf. Section 3.4.1):

$$\begin{eqnarray} &j_{\rm kep}^{-1} =&\left\Vert \frac{\partial (s_0, \alpha _0, s_z, \gamma _\parallel, \gamma _\perp, \gamma _z)}{\partial (e, a, t_{\rm peri}, \Omega _{\rm node}, i, \omega _{\rm peri})}\right\Vert =(D_l \theta _{{\rm E}})^{-6} s_0^{-4}\nonumber\\ &&\times\left\Vert \frac{\partial (r_x, r_y, r_z, v_x, v_y, v_z)}{\partial (e, a, t_{\rm peri}, \Omega _{\rm node}, i, \omega _{\rm peri})}\right\Vert = \frac{j_{{\rm ph}/{\rm kep}}}{(D_l \theta _{{\rm E}})^6 s_0^4}.\nonumber\\ && \end{eqnarray} \tag{ B6 }$$

To calculate j_ph/kep, we need to evaluate derivatives of every component of (r, v) with respect to every orbital element. To find ∂r/∂(e, a, t_peri), we will calculate ∂r'/∂(e, a, t_peri) in the plane of the orbit and then rotate it to the microlensing system of coordinates using the rotation matrix ${\mathbb {R}}$ given by the three orbital Euler angles (Ω_node—longitude of the ascending node, i—inclination, ω_peri—argument of periapsis):

$$\begin{equation} {\bf r}={\mathbb {R}}{\bf r}^{\prime },\qquad {\bf v}={\mathbb {R}}{\bf v}^{\prime }, \end{equation} \tag{ B7 }$$

where

$$\begin{equation} {\mathbb {R}}= {\mathbb {R}}_z(\Omega _{\rm node}) \cdot {\mathbb {R}}_x(i) \cdot {\mathbb {R}}_z(\omega _{\rm peri}) \end{equation} \tag{ B8 }$$

and ${\mathbb {R}}_x(\beta)$, ${\mathbb {R}}_z(\beta)$ are the matrix operators of rotation by an angle β around the first and third axes, respectively,

$$\begin{eqnarray} {\mathbb {R}}_x(\beta)=\left(\begin{array}{@{}ccc@{}} \cos \beta & -\sin \beta & 0 \\ \sin \beta & \cos \beta & 0 \\ 0 & 0 & 1 \end{array}\right)\nonumber\\ {\mathbb {R}}_z(\beta)=\left(\begin{array}{@{}ccc@{}} 1 & 0 & 0\\ 0 & \cos \beta & -\sin \beta \\ 0 & \sin \beta & \cos \beta \end{array}\right). \end{eqnarray} \tag{ B9 }$$

Figure 8 shows the orientation of the relative orbit with respect to the microlensing coordinates system. Conventions used for Euler angles (Ω_node, i, ω_peri) can also be read from this figure.

Zoom In Zoom Out Reset image size

From Equation (B5) we obtain

$$\begin{equation} \frac{\partial {\bf r}^{\prime }}{\partial t_{\rm peri}} = - {\bf v}^{\prime }, \qquad \frac{\partial {\bf v}^{\prime }}{\partial t_{\rm peri}} = - {\bf a}^{\prime } \equiv \frac{GM}{(r^{\prime })^3}{\bf r}^{\prime }, \end{equation} \tag{ B10 }$$

$$\begin{equation} \frac{\partial {\bf r}^{\prime }}{\partial a} = -\frac{3\Delta t}{2a}{\bf v}^{\prime } + \frac{{\bf r}^{\prime }}{a}, \qquad \frac{\partial {\bf v}^{\prime }}{\partial a} = -\frac{3\Delta t}{2a}{\bf a}^{\prime } -\frac{1}{2} \frac{{\bf v}^{\prime }}{a}, \end{equation} \tag{ B11 }$$

$$\begin{equation} \frac{\partial {\bf r}^{\prime }}{\partial e} = \frac{\sin {\mathcal {E}}}{n}{\bf v}^{\prime } - \left({\begin{array}{@{}c@{}} 1\\\ms \dsty\frac{e}{\sqrt{1-e^2}}\sin {\mathcal {E}}\\\ms 0\end{array}}\right)a, \end{equation} \tag{ B12 }$$

$$\begin{equation} \frac{\partial {\bf v}^{\prime }}{\partial e} = \frac{\sin {\mathcal {E}}}{n}{\bf a}^{\prime } +\frac{a\cos {\mathcal {E}}}{r^{\prime }}{\bf v}^{\prime } -\left({\begin{array}{@{}c@{}} 0\cr 1 \cr 0\end{array}}\right) \frac{a^2 ne\cos {\mathcal {E}}}{r^{\prime }\sqrt{1-e^2}}, \end{equation} \tag{ B13 }$$

where a' is the acceleration (not to be confused with a).

We then rotate these derivatives with ${\mathbb {R}}$:

$$\begin{equation} \frac{\partial {\bf r}}{\partial (e,a,t_{\rm peri})}= {\mathbb {R}}\frac{\partial {\bf r}^{\prime }}{\partial (e,a,t_{\rm peri})},\quad\!\! \frac{\partial {\bf v}}{\partial (e,a,t_{\rm peri})}= {\mathbb {R}}\frac{\partial {\bf v}^{\prime }}{\partial (e,a,t_{\rm peri})}, \end{equation} \tag{ B14 }$$

and evaluate ∂r/∂(Ω_node, i, ω_peri),

$$\begin{eqnarray} \frac{\partial {\bf r}}{\partial \Omega _{\rm node}} & = \frac{\partial {\mathbb {R}}}{\partial \Omega _{\rm node}} {\bf r}^{\prime } = \frac{\partial {\mathbb {R}}_z(\Omega _{\rm node})}{\partial \Omega _{\rm node}} {\mathbb {R}}_x(i) {\mathbb {R}}_z(\omega _{\rm peri}) {\bf r}^{\prime },\nonumber\\ \end{eqnarray} \tag{ B15 }$$

$$\begin{eqnarray} \frac{\partial {\bf r}}{\partial i} & = \frac{\partial {\mathbb {R}}}{\partial i} {\bf r}^{\prime } = {\mathbb {R}}_z(\Omega _{\rm node}) \frac{\partial {\mathbb {R}}_x(i)}{\partial i} {\mathbb {R}}_z(\omega _{\rm peri}) {\bf r}^{\prime }, \end{eqnarray} \tag{ B16 }$$

$$\begin{eqnarray} \frac{\partial {\bf r}}{\partial \omega _{\rm peri}} & = \frac{\partial {\mathbb {R}}}{\partial \omega _{\rm peri}} {\bf r}^{\prime } = {\mathbb {R}}_z(\Omega _{\rm node}) {\mathbb {R}}_x(i) \frac{\partial {\mathbb {R}}_z(\omega _{\rm peri})}{\partial \omega _{\rm peri}} {\bf r}^{\prime }, \nonumber\\ \end{eqnarray} \tag{ B17 }$$

and similarly for ∂v/∂(Ω_node, i, ω_peri).

From all derived derivatives we construct a 6 × 6 matrix and calculate its determinant, using LUP matrix decomposition (LU decomposition with partial pivoting) on the lower and upper triangular matrices. Then, the determinant is given by simple multiplication of all diagonal elements of the triangular matrices. The Jacobian (j_ph/kep) is given by the absolute value of the determinant.

It is useful to note that j_ph/kep goes to zero proportional to esin i, so for orbits close to circular or close to face-on, one must be careful about numerical problems when dividing by j_ph/kep.

In this appendix, we show how to translate phase-space parameters in the microlensing system coordinates (r, v) to orbital elements (e, a, t_peri, Ω_node, i, ω_peri).

The specific orbital energy (ε) and specific relative angular momentum (h) are conserved, so they can be calculated given the values of r and v at any time:

$$\begin{eqnarray} \varepsilon &= \frac{v^2}{2} - \frac{GM}{r} = -\frac{GM}{2a},\quad {\bf h}={\bf r} \times {\bf v}. \end{eqnarray} \tag{ B18 }$$

If we express r in the units of AU, v in the units of AU yr⁻¹, and M in M_☉, then G ≡ 4π². The semi-major axis can be calculated from Equation (B18), and the period can be derived from Kepler's third law:

$$\begin{eqnarray} a &= - \frac{GM}{2 \varepsilon }, \quad P = 2\pi \sqrt{\frac{a^3}{GM}}. \end{eqnarray} \tag{ B19 }$$

Next, we find the unit vectors (versors) of the system of coordinates relative to the orbit, with ${\hat{\bm {z}}}$ parallel to the angular momentum vector and ${\hat{\bm {x}}}$ parallel to the eccentricity vector,⁶⁴ e, which is a constant of motion and points in the direction of periastron,

$$\begin{equation} {\bf e} = \frac{{\bf v} \times {\bf h}}{GM} - \frac{{\bf r}}{r}. \end{equation} \tag{ B20 }$$

The eccentricity (e = ||e||) can also be calculated with formula 1 − e² = h²/(GMa).

The three versors are then given by

$$\begin{equation} {\hat{\bm {x}}}=\frac{{\bf e}}{e}, \quad {\hat{\bm {z}}} = \frac{{\bf h}}{h}, \quad {\hat{\bm {y}}}={\hat{\bm {z}}}\times {\hat{\bm {x}}}. \end{equation} \tag{ B21 }$$

From projections of the positional vector (r) onto the versors (${\hat{\bm {x}}}$ and ${\hat{\bm {y}}}$), we calculate the true anomaly (ν),

$$\begin{equation} \sin \nu = \frac{{\hat{\bm {y}}}\,{\bm \cdot}\, {\bf r}}{r}, \quad \cos \nu = \frac{{\hat{\bm {x}}}\,{\bm \cdot}\, {\bf r}}{r}, \end{equation} \tag{ B22 }$$

then the eccentric anomaly

$$\begin{equation} \cos {\mathcal {E}}= \frac{\cos \nu + e}{1+e \cos \nu }, \end{equation} \tag{ B23 }$$

with the phase ambiguity in ${\mathcal {E}}$ itself resolved by

$$\begin{equation} 0\le {\mathcal {E}}\le \pi \Longleftrightarrow 0\le \nu \le \pi,\; -\pi <{\mathcal {E}}<0 \Longleftrightarrow \pi <\nu <2\pi. \end{equation} \tag{ B24 }$$

This convention minimizes |t_peri − t_kep|.

The time of periastron can be derived from the inverted Equation (B3) written for the current time (t = t_{0, kep}):

$$\begin{equation} t_{\rm peri}= t_{0,{\rm kep}} - \frac{{\mathcal {E}}- e \sin {\mathcal {E}}}{n}. \end{equation} \tag{ B25 }$$

The coordinates of the three versors describing the orbital coordinate system (${\hat{\bm {x}}}$, ${\hat{\bm {y}}}$, ${\hat{\bm {z}}}$) give us the full rotation matrix ${\mathbb {R}}$, defined earlier by Equations (B8),

$$\begin{equation} {\mathbb {R}}= (\begin{array}{@{}ccc@{}} {\hat{\bm {x}}} & {\hat{\bm {y}}} & {\hat{\bm {z}}} \end{array}). \end{equation} \tag{ B26 }$$

We can use this matrix to find the angles Ω_node, i, and ω_peri from the set of nine Equations (B8) for every element of ${\mathbb {R}}$. The inclination, which is always 0 ⩽ i ⩽ π, is given by

$$\begin{equation} \cos i = {\mathbb {R}}_{33}. \end{equation} \tag{ B27 }$$

Then, the longitude of the ascending node can be calculated from

$$\begin{equation} \cos \Omega _{\rm node}= -\frac{{\mathbb {R}}_{23}}{\sin i}, \quad \sin \Omega _{\rm node}= \frac{{\mathbb {R}}_{13}}{\sin i}, \end{equation} \tag{ B28 }$$

and the argument of periapsis is

$$\begin{eqnarray} \cos \omega _{\rm peri}&= {\mathbb {R}}_{11} \cos \Omega _{\rm node}+ {\mathbb {R}}_{21} \sin \Omega _{\rm node}, \end{eqnarray} \tag{ B29 }$$

$$\begin{eqnarray} \sin \omega _{\rm peri}&= ({\mathbb {R}}_{21} \cos \Omega _{\rm node}- {\mathbb {R}}_{11} \sin \Omega _{\rm node})\cos i + {\mathbb {R}}_{31} \sin i. \nonumber\\ \end{eqnarray} \tag{ B30 }$$

(For the special case of face-on orbits (sin i = 0), we assume ω_peri ≡ 0 and calculate Ω_node from $\cos \Omega _{\rm node}= {\mathbb {R}}_{11}$ and $\sin \Omega _{\rm node}= {\mathbb {R}}_{21}$.)

The relative position of the two binary components at any time, r(t), can be calculated from Equation (B7), where ${\mathbb {R}}$ is a time-invariant matrix, and r'(t) is given by Equations (B2) and (B3).