MOA-2016-BLG-227Lb: A Massive Planet Characterized by Combining Light-curve Analysis and Keck AO Imaging

Gravitational microlensing is a powerful method for detecting extrasolar planets (Mao & Paczynski 1991; Gould & Loeb 1992; Gaudi 2012). Compared to other detection techniques, microlensing, which is complementary to other methods, is sensitive to low-mass planets (Bennett & Rhie 1996) orbiting beyond the snow line around relatively faint host stars like M dwarfs or brown dwarfs (Bennett et al. 2008; Sumi et al. 2016).

A difficulty with the microlensing method is the determination of the mass of a lens M_L and the distance to the lens system D_L. If we have an estimate for the angular Einstein radius θ_E and the microlens parallax π_E, the mass is directly determined by

$$\begin{eqnarray}&&{M}_{L}=\displaystyle \frac{{\theta }_{{\rm{E}}}}{\kappa {\pi }_{{\rm{E}}}},\end{eqnarray} \tag{ 1 }$$

where κ = 8.144 mas ${M}_{\odot }^{-1}$ (Gould 1992; Gaudi et al. 2008; Muraki et al. 2011). When the source distance, D_S ∼ 8 kpc, is known, the distance to the lens is given by

$$\begin{eqnarray}&&{D}_{L}=\displaystyle \frac{\mathrm{AU}}{{\pi }_{E}{\theta }_{E}+{\pi }_{S}},\end{eqnarray} \tag{ 2 }$$

where π_S ≡ AU/D_S. However, the microlens parallax can be observed for a fraction of planetary events, while the angular Einstein radius is observed for most planetary events.

One strategy to estimate M_L and D_L for events in which the microlens parallax cannot be detected is to use a Bayesian analysis based on probability distributions from a standard Galactic model (e.g., Beaulieu et al. 2006; Bennett et al. 2014; Koshimoto et al. 2014; Shvartzvald et al. 2014). However, such an analysis must necessarily make an assumption about the probability that stars of a given mass and distance will host a planet. The most common assumption is that all stellar microlens stars are equally likely to host a planet with the properties of the microlens planet in question. It may be that the probability of hosting a planet of the measured mass ratio and separation depends on the host mass or the distance from the Galactic center. But, without mass and distance measurements, these quantities are determined by our Bayesian prior assumptions. As a case in point, Bennett et al. (2014) analyzed MOA-2011-BLG-262 and found that a planetary mass host orbited by an Earth-mass "moon" model had almost the same likelihood as a star+planet model. But, since we have no precedent for such a rogue planet+moon system, they selected the more conventional star+planet system as the favored model. Also, the first discovered microlensing planet, OGLE-2003-BLG-235Lb, was at first thought to be a giant planet orbiting an M dwarf with a mass of M_* ∼ 0.36 M_⊙ from Bayesian analysis (Bond et al. 2004). Such a system is predicted to be rare according to the core accretion theory of planet formation (Laughlin et al. 2004; Kennedy & Kenyon 2008). Follow-up HST images revealed a more massive host star with mass of ${M}_{* }={0.63}_{-0.09}^{+0.07}{M}_{\odot }$ by detecting excess flux in multiple passbands (Bennett et al. 2006).

When we measure the lens flux with high angular resolution HST or adaptive optics (AO) images (e.g., Bennett et al. 2006, 2007, 2015; Batista et al. 2015), we can then calculate the lens mass M_L using a mass–luminosity relation combined with the mass–distance relation derived from the θ_E measurement. High angular resolution images are needed because microlensed source stars are generally located in dense Galactic bulge fields where there are usually multiple bright main sequence stars per ground-based seeing disk.

Because the size of the angular Einstein radius is ≲1 mas, and the lens–source relative proper motion is typically μ_rel ∼ 6 mas yr⁻¹, it is possible that the lens and source stars will remain unresolved even in high angular resolution images taken within a few years of the microlensing event. In such cases, there will be excess flux above that contributed by the source star, and this excess flux must include the lens star flux. Some studies (Batista et al. 2014; Fukui et al. 2015; Koshimoto et al. 2017), which detected an excess flux, have assumed that this excess flux is dominated by the lens flux, and they have derived the lens mass under this assumption.

With this method, it might seem that no assumptions are required regarding the probability of the microlens stars hosting planets, and there would be no biases due to any inadequacies of the Galactic model used. However, Bhattacharya et al. (2017) use HST imaging to show that the excess flux at the position of the MOA-2008-BLG-310 source is not due to the lens star, and N. Koshimoto et al. (2017, in preparation) have developed a Bayesian method to study the possibility of excess flux from stars other than the lens star. Possibilities include unrelated stars and companions to the source and lens stars. They find that it can be difficult to exclude all these contamination scenarios, especially for events with small angular Einstein radii. In those cases where we cannot exclude the contamination scenarios, we can again use a Bayesian analysis similar to the one described above to estimate the probability distribution of the lens properties. This means that we need to assume prior distributions for stellar binary systems and the stellar luminosity function even when we have detected excess flux in high angular resolution images. In cases where the lens properties are confirmed by a measurement of the lens–source relative proper motion (Batista et al. 2015; Bennett et al. 2015) or microlensing parallax measurements (Gaudi et al. 2008; Beaulieu et al. 2016; Bennett et al. 2016), this contamination can be ruled out. Attempts at lens–source relative proper motion measurements can also confirm contamination (Bhattacharya et al. 2017) in cases where the measured proper motion of the star responsible for the excess flux does not match the microlensing light-curve prediction.

In this paper, we report the discovery of the planetary microlensing event MOA-2016-BLG-227. Observations and data reduction are described in Sections 2 and 3. Our modeling results are presented in Section 4. In Section 5, we model the foreground extinction by comparing observed color–magnitude diagrams (CMDs) to different extinction laws, and compare the results from the different extinction laws. Then, we use the favored extinction law to determine the angular Einstein radius, θ_E. In Section 6, we describe our Keck AO observations and photometry, and we determine the excess flux at the position of the source. In Section 7, we describe our Bayesian method to determine the probability that this excess flux is due to lens star and various combinations of other "contaminating" stars. The posterior probabilities for this MOA-2016-BLG-227 planetary microlensing event are presented, and we consider the effect of different planet-hosting probability priors. Finally, we discuss and conclude the results of our work in Section 8.

The Microlensing Observations in Astrophysics (MOA; Bond et al. 2001; Sumi et al. 2003) group conducts a high cadence survey toward the Galactic bulge using the 2.2 deg² FOV MOA-cam3 (Sako et al. 2008) CCD camera mounted on the 1.8 m MOA-II telescope at the University of Canterbury Mt. John Observatory in New Zealand. The MOA group alerts about 600 microlensing events per year. Most observations are conducted in a customized MOA-Red filter, which is similar to the sum of the standard Cousins R-and I-band filters. Observations with the MOA V filter (Bessell V-band) are taken once every clear night in each MOA field.

The microlensing event MOA-2016-BLG-227 was discovered and announced by the MOA alert system (Bond et al. 2001) at (R.A., decl.)_J2000 = (18:05:53.70, −27:42:51.43) and (l, b) = (3303, −3240) on 5 May 2016 (HJD' ≡ HJD–2450000 ∼ 7514). This event occurred during the microlensing Campaign 9 of the K2 Mission (K2C9; Henderson et al. 2016), and it was located close to (but not in) the area of sky that was surveyed for the K2C9. This part of the K2 field that was downloaded at 30 minute intervals is known as the "superstamp." Because this event was so close to the superstamp, several other groups conducting observing campaigns coordinated with the K2C9 observations also observed this event.

The Wise group used the Jay Baum Rich telescope, a Centurion 28 inch telescope (C28) at the Wise Observatory in Israel, which is equipped with a 1 deg² camera. The group monitored the K2C9 superstamp during the campaign with six survey fields that were observed three to five times per night with the Astrodon Exo-Planet BB (blue-blocking) filter. Although the MOA-2016-BLG-227 target was just outside the K2C9 superstamp, it was still within the Wise survey footprint.

The event was also observed with the wide-field near-infrared (NIR) camera (WFCAM) on the UKIRT 3.8 m telescope on Maunakea, Hawaii, as part of an NIR microlensing survey in conjunction with the K2C9 (Shvartzvald et al. 2017) survey. The UKIRT survey covered 6 deg², including the entire K2C9 superstamp and extending almost to the Galactic plane, with a cadence of two to three observations per night. Observations were taken in the H-band, with each epoch composed of sixteen 5 s co-added dithered exposures (two co-adds, two jitter points, and 2 × 2 microsteps).

The Canada–France–Hawaii Telescope (CFHT), also on Maunakea, serendipitously observed the event during the CFHT-K2C9 Microlensing Survey. The CFHT operated a multicolor survey of the K2C9 superstamp using the Megacam Instrument (Boulade et al. 2003). The CFHT observations for the event were conducted through the g-, r-, and i-band filters.

The VLT Survey Telescope (VST) is a 2.61 m telescope installed at ESO's Paranal Observatory, and it carried out K2C9 observations as a 99 hr filler program (Arnaboldi et al. 1998; Kuijken et al. 2002). Observations for such a filler program could only be carried out whenever the seeing was worse than 1 arcsec or conditions were non-photometric. The main objective of the microlensing program was to monitor the K2C9 superstamp in an automatized mode to improve the event coverage and to secure color information in the SDSS r and Johnson V passbands. Due to weather conditions, Johnson V images were only taken in the second half of the K2C9 survey, and therefore MOA-2016-BLG-227 is only covered by SDSS r. The exact pointing strategy was adjusted to cover the superstamp with six pointings and to contain as many microlensing events from earlier seasons as possible. In addition, a two-point dither was obtained to reduce the impact of bad pixels and detector gaps. Consequently, some events, like MOA-2016-BLG-227, received more coverage and have been observed with different CCDs.

Figure 1 shows the observed MOA-2016-BLG-227 light curve. MOA announced the detection of a light-curve anomaly for this event on 2016 May 9 (HJD' = HJD −2450000 ∼ 7518), and identified the anomaly as a planetary signal 4.5 hr after the anomaly alert. Although MOA detected a strong planetary caustic exit, the observing conditions were poor at the MOA observing site both immediately before and after this strong light curve feature. Fortunately, the additional observations from the Wise, UKIRT, CFHT, and VST telescopes covered the other important features of the light curve.

Zoom In Zoom Out Reset image size

Photometry of the MOA, Wise, and UKIRT data was conducted using the offline difference image analysis pipeline of Bond et al. (2017) in which stellar images are measured using an analytical PSF model of the form used in the DoPHOT photometry code (Schechter et al. 1993).

Differential flux light curves of the CFHT data were produced from Elixir-calibrated images⁴⁹ using a custom difference imaging analysis pipeline based on ISIS version 2.2 (Alard & Lupton 1998; Alard 2000) and utilizing an improved interpolation routine⁵⁰ (Bertin & Arnouts 1996; Siverd et al. 2012). Further details of the CFHT data reduction will be presented in a future paper.

Since there is no public VST instrument pipeline, calibration images from ESO's archive were used and combined. Restrictive bad pixel masks were extracted to prevent the inclusion of flat-field pixels with >1% nightly variation or with >10% deviation from the average. The calibrated images were reduced with the difference imaging package DanDIA (Bramich 2008), which uses a numerical kernel for difference imaging and the routines from the RoboNet pipeline for photometry (Tsapras et al. 2009).

It is known that error bars estimated by crowded field photometry codes can be under- or overestimated depending on the specific details of the event. The error bars provided by the photometry codes are sufficient to find the best-fit models, but they do not allow a proper determination of the microlensing light-curve model parameter uncertainties. Therefore, we empirically normalize the error bars for each data set. We used the formula presented in Yee et al. (2012) for normalization, ${\sigma }_{i}^{\prime }=k\sqrt{{\sigma }_{i}^{2}+{e}_{\min }^{2}},$ where σ_i is the original error of the ith data point in magnitudes, and the parameters for normalization are k and e_min. The parameters k and e_min are adjusted so that the cumulative χ² distribution as a function of the number of data points sorted by each magnification of the preliminary best-fit model is a straight line of slope 1.

The data set used for our analysis and the obtained normalization parameters are summarized in Table 1.

The modeling of a binary lens event requires the following parameters: the time of the closest approach of the source to the lens center of mass, t₀; the impact parameter, u₀, of the source trajectory with respect to the center of mass of the lens system; the Einstein radius crossing time t_E = θ_E/μ_rel; the lens mass ratio, q ≡ M_p/M_host; the separation of the lens masses, s; and the angle between the trajectory and the binary lens axis, α, and the source size ρ ≡ θ_*/θ_E. The parameters u₀, s, and ρ are given in units of the Einstein radius, and M_host and M_p are the masses of the host star and its planetary companion. With these seven parameters, we can calculate the magnification as a function of time A(t). In the crowded stellar fields where most microlensing events are found, most source stars are blended with one or more other stars, so that we cannot determine the source star brightness directly from images where the source is not magnified. Therefore, we add another set of linear parameters for each data set, the source and blend fluxes, f_S and f_b, which are related to the observed flux by F(t) = f_S A(t) + f_b.

When we include the finite source effect, we must consider limb-darkening effects. We adopt a linear limb-darkening law with one parameter, u_λ, for each data set. From the intrinsic source color, discussed in Section 5.1, we choose the atmospheric parameters for stars with similar intrinsic color from Bensby et al. (2013). This yields an effective temperature of ${T}_{\mathrm{eff}}\sim 5500\,{\rm{K}}$, a surface gravity of $\mathrm{log}\,[g/(\mathrm{cm}\,{{\rm{s}}}^{-2})]=4.0$, a metallicity of [M/H] = 0.0, and a microturbulence velocity of ξ = 1.0 km s⁻¹. We select the limb-darkening coefficients from the ATLAS model by Claret & Bloemen (2011) using these atmospheric parameters. We have ${u}_{\mathrm{MOA} \mbox{-} \mathrm{Red}}=0.5585$ for MOA-Red, u_V = 0.6822 for MOA-V, u_R = 0.6015 for Wise Astrodon, u_H = 0.3170 for UKIRT H, u_i = 0.5360 for CFHT i, u_r = 0.6257 for CFHT r, VST-71 r, and VST-95 r, and u_g = 0.7565 for CFHT g. We used the mean of the u_I and u_R values for the limb-darkening coefficients for the MOA-Red passband. Here we adopted the R-band limb-darkening coefficient for the Wise Astrodon data. As the Wise Astrodon filter is non-standard, our choice is not perfect. However, we note that even if we adopt u = 0 for the limb-darkening coefficient used with the Wise data, our χ² value changes by only 1.5. The limb-darkening coefficients are also listed in Table 1.

To find the best-fit model, we conduct a global grid search using the method of Sumi et al. (2016), where we fit the light curves using the Metropolis algorithm (Metropolis et al. 1953), with magnification calculations from the image-centered ray-shooting method (Bennett & Rhie 1996; Bennett 2010). From this, we find a unique model in which the source crosses the resonant caustic. We show the model light curve in Figure 1, the caustic and the source trajectory in Figure 2, and the best-fit model parameters in Table 2 along with the parameter error bars, which are calculated with a Markov Chain Monte Carlo (MCMC; Verde et al. 2003).

Zoom In Zoom Out Reset image size

Table 2.  The Parameters for the Best-fit Binary Lens Model

Parameter	Unit	Value
t0	HJD–2450000	${7517.5078}_{-0.006}^{+0.007}$
tE	days	${17.03}_{-0.20}^{+0.08}$
u0	10−2	$-{8.33}_{-0.16}^{+0.08}$
q	10−3	${9.28}_{-0.11}^{+0.20}$
s		${0.9312}_{-0.0009}^{+0.0004}$
α	rad	${2.509}_{-0.004}^{+0.003}$
ρ	10−3	${3.01}_{-0.05}^{+0.09}$
χ2		2538.9
dof		2538

Notes. Superscripts and subscript indicates the the 84th and 16th percentile from the best-fit values, respectively.

Download table as:  ASCIITypeset image

We also model the light curve including the microlensing parallax effect due to the Earth's orbital motion (Gould 1992; Alcock et al. 1995), although this event is unlikely to reveal a significant microlensing parallax signal because of its relatively short timescale. We find that the inclusion of the parallax effect improves the fit by Δχ² ∼ 14. However, the parts of the light curve that contribute to this decrease in χ² have a scatter similar to the variability of the MOA baseline data, and the best-fit microlensing parallax parameter is abnormally large, ${\pi }_{E}={1.3}_{-0.3}^{+2.1}$, yielding a very small lens mass of M_L ∼ 0.02 M_⊙. Therefore, we conclude that the improvement of the fit by the parallax effect is due to systematic errors in the MOA baseline data.

Because we have measured the finite source size, ρ, to a precision of ∼2%, the determination of the angular source star radius θ_* will yield the angular Einstein radius θ_E = θ_*/ρ. This, in turn, provides the mass–distance relation (Bennett 2008; Gaudi 2012):

$$\begin{eqnarray}{M}_{L} & = & \displaystyle \frac{{c}^{2}}{4G}{\theta }_{E}^{2}\displaystyle \frac{{D}_{S}{D}_{L}}{{D}_{S}-{D}_{L}}=0.9823\,{M}_{\odot }{\left(\displaystyle \frac{{\theta }_{E}}{1\mathrm{mas}}\right)}^{2}\\ & & \times \left(\displaystyle \frac{x}{1-x}\right)\left(\displaystyle \frac{{D}_{S}}{8\,\mathrm{kpc}}\right),\end{eqnarray} \tag{ 3 }$$

where x = D_L/D_S. We can empirically derive θ_* from the intrinsic source magnitude and color (Kervella et al. 2004; Boyajian et al. 2014).

5.1. Calibration

Our fist step is to calibrate the source magnitude to a standard photometric system. We cross-referenced stars in the event field between our DoPHOT photometry catalog of stars in the MOA image and the OGLE-III catalog (Szymański et al. 2011) to convert MOA-Red and MOA V into standard magnitudes. Following the procedure presented in Bond et al. (2017), we find the relations

$$\begin{eqnarray}{I}_{\mathrm{OGLE}\mbox{--}\mathrm{III}}-{R}_{\mathrm{MOA}} & = & (28.186\pm 0.006)\\ & & -(0.247\pm 0.005){(V-R)}_{\mathrm{MOA}},\end{eqnarray} \tag{ 4 }$$

$$\begin{eqnarray}{V}_{\mathrm{OGLE}\mbox{--}\mathrm{III}}-{V}_{\mathrm{MOA}} & = & (28.391\pm 0.004)\\ & & -(0.123\pm 0.004){(V-R)}_{\mathrm{MOA}}.\end{eqnarray} \tag{ 5 }$$

Using these calibration formulae and the result of light-curve modeling, we obtain the source star magnitude I_S = 19.536 ± 0.019 and the color (V − I)_S = 1.60 ± 0.03.

We follow a similar procedure to cross-reference stars in our DoPHOT photometry catalog of stars in the UKIRT images to stars in the VVV (Minniti et al. 2010) catalog, which is calibrated to the Two Micron All Sky Survey (2MASS) photometric system (Carpenter 2001), thereby obtaining the relationship between these photometric systems. We use this same VVV catalog to plot CMDs in the next section and for the analysis of the Keck images in Section 6. Using the UKIRT H-band source magnitude obtained from the light-curve model and the calibration relation, we find H_S = 17.806 ± 0.017. We also measure the colors of the source star: (V − H)_S = 3.33 ± 0.03 and (I − H)_S = 1.730 ± 0.017.

5.2. Extinction and the Angular Einstein Radius

Next, we correct for extinction following the standard procedure (Yoo et al. 2004; Bennett et al. 2010) using the centroid of the red giant clump (RGC) in the CMD as a standard candle.

5.2.1. RGC Centroid Measurement

Figure 3 shows the (V − I, I) and (V − H, H) CMDs for stars within 2 arcmin of the source star. The V and I magnitudes are taken from the OGLE-III photometry catalog (Szymański et al. 2011), and the VVV (Minniti et al. 2010) catalog is calibrated to the 2MASS photometry scale for H-band magnitudes. To plot the V − H versus H CMD, we cross-referenced stars in the VVV catalog to stars in the OGLE-III catalog. For this cross-reference, we use only isolated stars that are cross-matched to within 1 arcsec of stars in the OGLE-III catalog to ensure one-to-one matching between the two catalogs. We note that the 1 arcsec limit corresponds to the average seeing in the VVV images. We find that the centroids of RGC in the (V − I, I) and (V − H, H) CMDs are I_cl = 15.33 ± 0.05, (V − I)_cl = 1.88 ± 0.02, (V − H)_cl = 4.03 ± 0.06 and (I − H)_cl = 2.11 ± 0.03.

Zoom In Zoom Out Reset image size

5.2.2. RGC Intrinsic Magnitude and Color

We use (V − I)_cl,0 = 1.06 ± 0.03 and I_cl,0 = 14.36 ± 0.05 for the intrinsic V − I color and I magnitude of the RGC (Bensby et al. 2013; Nataf et al. 2016) at the Galactic longitude of this event. Following Nataf et al. (2016), we calculate the intrinsic color of V − H and I − H in the photometric system we are using now (i.e., Johnson V, Cousins I, and 2MASS H) using the tool provided by Casagrande & VandenBerg (2014), which is based on a grid of MARCS model atmospheres (Gustafsson et al. 2008). Assuming the stellar atmospheric parameters [Fe/H] = −0.07 ± 0.10 (Gonzalez et al. 2013), $\mathrm{log}\,g\,=2.3\pm 0.1$, and [α/Fe] = 0.20 ± 0.05 (Hill et al. 2011; Johnson et al. 2014) for the RGC in the event field, we derive (V − H)_cl,0 = 2.36 ± 0.09 and (I − H)_cl,0 = 1.30 ± 0.06 by adjusting the last atmospheric parameter T_eff so that the (V − I) value is in the range of 1.03 < (V − I) < 1.09. We summarize the magnitude and colors for the RGC centroid and the source in Table 3.

Table 3.  The Source and RGC Magnitude and Colors

	I	V − I	V − H	I − H
RGC (measured from the CMD)	15.33 ± 0.05	1.88 ± 0.02	4.03 ± 0.06	2.11 ± 0.03
RGC (intrinsic)	14.36 ± 0.05	1.06 ± 0.03	2.36 ± 0.09	1.30 ± 0.06
Source (measured from the light curve)	19.54 ± 0.02	1.60 ± 0.03	3.33 ± 0.03	1.73 ± 0.02
Source (intrinsic)a	18.54 ± 0.09	0.78 ± 0.06	1.70 ± 0.11	0.92 ± 0.08

Note.

^aExtinction-corrected magnitudes using the Nishiyama et al. (2008) extinction model from Table 6.

Download table as:  ASCIITypeset image

5.2.3. Angular Einstein Radius

By subtracting the intrinsic RGC color and magnitude values from the measured RGC positions in our CMDs, we find an extinction value of A_I,obs = 0.98 ± 0.07 and color excess values of E(V − I)_obs = 0.82 ± 0.04, E(V − H)_obs = 1.67 ± 0.11, and E(I − H)_obs = 0.81 ± 0.07. Following the method of Bennett et al. (2010), we fit these values to the extinction laws of Cardelli et al. (1989) and Nishiyama et al. (2009, 2008) separately and compared the results. We present this analysis in the Appendix . From this comparison of models, we choose the Nishiyama et al. (2008) extinction law, which yields an H-band extinction of A_H = 0.19 ± 0.02 and a source angular radius of θ_* = 0.68 ± 0.02 μ as. This θ_* value implies an angular Einstein radius of ${\theta }_{E}={\theta }_{* }/\rho ={0.227}_{-0.009}^{+0.006}$ mas and a lens–source relative proper motion of ${\mu }_{\mathrm{rel}}={\theta }_{E}/{t}_{E}={4.88}_{-0.17}^{+0.14}\,\mathrm{mas}\,{\mathrm{yr}}^{-1}$.

On 2016 August 13 (HJD' = 7613.85), we observed MOA-2016-BLG-227 using the NIRC2 camera and the laser guide star (LGS) adaptive optics (AO) system mounted on the Keck II telescope at Maunakea, Hawaii. Observations were conducted in the H-band using the wide-field camera (004/pix). We took four dithered frames with 5 s exposures and three additional dithered frames with a total integration time of 90 s (6 co-adds of 15 s exposures). The first set of these images allows photometric calibration using unsaturated bright stars, and the second set provides the increased photometric sensitivity to provide a high signal-to-noise flux measurement of the target. Standard dark- and flat-field corrections were applied to the images, and sky subtraction was done using a stacked image from a nearby empty field. Each set of images was then astrometrically aligned and stacked. Finally, we use SExtracor (Bertin & Arnouts 1996) to extract the Keck source catalog from the stacked images.

A calibration catalog was extracted using an H-band image of the target area taken by the VISTA Variables in the Via Lactea survey (VVV; Minniti et al. 2010) reprocessed following the approach described in Beaulieu et al. (2016). We apply a zero point correction for the Keck source catalog using common VVV and Keck sources. The estimated zero point uncertainty is 0.05. Figure 4 shows the Keck II AO image of the field. It indicates a bright star close to the target. As a result, the dominant photometry error comes from the background flux in the wings of the PSF of the nearby star.

Zoom In Zoom Out Reset image size

We determine the source coordinates from a MOA difference image of the event while it was highly magnified. We then identify the position of the microlensing target (source+lens) on the Keck image (see Figure 4). The measured brightness of the target is H_Keck = 17.63 ± 0.06. Due to technical problems in the AO system, the stellar images display a sparse halo around each object. Thus, the FWHM of the Keck image is 0184 (measured as the average of isolated bright stars near the target). This sets a limit on our ability to exclude flux contribution from stars unrelated to the source and the lens, as we discuss below.

The light-curve analysis of the UKIRT  H-band data implies an (extinction-uncorrected) H-band source magnitude of H_S = 17.806 ± 0.017 (see Section 5.1). Because the Keck observations were taken after the event reached its baseline brightness (${t}_{\mathrm{obs},\mathrm{Keck}}-{t}_{0}=5.7\,{t}_{E}$), we can extract the excess flux by subtracting the source flux from the target flux. That is, ${H}_{\mathrm{ex},\mathrm{obs}}={H}_{S}-2.5\mathrm{log}({F}_{\mathrm{Keck}}/{F}_{S}-1)=19.7\pm 0.4$, where ${F}_{\mathrm{Keck}}/{F}_{S}={10}^{-0.4({H}_{\mathrm{Keck}}-{H}_{S})}$.

N. Koshimoto et al. (2017, in preparation) present a systematic Bayesian analysis for the identification of the star or stars producing excess flux at the position of the source seen in high angular resolution images. This analysis gives us the posterior probability distributions for the lens mass and the distance by combining the results of the light-curve modeling and the measured excess flux value. The method is summarized as follows.

• 1.  
  Determine prior probability distributions for four possibilities for the origin of the excess flux: the lens star, unrelated ambient stars, source companions, or lens companions. We denote these fluxes by F_L, F_amb, F_SC, and F_LC, respectively.
• 2.  
  Determine all combinations of the flux values for each type of star in the prior distribution that are consistent with the observed excess flux, F_excess = F_L + F_amb + F_SC + F_LC.

The extracted combinations at step 2 corresponds to the posterior probability distributions for the MOA-2016-BLG-227 event.

7.1. Prior Probability Distributions

Now, we must determine the prior probability distributions of the four types of stars that can contribute to the excess flux at the position of the source. We use all the information we have about this event—except for the value of excess flux—to create our prior probability distributions. This means that we include the FWHM of the Keck images, but not the measured magnitude of the object at the location of the microlensing event. Table 4 shows a summary of our assumptions.

Table 4.  Assumptions and Undetectable Limits Used for the Prior Probability Distributions

Priors of	Assumption	Closer Limit	Wider Limit	Used Observed Values
HL	Galactic model	...	...	tE, θE
Hamb	Luminosity function	...	0.8 FWHM	FWHM, Number density
HSC	Binary distributiona	${\theta }_{E}$/4	0.8 FWHM	FWHM, θE, HS
HLC	Binary distributiona	wLCb < u0	0.8 FWHM	FWHM, θE, u0, MLc

Notes.

^aThe binary distribution used by N. Koshimoto et al. (2017, in preparation), based on Duchêne & Kraus (2013). ^bThe caustic size created by the hypothetical companion to the lens, ${w}_{\mathrm{LC}}=4{q}_{\mathrm{LC}}/{({s}_{\mathrm{LC}}-{s}_{\mathrm{LC}}^{-1})}^{2}$. ^cThe M_L value extracted from the prior probability distributions to calculate the H_L value.

Download table as:  ASCIITypeset image

7.1.1. Lens Flux Prior

For the lens flux prior distribution, we conduct a Bayesian analysis using the observed t_E and ${\theta }_{E}$ values and the Galactic model, which has been used in a number of previous papers (Alcock et al. 1995; Beaulieu et al. 2006) to estimate lens properties for events with no microlensing parallax signal. We use the Galactic model of Han & Gould (1995) for the density and the velocity models and use the mass function presented in the Supplementary Information section of Sumi et al. (2011). Using this result and the mass–luminosity relation presented in N. Koshimoto et al. (2017, in preparation), we obtain the prior distribution for the lens apparent magnitude, H_L. We adopt the formula for the extinction to the lens, ${A}_{H,L}=(1-{e}^{-{D}_{L}/{h}_{\mathrm{dust}}})/(1-{e}^{-{D}_{S}/{h}_{\mathrm{dust}}})\,{A}_{H,S}$, following Bennett et al. (2015), where ${h}_{\mathrm{dust}}=(0.1\,\mathrm{kpc})/\sin | b|$ is the scale length of the dust toward the Galactic bulge, assuming a scale height of 0.1 kpc. Note that this Bayesian analysis gives us prior distributions for M_L, D_L, and D_S, in addition to the H_L prior distribution, but based on the observed t_E and θ_E values. These values are needed for the calculation of the probability distributions below.

7.1.2. Ambient Star Flux Prior

In order to determine the prior probability distribution for the flux of any unrelated ambient stars, we determine the number density of stars in Keck AO images, centered on the target, within a magnitude range selected to have high completeness, and divide that number by the area of the image. Then we use the luminosity function of Zoccali et al. (2003) to derive the number density of stars as a function of H magnitude, normalized to this measured number density in the Keck AO image. In this calculation, we correct for the differences in extinction and distance moduli between our field and that of Zoccali et al. (2003), using the distance moduli from Table 3 of Nataf et al. (2013) and extinction values for both fields.

When correcting for the extinction difference, we also consider the difference between the extinction laws used. Zoccali et al. (2003) derived an A_H value using the C89 extinction law with R_V = 3.1, whereas our preferred N08 extinction law implies a significantly different A_H value. To correct for this difference, we calculate the A_H value toward their field using the N08 extinction-law fit to the RGC centroid in the OGLE-III CMD and the R_JKVI value from Table 3 of Nataf et al. (2013) for their field. The A_H value we derived here is A_H = 0.122, which is different from the value of A_H = 0.265 used by Zoccali et al. (2003). Therefore, we convert their extinction-corrected H-band luminosity function to a luminosity function with our preferred extinction model by adding ΔA_H = 0.265–0.122 = 0.142 to their extinction-corrected magnitudes, and then add the extinction appropriate for our field, A_H = 0.19.

We assume that stars can be resolved only if they are separated from the source by ≥0.8 FWHM = 148 mas. Under this assumption, the expected number of ambient stars within the circle is derived by multiplying the area of this unresolvable region by the total number density derived above. We determine the number of stars following the Poisson distribution with the mean value of the expected number of stars. We use the corrected luminosity function to determine the magnitude of each star.

7.1.3. Source and Lens Companion Flux Priors

We calculate the source and lens companion flux priors with the stellar binary distribution described in N. Koshimoto et al. (2017, in preparation). The binary distribution is based on the summary in a review paper (Duchêne & Kraus 2013), which provides THE distributions of the stellar multiplicity fraction, and the mass ratio and semimajor axis distributions.

For the flux of the source companions, we calculate the source mass M_SC = q_SC M_S and then convert that into a source companion magnitude, H_SC, using a mass–luminosity relation. The mass ratio q_SC is derived from the binary distribution. We derive the source mass, ${M}_{S},$ from the combination of H_S, D_S, and using the mass–luminosity relation. Similarly, we calculate the lens companions magnitude, H_LC, from M_LC = q_LC M_L, where the lens mass M_L comes from the same distribution that was used to obtain the lens flux probability distribution.

We consider companions to the lens or source located in the same unresolvable regions in the vicinity of the source, just as in the case of ambient stars. Stellar companions have a separation distribution that is much closer to logarithmic than the uniform distribution expected for ambient stars. As a result, we must now exclude companions that are too close to the source and lens as well as companions that are so widely separated that they will be resolved. Companions that are too close to the source could be magnified themselves, and companions that are too close to the lens could serve as an additional lens star. Such a constraint would have no effect on the ambient star probability, because the probability of an ambient star very close to the source or lens is much smaller than that of a stellar companion. Following Batista et al. (2014), we adopt θ_E/4 as the close limit for source companions and ${w}_{\mathrm{LC}}\lt {u}_{0}$ as the close limit for lens companions, where ${w}_{\mathrm{LC}}=4{q}_{\mathrm{LC}}/{({s}_{\mathrm{LC}}-{s}_{\mathrm{LC}}^{-1})}^{2}$ (Chung et al. 2005), and q_LC and s_LC and are the stellar binary lens mass ratio and separation, respectively. We take 0.8 FWHM as the maximum unresolvable radius.

We also consider triple and quadruple systems when estimating the effect of companions to the source and lens, following N. Koshimoto et al. (2017, in preparation), but we find no significant difference from the case that only considers binary systems. We therefore do not include triple and quadruple systems in this analysis, for simplicity.

7.1.4. Excess Flux Prior

Figure 5 shows the prior probability distributions we derived following the procedure described above to calculate the flux of each type of star. In addition to the magnitude of the four types of stars that might contribute to the excess flux, we show the prior distributions for the total excess flux H_excess, the lens mass M_L, and the distance to the lens D_L. Some of the panels in this figure have total probabilities P_total < 1. This is because many stars do not have binary companions, and there is a large probability of no measurable flux from an ambient star. The H_excess prior indicates a high probability at the observed magnitude of ${H}_{\mathrm{ex},\mathrm{obs}}=19.7\pm 0.4$. The three panels for individual stars, H_L, H_amb, and H_SC, show similar probabilities at the observed excess flux value. This indicates that it will be difficult to claim that all of the excess comes from the lens itself.

Zoom In Zoom Out Reset image size

7.2. Posterior Probability Distributions

We generate the posterior probability distributions shown in Figure 6 by extracting combinations of parameters, which have values of H_excess consistent with the measured value of H_ex,obs = 19.7 ± 0.4 using a Gaussian distribution in fluxes (not magnitudes). The probability that H_L ≤ 20 is almost the same as the probability for H_SC ≤ 20 and slightly higher than, but competitive with, the probability that H_amb ≤ 20, which results in very loose constraints on H_L and M_L. This result is consistent with our expectation as discussed in Section 7.1.4.

Zoom In Zoom Out Reset image size

The third to sixth columns of Table 5 shows the median, the 1σ error bars, and the 2σ range for H_L, M_L, and D_L for both the prior and posterior distributions. This same table also shows the values of the planet mass M_p, the projected separation a_⊥, and the three-dimensional star–planet separation a_3d calculated from the probability distributions, where a_3d is statistically estimated assuming a uniform orientation for the detected planets. In the bottom three rows, we present the probabilities that the fraction of the excess flux due to the lens, f_L, is larger than 0.1, 0.5, and 0.9, which corresponds to the magnitude difference between the lens and the total flux excess of 2.5 mag, 0.75 mag, and 0.11 mag, respectively.

Table 5.  Lens Properties Calculated from the Prior and Posterior Probability Distributions

Parameters	Unit	Prior	2σ range	Posterior	2σ range	Posterior	Posterior	Posterior
		α = 0		α = 0		α = 1	α = 2	α = 3
HL	mag	${21.9}_{-1.6}^{+1.3}$	18.0-28.7	${21.5}_{-1.3}^{+1.4}$	19.5-27.5	${20.8}_{-0.9}^{+1.3}$	${20.4}_{-0.7}^{+1.1}$	${20.2}_{-0.6}^{+0.8}$
ML	M⊙	${0.24}_{-0.12}^{+0.27}$	0.06-0.90	${0.29}_{-0.15}^{+0.23}$	0.07-0.66	${0.42}_{-0.20}^{+0.17}$	${0.50}_{-0.18}^{+0.13}$	${0.54}_{-0.15}^{+0.10}$
Mp	MJup	${2.3}_{-1.2}^{+2.6}$	0.5-8.8	${2.8}_{-1.5}^{+2.2}$	0.6-6.4	${4.1}_{-1.9}^{+1.7}$	${4.8}_{-1.8}^{+1.2}$	${5.3}_{-1.5}^{+1.0}$
DL	kpc	6.4 ± 1.0	3.5-8.5	6.5 ± 1.0	3.9-8.6	${6.8}_{-0.9}^{+1.0}$	${6.9}_{-0.9}^{+1.0}$	${7.1}_{-0.9}^{+1.0}$
a⊥	au	1.37 ± 0.23	0.76-1.84	1.39 ± 0.22	0.84-1.86	${1.45}_{-0.20}^{+0.22}$	${1.49}_{-0.20}^{+0.22}$	${1.51}_{-0.20}^{+0.22}$
a3d	au	${1.64}_{-0.36}^{+0.93}$	0.89-6.49	${1.67}_{-0.35}^{+0.94}$	0.97-6.62	${1.74}_{-0.35}^{+0.99}$	${1.79}_{-0.35}^{+1.02}$	${1.82}_{-0.36}^{+1.03}$
P(fL > 0.1)a	%	72.2	...	78.1	...	90.8	96.5	98.7
$P({f}_{L}\gt 0.5)$ a	%	48.0	...	41.4	...	56.7	69.6	78.4
P(fL > 0.9)a	%	33.8	...	24.0	...	29.9	38.0	44.5

Note. The values of the posterior probability distributions are shown also for different α values, the slope of the probability of hosting planets ${P}_{\mathrm{host}}\propto {M}^{\alpha }$. The values given in the form of the median with the 1σ uncertainty. The 2σ range is given for α = 0.

^aThe probabilities that the fraction of the lens flux to the excess flux, f_L ≡ F_L/F_excess, is larger than the indicated values. The fractions of 0.1, 0.5, and 0.9 correspond to the difference of magnitude, ${H}_{L}-{H}_{\mathrm{excess}}=-2.5\,\mathrm{log}({F}_{L}/{F}_{\mathrm{excess}})$, of 2.5 mag, 0.75 mag, and 0.11 mag.

Download table as:  ASCIITypeset image

The posterior distributions for the lens system properties are remarkably similar to the prior distributions. When we compare the 1σ ranges of the prior and posterior distributions, we see that the lens system is most likely to be composed of an M or K-dwarf star host and a gas-giant planet. However, the prior and posterior distributions differ from each other when we consider the 2σ ranges and the tails of the distributions. The possibility of a G dwarf host star is ruled out by the posterior distribution while the host star can be a G dwarf according to the prior distribution. This implies that the host star is likely to be an M or K dwarf.

7.3. Comparison of Different Planetary Host Priors

One assumption that we have made implicitly is that the properties of the lens star do not depend on the fact that we have detected a planet orbiting the star. This assumption could be false. Perhaps more massive stars are more likely to host planets of the measured mass ratio, or perhaps disk stars are more likely to host planets than bulge stars. The microlensing method can be used to address these questions, but we must be careful not to assume the answer to them.

We have assumed that this detection of the planetary signal does not bias any other property of the lens star, such as its mass or distance. If there was a strong dependence of the planet-hosting probability at the measured mass ratio of ${9.3}_{-0.1}^{+0.2}\times {10}^{-3}$, then this implicit prior could lead to incorrect conclusions. Some theoretical papers based on core accretion (Laughlin et al. 2004; Kennedy & Kenyon 2008) and analyses of exoplanets found by radial velocities (Johnson et al. 2010) have argued that gas giants are less frequently orbiting low-mass stars; however, the difference disappears when the planets are classified by their mass ratio, q, instead of their mass. Nevertheless, since the host-mass dependence of the planet-hosting probability is not well measured, we investigate how our results depend on the choice of this prior.

We consider a series of prior distributions where the planet-hosting probability follows a power law of the form P_host ∝ M^α, and we conduct a series of Bayesian analyses with α = 1, α = 2, and α = 3 in addition to the calculation with α = 0, presented above. Figure 7 shows both the prior and posterior probability distributions for the lens mass, M_L, with these different values of α. The lens property values for each posterior distribution are shown in Table 5. The median of expected lens flux approaches the measured excess flux as α increases (i.e., the power law becomes steeper), and consequently the median of the lens mass also increases and the parameter uncertainties decrease. Thus, larger α values imply that more of the excess flux is likely to come from the lens. Nevertheless, our basic conclusion that the host is an M or K dwarf hosting a gas-giant planet remains for all of the 1 ≤ α ≤ 3 priors.

Zoom In Zoom Out Reset image size

We have analyzed the planetary microlensing event MOA-2016-BLG-227, which was discovered next to the field observed by the microlensing campaign (Campaign 9) of the K2 Mission. The event and planetary signal were discovered by the MOA collaboration and a significant portion of the planet signal was covered by the data from the Wise, UKIRT, CFHT, and VST surveys, which observed the event as part of the K2C9 program. Analysis of these data yields a unique microlensing light-curve solution with a relatively large planetary mass ratio of $q={9.28}_{-0.11}^{+0.20}\times {10}^{-3}$. We considered several different extinction laws and decided that the N08 (Nishiyama et al. 2008) law was the best fit to our data, although our results would not change significantly with a different law. With this extinction law, we derive an angular Einstein radius of ${\theta }_{E}={0.227}_{-0.009}^{+0.006}\,\mathrm{mas}$, which yields the mass–distance relation given in Equation (3). We detected excess flux at the location of the source in a Keck AO image, and we performed a Bayesian analysis to estimate the relative probability of the different sources of this excess flux, such as the lens, an ambient star, or a companion to the host or source. Our analysis excludes the possibility that the host star is a G dwarf, leading us to the conclusion that the planet MOA-2016-BLG-227Lb is a super-Jupiter mass planet orbiting an M- or K-dwarf star likely located in the Galactic bulge. Such systems are predicted to be rare by the core accretion theory of planet formation. It is also thought that such a planet orbiting a white dwarf host at ${a}_{3d}\sim 2\,\mathrm{au}$ is unlikely (Batista et al. 2011).

If the planet frequency does not depend on the host star mass or distance, our Bayesian analysis indicates the system consists of a host star with mass of ${M}_{L}={0.29}_{-0.15}^{+0.23}{M}_{\odot }$ orbited by a planet with mass of ${M}_{{\rm{p}}}={2.8}_{-1.5}^{+2.2}{M}_{\mathrm{Jup}}$ with a three-dimensional star–planet separation of ${a}_{3d}={1.67}_{-0.35}^{+0.94}\,\mathrm{au}$. The system is located at D_L = 6.5 ± 1.0 kpc from the Sun. We also considered different priors for the planet-hosting probability as a function of host star mass. We consider planet-hosting prior probabilities that scale as P_host ∝ M^α with α = 1, 2, 3, in addition to the α = 0 prior that we use for our main results. As α increases, the median value of the lens mass also increases and the probability for the lens to be responsible for the excess H-band flux increases as well. Johnson et al. (2010) found a linear (i.e., α = 1) relationship between host mass and planet occurrence from 0.5 M_⊙ to 2.0 M_⊙ for giant planets within ∼2 au around host stars discovered by the radial velocity (RV) method. However, this analysis used a fixed minimum mass instead of a fixed mass ratio, and it does not appear that Johnson et al. (2010) did a detailed calculation of their detection efficiencies. Another result using RV planet data by Montet et al. (2014) gives $\alpha ={0.8}_{-0.9}^{+1.1}$, using a sample more similar to the microlensing planets, i.e., gas giants orbiting at 0 < a < 20 au around M-dwarf stars. However, our basic conclusion that the MOA-2016-BLG-227L host star is an M or K dwarf with a gas-giant planet located in the Galactic bulge would not change with a different α value, as indicated in Figure 7 and Table 5.

The probability that more than 90% of the excess flux seen in the Keck AO images comes from the lens is still 24.0% even assuming α = 0. This is significant enough that we cannot ignore the possibility that most of the excess flux comes from the lens star. One approach for obtaining further constraints is to get the color of the excess flux. If the excess flux is not from the lens, the derived lens mass and distance with H_excess may be inconsistent with the value derived using the excess flux in a different passband, if we assume that all of the excess flux comes from the lens. However, such a measurement could also yield ambiguous results. Another, more definitive, approach is to observe this event in the future when we can expect to detect the lens–source separation through precise PSF modeling with high-resolution space-based data (Bennett et al. 2007, 2015) or direct resolution with AO imaging (Batista et al. 2015). The lens–source relative proper motion value of ${\mu }_{\mathrm{rel}}={4.88}_{-0.17}^{+0.14}\,\mathrm{mas}\,{\mathrm{yr}}^{-1}$ indicates that we can expect to be able to resolve the lens, if it provides a large fraction of the excess flux in ∼2022 using HST (Bhattacharya et al. 2017) and in 2026 using Keck AO (Batista et al. 2015). Observations by the James Webb Space Telescope (Gardner et al. 2006), the Giant Magellan Telescope (Johns et al. 2012), the Thirty Meter Telescope (Nelson & Sanders 2008), and the Extremely Large Telescope (Gilmozzi & Spyromilio 2007) could detect the lens–source relative proper motion much sooner. If the separation of the excess flux from the source is different from the prediction of the microlensing model in these future high angular resolution observations, it would indicates that the lens is not the main cause of the excess flux, implying a lower mass planetary host star.

Work by N.K. is supported by JSPS KAKENHI Grant Number JP15J01676. The MOA project is supported by grants JSPS25103508 and 23340064. D.P.B., A.B., and D.S. were supported by NASA through grants NNX13AF64G, NNX16AC71G, and NNX16AN59G. Work by Y.S. and C.B.H. was supported by an appointment to the NASA Postdoctoral Program at the Jet Propulsion Laboratory, California Institute of Technology, administered by Universities Space Research Association through a contract with NASA. N.J.R. is a Royal Society of New Zealand Rutherford Discovery Fellow. Work by A.F. is supported by JSPS KAKENHI grant No. JP17H02871. Work by C.R. was supported by an appointment to the NASA Postdoctoral Program at the Goddard Space Flight Center, administered by USRA through a contract with NASA. Work by S.C.N. was supported by NExScI. A.C. acknowledges financial support from Université Pierre et Marie Curie under grants Émergence@Sorbonne Universités 2016 and Émergence-UPMC 2012. This work was supported by a NASA Keck PI Data Award, administered by the NASA Exoplanet Science Institute. Data presented herein were obtained at the W. M. Keck Observatory from telescope time allocated to the National Aeronautics and Space Administration through the agency's scientific partnership with the California Institute of Technology and the University of California. The Observatory was made possible by the generous financial support of the W.M. Keck Foundation. The authors wish to recognize and acknowledge the very significant cultural role and reverence that the summit of Maunakea has always had within the indigenous Hawaiian community. We are most fortunate to have the opportunity to conduct observations from this mountain. The United Kingdom Infrared Telescope (UKIRT) is supported by NASA and operated under an agreement among the University of Hawaii, the University of Arizona, and Lockheed Martin Advanced Technology Center; operations are enabled through the cooperation of the Joint Astronomy Centre of the Science and Technology Facilities Council of the U.K. We acknowledge the support from NASA HQ for the UKIRT observations in connection with K2C9. This research uses data obtained through the Telescope Access Program (TAP), which has been funded by the National Astronomical Observatories of China, the Chinese Academy of Sciences (the Strategic Priority Research Program "The Emergence of Cosmological Structures" Grant No. XDB09000000), and the Special Fund for Astronomy from the Ministry of Finance. This work is partly based on observations obtained with MegaPrime/MegaCam, a joint project of CFHT and CEA/DAPNIA, at the Canada–France–Hawaii Telescope (CFHT), which is operated by the National Research Council (NRC) of Canada, the Institut National des Science de l'Univers of the Centre National de la Recherche Scientifique (CNRS) of France, and the University of Hawaii. This work was performed in part under contract with the California Institute of Technology (Caltech)/Jet Propulsion Laboratory (JPL) funded by NASA through the Sagan Fellowship Program executed by the NASA Exoplanet Science Institute. Work by M.T.P. and B.S.G. was supported by NASA grant NNX16AC62G. This work was partly supported by the National Science Foundation of China (grant No. 11333003, 11390372 to S.M.). Based on observations made with ESO Telescopes at the La Silla Paranal Observatory under programme ID 097.C-0261.

In Section 5.2, we obtained the observed extinction value A_I,obs = 0.98 ± 0.07, and color excess values of E(V − I)_obs = 0.82 ± 0.04, E(V − H)_obs = 1.67 ± 0.11, and E(I − H)_obs = 0.81 ± 0.07. Then, we fit these values to the extinction laws of Cardelli et al. (1989) and Nishiyama et al. (2009, 2008) separately and compared the results. This was motivated by the fact that Nataf et al. (2016) reported a clear difference of their extinction law toward the Galactic bulge from the standard law of Cardelli et al. (1989). Hereafter, we refer to these papers as C89, N09, and N08, respectively. Note that the four observed extinction parameters (1 extinction and 3 color excess) are not independent. They can be derived from the three independent extinction values: ${A}_{I,\mathrm{obs}}$, A_V,obs, and A_H,obs.

The C89 law is given by Equations (1)–(3b) in their paper, and A_V and R_V serve as the parameters of their model.

Unlike C89, N09 does not provide a complete extinction model. They provide only ratios of extinctions for wavelengths longer than the J-band. So, we need additional information relating A_V or A_I and A_J, A_H, or A_K in order to calculate the values that we need for this paper: A_I, A_V, and A_H. Therefore, we used the ${R}_{{JKVI}}\equiv E(J-{K}_{s})/E(V-I)$ values from Nataf et al. (2013) in addition to the N09 extinction law. The R_JKVI value at the nearest grid point to the MOA-2016-BLG-227 event in Table 3 of Nataf et al. (2013) is 0.3089. However, the quality flag for this value is 1, which indicates an unreliable measurement, so we use a conservative uncertainty of R_JKVI = 0.31 ± 0.03. We adjust A_I and E(V − I) to minimize the χ² value between the observed A_I,obs, E(V − I)_obs, E(V − H)_obs, and E(I − H)_obs values and those values derived using the ratio ${A}_{H,2\mathrm{MASS}}/E{(J-{K}_{s})}_{2\mathrm{MASS}}=0.89$ from N09, in conjunction with the R_JKVI value from Nataf et al. (2013). We explicitly use the _2MASS subscript because N09 also provide their result in the IRSF/SIRIUS photometric system (Nagashima et al. 1999; Nagayama et al. 2003). Note that we calculate this ${A}_{H}/E(J-{K}_{s})$ value using their result for the field S+ (0° < l < 3°, −1° < b < 0°), which is the nearest of their fields to the MOA-2016-BLG-227 event position.

N08 also provide the ratio of extinctions toward the Galactic bulge (l ∼ 0°, b∼−2°). They find A_J/A_V = 0.183 ± 0.015, A_H/A_V = 0.103 ± 0.008, and ${A}_{{K}_{s}}/{A}_{V}=0.064\pm 0.005$. (These values are slightly different from the original values given by N08 because the values used in N08 were in the OGLE II and IRSF/SIRIUS photometric systems, so we converted them into the standard systems that we use here.) These values are well fit by a single power law, A_λ/A_V ∝ λ⁻². Nevertheless, we use the ratios themselves, instead of the single power law, because N08 does not test that their power law accurately reproduces A_I/A_V, which we have now. As in the case of N09, we keep these ratios fixed, and adjust A_I and E(V − I) to minimize the χ² between these relations and the observed A_I,obs, E(V − I)_obs, E(V − H)_obs, and E(I − H)_obs values. Notice that N08 had V-band data, and it was not necessary to use R_JKVI as a constraint. Therefore, we used the R_JKVI value as the additional observed data instead here in addition to ${A}_{I,\mathrm{obs}}$, A_V,obs, and A_H,obs to increase the number of degrees of freedom (dof).

Table 6 shows the results of fitting our extinction measurements to these three different extinction laws. This table also shows the angular source radius calculated from the extinction-corrected source magnitudes and colors using formulae from the analysis of Boyajian et al. (2014). We determine ${\theta }_{* ,{IH}}$ using Equations (1)–(2) and Table 1 of Boyajian et al. (2014), but the other relations were provided by private communications from Boyajian with a special analysis restricted to stellar colors that are relevant to the Galactic bulge sources observed in microlensing events. We use Equation (4) of Fukui et al. (2015) to determine ${\theta }_{* ,{VI}}$, and we use Equation (4) of Bennett et al. (2015) to determine ${\theta }_{* ,{VH}}$. Those formulae are

$$\begin{eqnarray}&&\mathrm{log}\,[2{\theta }_{* ,{VI}}/(1\,\mathrm{mas})]\\ &&\,=\,0.5014+0.4197{(V-I)}_{S,0}-0.2{I}_{S,0},\end{eqnarray} \tag{ 6 }$$

$$\begin{eqnarray}&&\,\mathrm{log}\,[2{\theta }_{* ,{VH}}/(1\,\mathrm{mas})]\\ &&\,=\,0.5367+0.0727{(V-H)}_{S,0}-0.2{H}_{S,0},\end{eqnarray} \tag{ 7 }$$

$$\begin{eqnarray}&&\mathrm{log}\,[2{\theta }_{* ,{IH}}/(1\,\mathrm{mas})]\\ &&\,=\,0.5303+0.3660{(I-H)}_{S,0}-0.2{I}_{S,0}.\end{eqnarray} \tag{ 8 }$$

Table 6.  Comparison of the Extinction and Angular Einstein Radius Based on Different Extinction Laws

Extinction Law	Nonea	Cardelli et al. (1989)	Nishiyama et al. (2009)	Nishiyama et al. (2008)
Relation	...	$\tfrac{{A}_{\lambda }}{{A}_{V}}=a(x)+\tfrac{b(x)}{{R}_{V}}$ b	$\tfrac{{A}_{H}}{E(J-{K}_{s})}$, RJKVIc	$\tfrac{{A}_{J}}{{A}_{V}}$, $\tfrac{{A}_{H}}{{A}_{V}}$, $\tfrac{{A}_{{K}_{s}}}{{A}_{V}}$
AV	1.80 ± 0.08	1.87 ± 0.12	1.83 ± 0.12	1.82 ± 0.12
AI	0.98 ± 0.07	1.04 ± 0.10	1.01 ± 0.08	1.00 ± 0.08
AH	0.17 ± 0.10	0.30 ± 0.04	0.23 ± 0.04	0.19 ± 0.02
E(V − I)	0.82 ± 0.04	0.83 ± 0.05	0.82 ± 0.05	0.82 ± 0.05
E(V − H)	1.67 ± 0.11	1.57 ± 0.10	1.60 ± 0.11	1.63 ± 0.11
E(I − H)	0.81 ± 0.07	0.73 ± 0.06	0.78 ± 0.09	0.81 ± 0.08
χ2/dofd	...	2.39/1	0.56/1	1.19/2e
θ*,VI (μas)	0.65 ± 0.04	0.67 ± 0.06	0.66 ± 0.05	0.66 ± 0.05
θ*,VH (μas)	0.67 ± 0.04	0.73 ± 0.02	0.70 ± 0.02	0.68 ± 0.02
θ*,IH (μas)	0.71 ± 0.07	0.79 ± 0.06	0.74 ± 0.06	0.72 ± 0.06
SDθ*f	0.035	0.061	0.042	0.030
θE (mas)g	...	...	...	${\bf{0}}.{{\bf{227}}}_{-0.009}^{+0.006}$
μrel (mas yr−1)g	...	...	...	${\bf{4}}.{{\bf{88}}}_{-0.17}^{+0.14}$

Notes. The values in boldface are used as final values.

^aResult without using an extinction law. The A_I, E(V − I), E(V − H), and E(I − H) values are determined directly from the data. ^bEquation (1) of C89; see the paper for the detailed model. ^cThe R_JKVI value comes from Table 3 of Nataf et al. (2013). ^dWhen calculating the total χ², we multiply each of the contributions from E(V − I), E(V − H), and E(I − H) by 2/3, because these values are not independent. ^eThe dof = 2 is because we used the R_JKVI value from Nataf et al. (2013) as an observed data point. ^fStandard deviation of the three θ_* values. ^gCalculations conducted only for the adopted θ_* value (θ_*,VH with N08).

Download table as:  ASCIITypeset image

If we compare the χ² value for each model fit in Table 6, we see that the ${\chi }^{2}/\mathrm{dof}$ for the N09 and N08 laws are smaller than the value from the C89 extinction law, although the C89 is not disfavored by a statistically significant amount. (The p-value of χ² = 2.39 for dof = 1 is still ∼0.12.) Note that a contribution of 0.96 to the total value of χ² = 1.19 arises from fitting the R_JKVI value to the N08 extinction law. So, the remaining contribution of 0.23 to ${\chi }^{2}$ arises from fitting the N08 model to our measurements of the RGC centroids. This indicates that the extinction law of N08 agrees with our measurement of the red clump centroids very well, but not quite so well with the R_JKVI value, which comes from Nataf et al. (2013).

From the point of view of consistency between the three θ_* values, the standard deviation of the three values (SD${}_{{\theta }_{* }}$ in the table) is smallest using the N08 extinction laws. The N08 extinction law also yields the smallest error bars for A_H and ${\theta }_{* ,{VH}}$.

Based on this analysis, we have decided to use the results from the N08 extinction laws in our analysis. We use ${\theta }_{* ,{VH}}$ for the final angular source radius, which is θ_* = 0.68 ± 0.02 μas. We show the source magnitudes and colors corrected for extinction using the N08 extinction laws in Table 3.