OGLE-2016-BLG-0613LABb: A Microlensing Planet in a Binary System

Since the sharp spikes in the lensing light curve are characteristic features of caustic-crossing binary-lens events, we start the modeling of the light curve based on the assumption that the lens is composed of two masses, M₁ and M₂. For the simple case in which the relative lens-source motion is rectilinear, the light curve of a binary-lens event is described by seven geometric parameters plus two parameters representing the fluxes from the source, F_s, and blend, F_b, for each data set. The geometric lensing parameters include three of a single-lens event (${t}_{0},{u}_{0}$ and ${t}_{{\rm{E}}}$), another three parameters describing the binary lens ($s,q$, and α), and the ratio of the angular source radius ${\theta }_{* }$ to the angular Einstein radius ${\theta }_{{\rm{E}}}$, i.e., $\rho ={\theta }_{* }/{\theta }_{{\rm{E}}}$ (normalized source radius). The lengths of u₀ and s are normalized to ${\theta }_{{\rm{E}}}$.

In the preliminary binary-lens modeling, we first conduct a grid search over the parameter space of $s,q$, and α, while the remaining parameters (${t}_{0},{u}_{0},{t}_{{\rm{E}}}$, and ρ) are searched for using a downhill approach. We choose $s,q$, and α as the grid parameters because lensing magnifications vary sensitively to the small changes of these parameters. For the downhill approach, we use the Markov Chain Monte Carlo (MCMC) method. The initial values of the MCMC parameters are roughly guessed considering the characteristics of the light curve such as the duration, caustic-crossing times, etc. From the preliminary search, we find that (1) there exist multiple local minima that can describe the overall feature of the light curve but (2) none of these solutions can explain the short-term anomaly at $\mathrm{HJD}^{\prime} \sim 7493$.

It is a well-known fact that there can exist multiple degenerate solutions in binary lensing modeling due to deep symmetries in the lens equation (Griest & Safizadeh 1998; Dominik 1999b; An 2005). There is also at least one well-studied "accidental" degeneracy between binary light curves due to planetary-mass and roughly equal-mass binaries, which does not occur for any deep reason (Choi et al. 2012; Bozza et al. 2016). However, the general problem of "accidental" degeneracies is poorly studied. One general principle, however, is that the lower the quality and quantity of data, the wider is the range of possible binary geometries that may fit the data reasonably well. As seen from Figure 2, OGLE-2016-BLG-0613 was exceptionally faint for a planet-bearing microlensing event, never getting brighter than $I\sim 17.5$. Moreover, although it was overall densely covered, by chance, both the binary caustic entrance (${\mathrm{HJD}}^{{\prime} }\sim 7483$) and exit (${\mathrm{HJD}}^{{\prime} }\sim 7495$) fell in the gap between KMTC and KMTA coverage. Hence, there is a total of only one data point (from UKIRT) on these features. Binary events whose caustics lack such coverage are rarely, if ever, intensively studied, so that little is known about how they are impacted by "accidental" degeneracies. Therefore, one must pay particular attention to identifying all degenerate topologies.

Although several local solutions are found from the preliminary grid search, some local solutions might be still missed possibly due to a poor guess of the initial values of the MCMC parameters or some other reasons. We, therefore, check the existence of additional local solutions using two systematic approaches.

In the first approach, we conduct a series of additional grid searches in which we provide various combinations of MCMC parameters as initial values. For a single lensing event, the values of the MCMC parameters are well characterized by the peak time (for t₀), peak magnification (for u₀), and the duration of the event (for ${t}_{{\rm{E}}}$). For OGLE-2016-BLG-0613, however, it is difficult to estimate these values from the light curve and thus it might be that some local solutions have been missed if the given initial parameters were too far away from the correct ones. From these searches, we identify five local solutions. Among them, two local solutions were missed in the initial search mainly due to the large difference between the Einstein timescales given as an initial value and the recovered value from modeling.

In the second approach, we directly consider each of 21 different caustic topologies that are consistent with the overall morphology of the light curve. As seen in Figure 2, the caustic structures appear near the peak of the pre-caustic plus post-caustic light curve. The caustic structure must therefore be either a four-sided ("central") caustic or a six-sided ("resonant") caustic. Allowing for the symmetry of these caustic structures around the binary axis, there are $(4/2)\times (4-1)=6$ topologically distinct source paths for the central caustic and $(6/2)\times (6-1)=15$ for the resonant caustic. We seed each of these topologies with an arbitrary (s,q) geometry (that permits such a path) and set initial values of $({t}_{0},{u}_{0},{t}_{{\rm{E}}},\alpha )$ such that the caustic entrance and exit occur at approximately the correct times. We then allow all parameters to vary using ${\chi }^{2}$ minimization. Although the initial seed solutions generally provide extremely poor matches to the data, the derived local minima are always in rough accord with the data, showing that the approach is working. Nevertheless, only four of these solutions have ${\chi }^{2}$ within a few hundred of the global minimum. These yield the same five solutions found by the grid searches above, with one of the four topology solutions corresponding to a close-wide pair of grid-search solutions (see below).

In Figure 4, we present the locations of the local solutions in the ${\rm{\Delta }}{\chi }^{2}$ map of the log s–log q parameter space. The individual local solutions are marked by circles and labeled as "Sol A," "Sol B," "Sol C," "Sol D," and "Sol E," respectively. In Table 2, we present the lesing parameters of the individual local binary-lensing solutions. In Figure 5, we also present the geometry of the lens systems, which show the source trajectory with respect to the positions of the lens components and caustics. One finds that the lensing parameters of the local solutions span over wide ranges. For example, the range of the binary mass ratio is $0.03\lesssim q\lesssim 6.0$.²⁷ This implies that the observed light curve is serendipitously described by multiple local solutions with widely different lensing parameters.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Table 2.  Local Binary-lens Solutions

Parameters	"Sol A"	"Sol B"	"Sol C"	"Sol D"	"Sol E"
${\chi }^{2}$	9212.4	9185.7	9179.2	9183.6	9195.1
t0 (HJD')	7493.443 ± 0.055	7490.032 ± 0.053	7493.872 ± 0.054	7490.004 ± 0.081	7486.311 ± 0.109
u0	0.069 ± 0.001	0.046 ± 0.001	0.023 ± 0.001	0.039 ± 0.001	−0.006 ± 0.007
${t}_{{\rm{E}}}$ (days)	44.63r ± 0.39	52.87 ± 0.66	74.09 ± 0.20	53.39 ± 0.40	17.15 ± 0.20
s	1.011 ± 0.006	0.730 ± 0.006	1.393 ± 0.003	1.926 ± 0.008	1.228 ± 0.006
q	0.050 ± 0.002	0.382 ± 0.005	0.026 ± 0.001	1.002 ± 0.081	6.032 ± 0.002
α (rad)	3.157 ± 0.010	3.681 ± 0.009	2.915 ± 0.008	3.611 ± 0.012	−0.348 ± 0.016
ρ (10−3)	0.21 ± 0.07	0.21 ± 0.04	0.31 ± 0.04	0.31 ± 0.05	0.31 ± 0.02

Note. HJD' = HJD–2,450,000.

Download table as:  ASCIITypeset image

Although the overall shape of the observed light curve is described by multiple solutions, none of the solutions can explain the short-term anomaly at $\mathrm{HJD}^{\prime} \sim 7493$. The anomaly is unlikely to be caused by systematics in the data because the feature was covered commonly by the OGLE, KMTS, KMTA, and UKIRT data. The region between caustic spikes can deviate from a smooth U shape if the source within the caustic asymptotically approaches the caustic curve. In such a case, however, the deviation occurring during the caustic approach is smooth, while the observed short-term anomaly appears to be a small-scale caustic-crossing feature composed of both caustic entrance (at $\mathrm{HJD}^{\prime} \sim 7292.5$) and exit ($\mathrm{HJD}^{\prime} \sim 7293.1$) and the U-shape trough between them.

Binary lenses form closed curves, it is mathematically impossible for a binary lens to generate two successive caustic entrances without an intervening caustic exit, nor similarly two successive caustic exits without an intervening caustic entrance. In the present case, we have both two such successive entrances followed by two such successive exits. This suggests that one should consider a new interpretation of the event other than the binary-lens interpretation.

It is known that a third body of a lens can cause caustic curves to be self-intersected (Schneider et al. 1992; Petters et al. 2001; Danék & Heyrovský 2015; Luhn et al. 2016), and the light curve resulting from the source trajectory passing over the intersected part of the caustic can result in an additional caustic-crossing feature within the major caustic feature. We therefore conduct a triple-lens modeling of the observed light curve in order to check whether the short-term anomaly can be explained by a third body.

The three-body lens modeling is conducted in two steps. In the first step, we conduct a grid search for the parameters related to the third body (${s}_{3},{q}_{3}$, and ψ) by fixing the binary-lens parameters at the values obtained from the binary-lens modeling. Once approximate values of the third-body parameters are found, we then refine the three-body solution by allowing all parameters to vary. The first step is based on the assumption that the overall light curve is well described by a binary-lens model and the signal of the third body can be treated as a perturbation to the binary-lens curve. For OGLE-2016-BLG-0613, this assumption is valid because of the good binary-lens fit to the overall light curve and the short-term nature of the anomaly. The lensing parameters related to the third body include the separation s₃ and the mass ratio ${q}_{3}={M}_{3}/{M}_{1}$ between the third and the primary of the binary lens and the position angle of the third body measured from the binary axis, ψ. According to our definition of the position angle ψ, the third body is located at $({x}_{{\rm{L}},3},{y}_{{\rm{L}},3})\,=({x}_{{\rm{L}},1}+{s}_{3}\cos \psi ,{y}_{{\rm{L}},1}+{s}_{3}\sin \psi )$, where $({x}_{{\rm{L}},1},{y}_{{\rm{L}},1})$ represents the position of the primary of the binary lens, i.e., M₁.

From triple-lens modeling, we find that the short-term anomaly can be explained by introducing a low-mass third body to the binary-lens solutions "Sol A" through "Sol D." For the case of "Sol E," we find no triple-lens solution that can explain the short-term anomaly. For each of the solutions that can explain the short-term anomaly, we find a pair of triple-lens solutions resulting from the close/wide degeneracy of the third-body, i.e., s₃ versus (approximately) ${s}_{3}^{-1}$ (Han et al. 2013; Song et al. 2014). We designate the solutions with ${s}_{3}\gt 1$ and ${s}_{3}\lt 1$ as "wide" and "close" solutions, respectively.²⁸

In Table 3, we list the lensing parameters of the triple-lens solutions along with ${\chi }^{2}$ values. Also presented are the source and blend magnitudes, I_s and I_b, respectively.²⁹ From the mass ratios of the third body, one finds that $1.8\times {10}^{-3}\lesssim {q}_{3}\lesssim 6.4\times {10}^{-3}$, indicating that the third body is a planetary-mass object regardless of the models. From the mass ratios between the binary components, one finds that ${q}_{2}\lesssim 0.06$ for "Sol A" and "Sol C," while ${q}_{2}\gtrsim 0.35$ for "Sol B" and "Sol D." This indicates that the lens of "Sol A" and "Sol C" would be composed of a star, a brown dwarf (BD), and a planet, while the lens of "Sol B" and "Sol D" would consist of a stellar binary plus a planet.

Table 3.  Local Triple-lens Solutions

Parameters	"Sol A"		"Sol B"		"Sol C"		"Sol D"
	Wide	cClose	Wide	Close	Wide	Close	Wide	Close
${\chi }^{2}$	4869.9	4881.05	4799.2	4812.9	4789.2	4805.6	4802.4	4816.8
t0 (HJD')	7493.551 ± 0.038	7493.572 ± 0.043	7490.135 ± 0.063	7490.489 ± 0.120	7494.153 ± 0.048	7494.177 ± 0.017	7490.095 ± 0.075	7489.920 ± 0.039
u0	0.068 ± 0.001	0.072 ± 0.001	0.048 ± 0.001	0.047 ± 0.001	0.021 ± 0.001	0.022 ± 0.001	0.038 ± 0.002	0.038 ± 0.001
${t}_{{\rm{E}}}$ (days)	44.48 ± 0.71	42.13 ± 0.19	51.45 ± 0.41	44.94 ± 1.33	74.62 ± 1.69	71.90 ± 0.57	53.53 ± 0.43	52.01 ± 0.10
s2	1.009 ± 0.004	1.025 ± 0.004	0.743 ± 0.005	0.802 ± 0.015	1.396 ± 0.010	1.386 ± 0.004	1.941 ± 0.019	1.932 ± 0.003
q2	0.053 ± 0.001	0.055 ± 0.001	0.386 ± 0.010	0.359 ± 0.011	0.029 ± 0.002	0.027 ± 0.001	1.051 ± 0.052	1.114 ± 0.019
α (rad)	3.169 ± 0.009	3.172 ± 0.012	3.659 ± 0.005	3.663 ± 0.004	2.948 ± 0.010	2.954 ± 0.003	3.572 ± 0.014	3.587 ± 0.005
s3	1.111 ± 0.001	0.971 ± 0.001	1.064 ± 0.001	0.872 ± 0.002	1.168 ± 0.008	0.872 ± 0.005	0.833 ± 0.012	0.677 ± 0.004
q3 (10−3)	2.44 ± 0.07	2.07 ± 0.01	5.54 ± 0.13	4.43 ± 0.13	3.27 ± 0.24	3.24 ± 0.18	6.39 ± 0.49	1.77 ± 0.05
ψ (rad)	5.079 ± 0.010	5.079 ± 0.007	4.672 ± 0.011	4.774 ± 0.022	5.332 ± 0.016	5.350 ± 0.004	4.519 ± 0.011	1.852 ± 0.011
ρ (10−3)	0.35 ± 0.03	0.39 ± 0.03	0.30 ± 0.03	0.36 ± 0.03	0.22 ± 0.02	0.23 ± 0.02	0.32 ± 0.03	0.25 ± 0.02
${I}_{s,\mathrm{OGLE}}$ (mag)	21.84 ± 0.01	21.81 ± 0.01	22.25 ± 0.02	22.01 ± 0.04	23.00 ± 0.05	22.90 ± 0.02	21.96 ± 0.02	21.96 ± 0.01
Hs (mag)	18.55 ± 0.04	18.52 ± 0.04	18.96 ± 0.04	18.72 ± 0.06	19.71 ± 0.06	19.61 ± 0.04	18.67 ± 0.04	18.67 ± 0.04
${I}_{b,\mathrm{OGLE}}$ (mag)	19.58	19.59	19.55	19.56	19.50	19.51	19.57	19.57

Note. HJD' = HJD–2,450,000.

Download table as:  ASCIITypeset image

The model light curve of the triple-lens solution and the residual from the observed data are presented in Figure 6 for the "Sol B (wide)" model as a representative model. We note that ${\chi }^{2}$ differences among the models "Sol B," "Sol C," and "Sol D" are $\lesssim 30$ and thus the fits of "Sol C" and "Sol D" are similar to the fit of the presented model. For "Sol A," on the other hand, the fit is worse than the presented model by ${\rm{\Delta }}{\chi }^{2}\sim 80$. Figure 7 shows the caustic geometry of the individual triple-lens solutions. In all cases, it is found that the short-term anomaly is produced by the deformation of the caustic caused by the third body.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

We check whether the model further improves by additionally considering the parallax effect induced by the orbital motion of the Earth around the Sun (Gould 1992). We find that the microlens parallax ${\pi }_{{\rm{E}}}$ can be neither reliably measured nor meaningfully constrained. We note that the event was in the field of the space-based lensing survey using the Kepler space telescope (K2C9: Henderson et al. 2016). Since the Kepler telescope is in a heliocentric orbit, the space-based observations could have led to the measurement of the microlens parallax (Refsdal 1966; Gould 1992). The K2C9 campaign was planned to start at $\mathrm{HJD}^{\prime} \sim 7486$, when the event was in progress. However, the campaign could start only at $\mathrm{HJD}^{\prime} \sim 7501$ because of an emergency mode and thus the event was missed.

We characterize the source star by measuring the dereddened color and brightness, which are calibrated using the centroid of the giant clump (GC) in the color–magnitude diagram (CMD; Yoo et al. 2004). Although the event was observed in the V band by both the OGLE and KMTNet surveys, the V-band photometry quality is not good enough for a reliable measurement of the V-band baseline source flux due to the faintness of the source star combined with the high extinction toward the field. We therefore use the OGLE I-band data and UKIRT H-band data.

In Figure 8, we mark the location of the source star with respect to the GC centroid in the $(I-H)/I$ CMD of neighboring stars around the source. We measure the GC centroid ${(I-H,I)}_{\mathrm{GC}}=(3.58,16.57)\pm (0.01,0.07)$, where the errors are derived from the standard error of the mean of clump stars. As is virtually always the case, for each model, the source color is essentially identical ${(I-H)}_{s}=3.29\pm 0.04$ to much greater precision than the measurement error, which is dominated by the H-band photometric errors. Thus ${\rm{\Delta }}(I-H)=-0.29\pm 0.04$. Transformation to $(V-I)$ of both the value and error using the color–color relation of Bessell & Brett (1988) yields ${\rm{\Delta }}(V-I)=-0.21\pm 0.03$. Based on the well-defined dereddened color of GC centroid ${(V-I)}_{0,\mathrm{GC}}=1.06$ (Bensby et al. 2011), it is estimated that the dereddened color of the source is

$$\begin{eqnarray}&&{(V-I)}_{0,s}={(V-I)}_{0,\mathrm{GC}}+{\rm{\Delta }}(V-I)=0.85\pm 0.03.\end{eqnarray} \tag{ 2 }$$

This color corresponds to that of a late G-type main-sequence star with an absolute magnitude of ${M}_{I,s}\sim 4.9$. From the absolute magnitude of the GC centroid of ${M}_{I,\mathrm{GC}}=-0.12$ (Nataf et al. 2013) combined with its apparent brightness ${I}_{\mathrm{GC}}=16.57$ measured on the CMD, then, the apparent magnitude of the source star should be

$$\begin{eqnarray}&&{I}_{s}={I}_{\mathrm{GC}}-({M}_{I,\mathrm{GC}}-{M}_{I,s})\sim 21.6.\end{eqnarray} \tag{ 3 }$$

Zoom In Zoom Out Reset image size

It is found that the observed light curve can be explained by $(4\times 2=)8$ degenerate triple-lens solutions. In this subsection, we further investigate the individual solutions in order to check the feasibility of resolving the degeneracy among the solutions.

First, we exclude "Sol A" due to its relatively poor fit to the observed light curve compared to the other solutions. From the comparison of ${\chi }^{2}$ values presented in Table 3, it is found that "Sol A" provides a fit that is poorer than "Sol B," "Sol C," and "Sol D" by ${\rm{\Delta }}{\chi }^{2}=70.7$, 80.7, and 67.5, respectively. These ${\chi }^{2}$ differences are statistically significant enough to exclude "Sol A."

This leaves three pairs of solutions, one pair for each of "Sols B, C, D." "Sol C" (wide) is favored by ${\rm{\Delta }}{\chi }^{2}=10$ over any other solutions. However, the source position on the CMD is a priori substantially less likely than for "Sol B" and "Sol D." The source brightness determined from the source flux F_s, presented in Table 3, are ${I}_{s}=23.00\pm 0.05$ and 22.90 ± 0.02 for the wide and close solutions, respectively. These are $\gtrsim 1.3$ magnitude fainter than I_s = 21.6, i.e., Equation (3), estimated based on the relative source position with respect to the GC centroid in the CMD. See Figure 8. To explain such a faint source would require either that the source is very distant (at, e.g., ${D}_{{\rm{S}}}\sim 15\,$ kpc), or that it is intrinsically dim for its color due to low metallicity (e.g., [Fe/H]$\sim -1.3$). Either of these is possible, although with low probability. For example, Bensby et al. (2017a) find a total of three stars with [Fe/H]$\lt -1.2$ out of their sample of 90 microlensed bulge dwarfs and subgiants. Similarly, Ness et al. (2013) found about 2% of stars with $[\mathrm{Fe}/{\rm{H}}]\lt -1.2$ in a much larger sample, albeit substantially farther from the Galactic plane than typical microlensing events. Compare Figures 1 and 10 of Bensby et al. (2017b) with Figures 1 and 6 of Ness et al. (2013), respectively. The Bensby et al. (2017a) sample is appropriate for comparison because it faces qualitatively similar selection biases to microlensing surveys that lead to planet detection. On the other hand, Bensby et al. (2017a) found no clear evidence for microlensed sources that lay substantially behind the Galactic bulge.

In brief, "Sol C" is mildly favored by ${\chi }^{2}$ but requires a somewhat unlikely source. If its ${\rm{\Delta }}{\chi }^{2}=10$ advantage could be interpreted at face value according to Gaussian statistics, this solution would be mildly preferred by ${\rm{\Delta }}\,\mathrm{ln}\,L\,=\mathrm{ln}(3/90)+10/2=1.6$. However, it is well-known that microlensing light curves are affected by subtle systematics at the ${\rm{\Delta }}{\chi }^{2}\,=$ few level, and so cannot be judged according to naive Gaussian statistics. Therefore, we consider that "Sol C" is viable but somewhat disfavored.

On the other hand, "Sol B" and "Sol D" not only provide good fits to the observed light curve but also meet the source brightness constraint. The apparent source magnitudes estimated from F_s values of the models are ${I}_{s}\,=22.25\pm 0.02$/22.01 ± 0.04 for the wide/close solutions of "Sol B" and ${I}_{s}=21.96\pm 0.02$/21.96 ± 0.01 for the wide/close solutions of "Sol D." These are in accordance with ${I}_{s}\sim 21.6$ estimated from the CMD. We find that the lensing parameters of "Sol B" and "Sol D" models are in the relation of

$$\begin{eqnarray}&&{s}_{B}\sim \displaystyle \frac{1+{q}_{B}}{1-{q}_{B}\sqrt{1+{q}_{D}}}{s}_{D}^{-1},\end{eqnarray} \tag{ 4 }$$

where $({s}_{B},{q}_{B})$ represent the binary separation and mass ratio of "Sol B" and $({s}_{D},{q}_{D})$ represent those of "Sol D." This indicates that "Sol B" and "Sol D" are the pair of solutions resulting from the well-known close/wide binary-lens degeneracy, which is rooted in the symmetry of the lens equation and thus can be severe (Dominik 1999b; An 2005). For OGLE-2016-BLG-0613, the degeneracy is quite severe with ${\rm{\Delta }}{\chi }^{2}\lesssim 5$. Note in particular from Figure 5 that the topologies of "Sol B" and "Sol D" are essentially identical and are also distinct from the other three topologies shown.

The angular Einstein radius is estimated from the combination of the normalized source radius ρ and the angular source radius ${\theta }_{* }$, i.e., ${\theta }_{{\rm{E}}}={\theta }_{* }/\rho$. The normalized source radius is measured by analyzing the caustic-crossing part of the light curve, where the lensing magnifications are affected by finite-source effects. Although three of the four caustic crossings were covered by either zero or one point (and so yield essentially no information about ρ), the planetary-caustic entrance was very well covered by KMTS data. See the lower middle panel of Figure 6.

The angular source radius is estimated based on the dereddened color and brightness. The dereddeded V-band brightness is estimated by

$$\begin{eqnarray}&&{V}_{0,s}={I}_{0,s}+{(V-I)}_{0,s},\end{eqnarray} \tag{ 5 }$$

where ${I}_{0,s}={I}_{\mathrm{GC},0}-({I}_{\mathrm{GC}}-{I}_{s}),{I}_{\mathrm{GC},0}=14.38$ (Nataf et al. 2013), ${I}_{\mathrm{GC}}\,=\,16.57\,\pm \,0.07$, ${(V-I)}_{0,s}=0.85\,\pm \,0.03$ (Equation (2)), and I_s is given for each solution in Table 3. We convert ${(V-I)}_{0}$ into ${(V-K)}_{0}$ using the VI/VK relation of Bessell & Brett (1988). Then, the angular source radius ${\theta }_{* }$ is determined from the relation between ${\theta }_{* }$ and ${(V,V-K)}_{0}$ provided by Kervella et al. (2004).

In Table 4, we present the angular Einstein radii for the viable models "Sol B," "Sol C," and "Sol D." Also presented are the relative lens-source proper motion determined by

$$\begin{eqnarray}&&\mu =\displaystyle \frac{{\theta }_{{\rm{E}}}}{{t}_{{\rm{E}}}},\qquad {\theta }_{{\rm{E}}}\equiv \sqrt{\kappa M{\pi }_{\mathrm{rel}}},\end{eqnarray} \tag{ 6 }$$

where ${\pi }_{\mathrm{rel}}=\mathrm{au}({D}_{L}^{-1}-{D}_{S}^{-1})$ is the lens-source relative parallax and $\kappa \equiv 4G/({c}^{2}\mathrm{au})\simeq 8.14\ \mathrm{mas}\ {M}_{\odot }^{-1}$. The inferred angular Einstein radii for the ensemble of solutions are in the range of $1.1\lesssim {\theta }_{{\rm{E}}}/\mathrm{mas}\lesssim 1.6$. These large values of ${\theta }_{{\rm{E}}}$ virtually rule out bulge lenses and so directly imply that the lens is in the Galactic disk. That is, if we hypothesized that the lens were in the bulge, then (since the bulge is a relatively old population), we could infer $M\lesssim 1.3\,{M}_{\odot }$. Hence, for ${\theta }_{{\rm{E}}}\gtrsim 1.1\,\mathrm{mas}$, Equation (6) implies ${\pi }_{\mathrm{rel}}\gtrsim 0.11\,\mathrm{mas}$, which would contradict the hypothesis that the lens was in the bulge.

Table 4.  Einstein Ring Radius and Proper Motion

Parameters	"Sol B"		"Sol C"		"Sol D"
	Wide	Close	Wide	Close	Wide	Close
${\theta }_{{\rm{E}}}$ (mas)	1.20 ± 0.24	1.12 ± 0.22	1.15 ± 0.22	1.15 ± 0.22	1.28 ± 0.25	1.63 ± 0.32
μ (mas yr−1)	8.50 ± 1.71	9.06 ± 1.75	5.64 ± 1.11	5.86 ± 1.14	8.70 ± 1.72	11.5 ± 2.20

Download table as:  ASCIITypeset image

For the unique determinations of the lens mass M and distance ${D}_{{\rm{L}}}$, one needs to measure both the angular Einstein radius ${\theta }_{{\rm{E}}}$ and the microlens parallax ${\pi }_{{\rm{E}}}$:

$$\begin{eqnarray}&&M=\displaystyle \frac{{\theta }_{{\rm{E}}}}{\kappa {\pi }_{{\rm{E}}}};\qquad {D}_{{\rm{L}}}=\displaystyle \frac{\mathrm{au}}{{\pi }_{{\rm{E}}}{\theta }_{{\rm{E}}}+{\pi }_{{\rm{S}}}},\end{eqnarray} \tag{ 7 }$$

where ${\pi }_{{\rm{S}}}=\mathrm{au}/{D}_{S}$ is the source parallax. For OGLE-2016-BLOG-0613, ${\theta }_{{\rm{E}}}$ is measured but ${\pi }_{{\rm{E}}}$ is not measured and thus the values of M and ${D}_{{\rm{L}}}$ cannot be uniquely determined. However, one can still constrain the physical lens parameters based on the measured values of the event timescale ${t}_{{\rm{E}}}$ and the angular Einstein radius ${\theta }_{{\rm{E}}}$.

In order to estimate the mass and distance to the lens, we conduct a Bayesian analysis of the event based on the mass function combined with the models of the physical and dynamical distributions of objects in the Galaxy. We use the initial mass function of Chabrier (2003a) for the mass function of Galactic bulge objects, while we use the present day mass function of Chabrier (2003b) for disk objects. In the mass function, we do not include stellar remnants, i.e., white dwarfs, neutron stars, and black holes, because it would be difficult for planets to survive the AGB/planetary-nebular phase of stellar evolution and no planet belonging to a remnant host is so far known, e.g., Kilic et al. (2009). For the matter distribution, we adopt the Galactic model of Han & Gould (2003). In this model, the matter density distribution is constructed based on a double-exponential disk and a triaxial bulge. We use the dynamical model of Han & Gould (1995) to construct the velocity distribution. In this model, the disk velocity distribution is assumed to be Gaussian about the rotation velocity of the disk and the bulge velocity distribution is modeled to be a triaxial Gaussian with velocity components deduced from the flattening of the bulge via the tensor virial theorem. Based on these models, we generate a large number of artificial lensing events by conducting a Monte Carlo simulation. We then estimate the ranges of M and ${D}_{{\rm{L}}}$ corresponding to the measured event timescale and the angular Einstein radius.

In Table 5, we list the physical parameters of the lens system estimated from the the Bayesian analysis. In Figure 9, we also present the distributions of the primary mass M₁ and the distance to the lens ${D}_{{\rm{L}}}$ obtained from the Bayesian analysis for the individual solutions. The values denoted by M_i represent the masses of the individual lens components, and ${a}_{\perp ,1-2}$ and ${a}_{\perp ,1-3}$ denote the projected ${M}_{1}-{M}_{2}$ and ${M}_{1}-{M}_{3}$ separations, respectively. We note that the unit of M₁ and M₂ is the solar mass, ${M}_{\odot }$, while the unit of M₃ is the Jupiter mass, M_J. The presented physical parameters are the median values of the corresponding distributions and the uncertainties are estimated as the standard deviations of the distributions.

Zoom In Zoom Out Reset image size

Table 5.  Physical Lens Parameters

Parameters	"Sol B"		"Sol C"		"Sol D"
	Wide	Close	Wide	Close	Wide	Close
${D}_{{\rm{L}}}$ (kpc)	${3.41}_{-1.42}^{+1.34}$	${3.56}_{-1.51}^{+1.30}$	${3.19}_{-1.27}^{+1.40}$	${3.21}_{-1.20}^{+1.35}$	${3.84}_{-1.63}^{+1.44}$	${2.93}_{-1.44}^{+1.59}$
M1 (${M}_{\odot }$)	${0.72}_{-0.42}^{+0.55}$	${0.65}_{-0.37}^{+0.50}$	${0.80}_{-0.43}^{+0.65}$	${0.81}_{-0.44}^{+0.64}$	${0.66}_{-0.38}^{+0.58}$	${0.63}_{-0.41}^{+0.72}$
M2 (${M}_{\odot }$)	${0.28}_{-0.16}^{+0.21}$	${0.23}_{-0.13}^{+0.18}$	${0.023}_{-0.012}^{+0.019}$	${0.023}_{-0.012}^{+0.017}$	${0.69}_{-0.40}^{+0.61}$	${0.70}_{-0.46}^{+0.80}$
M3 (MJ)	${4.18}_{-2.43}^{+3.19}$	${3.01}_{-1.71}^{+2.32}$	${2.74}_{-1.51}^{+2.23}$	${2.75}_{-1.49}^{+2.17}$	${4.42}_{-2.54}^{+3.88}$	${1.17}_{-0.76}^{+1.33}$
${a}_{\perp ,1-2}$ (au)	${2.96}_{-1.23}^{+1.16}$	${3.02}_{-1.28}^{+1.10}$	${5.00}_{-1.99}^{+2.19}$	${5.00}_{-1.87}^{+2.10}$	${9.09}_{-3.86}^{+3.41}$	${8.66}_{-0.38}^{+0.42}$
${a}_{\perp ,1-3}$ (au)	${6.40}_{-2.63}^{+2.51}$	${3.28}_{-1.39}^{+1.20}$	${4.19}_{-1.67}^{+1.84}$	${3.15}_{-1.18}^{+1.32}$	${3.90}_{-1.66}^{+1.46}$	${3.04}_{-0.13}^{+0.15}$

Download table as:  ASCIITypeset image

All three pairs of solutions shown in Table 5 are comprised of a super-Jupiter planet in a binary-star system, whose primary is a K dwarf. They are all located in the Galactic disk, about halfway toward the bulge. The major difference among these solutions is that for "Sol B" and "Sol D," the components of the binary are of equal (B) or comparable (D) mass, whereas for "Sol C," the second component of the binary is a low-mass brown dwarf. We note that in all cases, the projected separations of the secondary and the planet are roughly comparable, so that if the system lay in the plane of the sky it would be unstable. This is a result of a generic selection bias of microlensing. Binary lenses, and in particular planets, are most easily discovered if they lie separated in projection by roughly one Einstein radius. See, e.g., Figure 7 of Mróz et al. (2017). This bias affects multilens systems by the square. OGLE-2012-BLG-0026 (Han et al. 2013) provides an excellent example of such bias. Hence, stable, hierarchical systems are preferentially seen at an angle such that the projected separations are comparable.