MOA-2009-BLG-319Lb: A Sub-Saturn Planet inside the Predicted Mass Desert

Gravitational microlensing has the unique ability to detect cold exoplanets beyond the snow line (Mao & Paczynski 1991; Gould & Loeb 1992) and down to Earth masses (Bennett & Rhie 1996). So far microlensing has detected ∼100 planets at distances up to the Galactic bulge. One drawback of this method is that for most light curves only the mass ratio of the lens system is measured, which leaves some physical parameters of the system significantly unconstrained. This results in large estimated uncertainties, particularly in the inferred stellar host and companion masses, due to uncertain priors used in the standard Bayesian modeling approach. One can mitigate this limitation by resolving the source and lens independently with high angular resolution imaging (i.e., Hubble Space Telescope (HST), Keck adaptive optics (AO), Subaru AO) several years after peak magnification, for which Bennett et al. (2006, 2007) laid the theoretical groundwork. This high angular resolution imaging allows us to further constrain the lens–source separation, relative proper motion between the targets, and lens flux, which can then be used with mass–luminosity relations (Henry & McCarthy 1993; Henry et al. 1999; Delfosse et al. 2000) to infer a direct mass for the host.

Several microlensing source and lens stars have now been measured with these techniques, beginning with OGLE-2005-BLG-169 (Batista et al. 2015; Bennett et al. 2015). These follow-up observations from Keck II and HST confirmed, for the first time, the planetary interpretation from the light curve by verifying the lens–source relative proper motion as predicted by the original light-curve measurement. The host star mass was precisely determined to be 0.69 ± 0.02M_⊙, with a planetary companion of mass 14.1 ± 0.9M_⊕.

This current analysis is part of the NASA Keck Key Strategic Mission Support (KSMS) program, "Development of the WFIRST Exoplanet Mass Measurement Method" (Bennett 2018), which is a pathfinder project for the Nancy Grace Roman Space Telescope (formerly known as WFIRST; Spergel et al. 2015). A large fraction of the Roman Telescope observing time will be devoted to the Roman Galactic Exoplanet Survey (RGES), which is a dedicated microlensing survey (Bennett & Rhie 2002; Bennett et al. 2010a; Penny et al. 2019; Johnson et al. 2020) that will complement previous large statistical studies of transiting planets from the Kepler telescope (Borucki et al. 2011), among others. The KSMS program has already measured the masses of several microlensing host stars and their planetary companions (Bhattacharya et al. 2018; Bennett et al. 2020; Vandorou et al. 2020). Several more lens system mass measurements from the KSMS program are in preparation (Bhattacharya et al. 2021; C. Ranc et al. 2020, in preparation; Blackman et al. 2021). A majority of the targets observed in this program were included in the statistical sample of Suzuki et al. (2016, 2018), which shows a break and likely peak in the mass ratio function for wide-orbit planets at about a Neptune mass. This study is the most complete statistical sample of microlensing planets to date, and the results are seemingly at odds with the runaway gas accretion scenario of the leading core accretion theory of planet formation (Lissauer 1993; Pollack et al. 1996), which predicted a planet desert at sub-Saturn masses (Ida & Lin 2004) for gas giants at wide orbits. Suzuki et al. (2018) studied only the exoplanet mass ratio, q, so they could not determine whether there was a gap over part of the host mass range. For example, since the core accretion theory was primarily developed with solar-type host stars in mind, the gap expected from the runaway gas accretion scenario might exist for solar-type stars but be washed out with the low-mass M-dwarf hosts that are also included in the microlens sample. Mass measurements like the one presented in this paper can probe this possibility.

This paper is organized as follows: Section 2 describes the original observations for MOA-2009-BLG-319. In Section 3 we perform improved photometry of the light curve and present an updated analysis of the light curve. In Section 4, we describe the Keck AO follow-up analysis and a new MCMC routine for precise astrometry in Keck AO imaging. Section 5 details our lens–source relative proper-motion measurements. Section 6 describes the lens system properties with new constraints from Keck high-resolution imaging. Finally, we discuss the results and conclude the paper in Section 7.

MOA-2009-BLG-319, located at R.A. = 18:06:58.026, decl. = −26:49:10.945 and with Galactic coordinates (l, b) = (4.202, −3.014), was first alerted by the Microlensing Observations in Astrophysics (MOA; Bond et al. 2001; Sumi et al. 2003) collaboration on 2009 June 20. MOA initially reported "low-level systematics" in their observations shortly after continuous monitoring began. This light-curve feature turned out to be the first of several planetary caustic crossings throughout the duration of this high-magnification event. At the time of publication, MOA-2009-BLG-319 (Miyake et al. 2011) had the best-sampled light curve of all observed microlensing events.

Our photometry methods have improved since the Miyake et al. (2011) analysis, so we have re-reduced the photometry for a number of the data sets. We have used the method of Bond et al. (2001, 2017) to reduce the data from the MOA-II telescope, the Mt. John Observatory Boller and Chivens 0.61 m telescope (operated by the MOA group), and the SMARTS telescope at CTIO. The MOA-II data were corrected for systematic errors due to chromatic differential refraction (Bennett et al. 2012). The SMARTS-CTIO data were previously reduced with DoPHOT (Schechter et al. 1993), but the difference imaging photometry that we provide (Bond et al. 2001, 2017) is well known to be a substantial improvement. New reductions are also needed to provide a Markov Chain Monte Carlo (MCMC) distribution to understand the distribution of models that are consistent with the data.

While more than 20 data sets were used for the original paper, many of these do not actually constrain the light-curve model. Therefore, we fit only to the following data sets: the MOA-II red band; the MOA 0.61 m Boller and Chivens V and I bands; SMARTS-CTIO V, I, and H bands; the Robonet Faulkes telescope (north and south) I band; the Liverpool telescope I band; and the Bronberg Observatory unfiltered data. Figure 1 shows the best-fit model with the data used in this paper, except for the sparsely sampled V-band data. The CTIO data were taken with the ANDICAM instrument of the SMARTS-CTIO telescope, which takes optical and infrared data simultaneously. The infrared data from this telescope are known to occasionally display systematic errors between images taken at the five different dither positions, which are apparently due to subpixel-scale sensitivity variations (Dong et al. 2009a). Therefore, we treat the data from these different dither positions as independent data sets, shown in different shades of green in Figure 1 as CTIO-H0 through CTIO-H4.

Zoom In Zoom Out Reset image size

The light-curve modeling follows the image-centered ray-shooting method of Bennett & Rhie (1996) and Bennett (2010). Figure 1 shows our best-fit planetary model for this event, and Table 1 shows the parameters of our best-fit model, as well as the MCMC averages of models consistent with the data. These are also compared to the distribution from the original study of Miyake et al. (2011).

A follow-up light-curve analysis by Shin et al. (2015) considered two-planet models for MOA-2009-BLG-319 and a number of other planetary microlensing events, and their analysis found a significant χ² improvement, Δχ² > 100, for their best two-planet model for this event. However, this analysis was incomplete, as they did not consider other triple-lens models for this event. The analysis of planetary microlensing event OGLE-2007-BLG-349 indicates that circumbinary models can describe deviations that are also consistent with two-planet models (Bennett et al. 2016), and there can also be degeneracies between circumbinary planet models and circumstellar planet models in binary systems (Gould et al. 2014). We will not consider these triple-lens models further in this paper, as the analysis of these triple-lens models is not complete. We should note, however, that if the two-planet model is correct, then the conclusions of this paper will be unchanged except that there will be an additional, lower-mass planet. Also, these triple-lens models are relevant for the consideration of a microlensing parallax signal. While the MOA-2009-BLG-319 Einstein radius crossing time is too short to expect a microlensing parallax signal due to the orbital motion of Earth, the dense coverage of the light-curve peak by widely separated observatories suggests the possibility of a terrestrial microlensing parallax signal (Hardy & Walker 1995; Holz & Wald 1996; Gould et al. 2009), as pointed out by Miyake et al. (2011). However, the triple-lens models will affect the same part of the light curve. Thus, it would not be useful to investigate any microlensing parallax solution without also considering a third lens mass.

In order to determine the source radius, we need to determine the extinction-corrected source magnitude and color. Miyake et al. (2011) used the SMARTS-CTIO V- and I-band data for this. However, these SMARTS-CTIO data were reduced with DoPHOT, and this has occasionally led to magnitude and color measurements that led to spurious conclusions about the properties of planetary microlens systems (Bennett et al. 2017). This is why it was necessary to use the difference imaging code and calibration method of Bond et al. (2017) for this reanalysis of the SMARTS-CTIO V and I-band data. Also, predicted properties of the bulge red clump giant stars that are used to determine the extinction have changed since the Miyake et al. (2011) analysis. We have calibrated the SMARTS-CTIO V- and I-band data to the OGLE-III catalog (Szymański et al. 2011), and then we located the red clump centroid at V_rc − I_rc = 1.98, I_rc = 15.44, following the method of Bennett et al. (2010b). Using the bulge red clump giant magnitude, color, and distance from Nataf et al. (2013), we find I- and V-band extinction of A_I = 1.116 and A_V = 2.036. Using the source magnitudes from Table 1, we find extinction-corrected magnitudes of I_S0 = 18.878 ± 0.069 and V_S0 = 19.678 ± 0.069. This allows us to use the surface brightness relation from the analysis of Boyajian et al. (2014), but we use the following custom formula (Bhattacharya et al. 2016) using stars spanning the range in colors that are relevant for microlensing events:

$$\begin{eqnarray}&&\mathrm{log}(2{\theta }_{* })=0.5014+0.4197({V}_{s0}-{I}_{s0})-0.2{I}_{s0}.\end{eqnarray} \tag{ 1 }$$

This yields θ_* = 0.576 ± 0.077 μas, which is smaller than the Miyake et al. (2011) value of θ_* = 0.66 ± 0.06 μas. Our measurement is consistent with the μ_rel measurement from Keck. This difference from the value that Miyake et al. (2011) find is due in part to the combination of the error in magnitude from DoPHOT and an improved knowledge of the red clump from Nataf et al. (2013) as described earlier.

To measure our new lens system parameters, we sum over the MCMC results using a Galactic model (Bennett et al. 2014) with weights for the microlensing rate and our μ_rel,H value from Keck (described in Section 5). We constrain the possible source distances to follow the weighted distribution from the microlensing event rate in our Galactic model, which results in a best-fit source distance of D_S = 8.25 ± 0.86 kpc. These new light-curve modeling results produce smaller best-fit values for the mass ratio, q, and angular Einstein radius, θ_E, and a larger t_E value as can be seen in Table 1. This difference is due to the new detrended MOA-R and CTIO difference imaging photometry.

Since we do not have a measurement of the microlensing parallax π_E, we use the Keck lens flux and mass–luminosity relations (Henry & McCarthy 1993; Henry et al. 1999; Delfosse et al. 2000) in order to constrain the lens distance. The extinction in the foreground of the lens is calculated assuming a dust scale height of h_dust = 0.10 ± 0.02 kpc.

The target MOA-2009-BLG-319 was observed with the NIRC2 instrument on Keck II in the H and Kshort (hereafter K) bands on 2018 May 25 and K band on 2019 May 28. The 2018 K-band data have a point-spread function (PSF) FWHM of ∼70 mas. The 2018 K-band data have somewhat poorer quality than the 2019 K-band data, and the 2018 H-band data are even more problematic, with a larger PSF (FWHM ∼ 120 mas). In Section 4.2 we discuss the analysis of the 2018 K-band data, and in Section 4.3 we test the limits of our detection capabilities with the very marginal 2018 H-band signal.

For the 2018 and 2019 observations, both the NIRC2 wide and narrow cameras were used. The pixel scales for the wide and narrow cameras are 39.69 and 9.942 mas pixel⁻¹, respectively. All of the images were taken using the Keck II laser guide star AO system.

As we discuss below in Sections 4.2 and 4.3, our highest-precision measurements come from the 2019 data, so we will focus on the analysis of those data. For the 2019 data, a co-add of nine dithered wide camera images was used for photometric calibration to images from the Vista Variables in the Via Lactea (VVV) survey (Minniti et al. 2010) following the procedure of Beaulieu et al. (2018). The wide camera images were flat-field and dark current corrected using standard methods, and stacked using the SWarp software (Bertin 2010). We performed astrometry and photometry on the co-added wide camera image using SExtractor (Bertin & Arnouts 1996) and subsequently calibrated the narrow camera images to the wide camera image by matching two dozen bright isolated stars in the frames. This calibration analysis results in uncertainties of 0.06 mag.

For the 2019 K-band narrow data, we combined 30 flat-field frames, 10 dark frames, and 15 sky frames for calibrating our science images. Following the methods of Service et al. (2016) and Yelda et al. (2010), we then combined nine K-band narrow camera science frames with an integration time of 60 s per frame. The combined frame can be seen in the left panel of Figure 2, which has a PSF FWHM of ∼73 mas. The reduction of the 2018 H-band data follows the same pipeline as the K band described in this section.

Zoom In Zoom Out Reset image size

Lastly, there were 10 K-band images of the target taken on 2015 July 26 with the NIRC2 narrow camera that were combined to make one co-added science frame. There were no sky frames taken for the 2015 data, which contributes to the lower signal-to-noise ratio seen in these data. The much smaller lens–source separation at the time of these images also implies that our lens–source relative proper motion, μ_rel,H, and lens brightness measurements will be less precise than the later images. Even with these lower-signal data, a careful DAOPHOT reduction successfully detects the lens. Further details of the 2015 analysis are given in Section 4.4. The main benefit of these early images is that they allow us to verify the identification of the lens star, by showing that it is moving away from the source at a rate consistent with the occurrence of the microlensing event in 2009 June.

4.1. PSF Fitting Photometry

Because the two stars in the blend have a separation in 2019 of ∼FWHM, it is necessary to use a PSF fitting routine to measure both targets independently. Following the methods of Bhattacharya et al. (2018) and references therein, we use the photometry routine DAOPHOT-II (Stetson 1987) to generate and fit an empirical PSF to the source+lens blend. The AO corrections for observations of our Galactic bulge fields using the instruments currently on the Keck telescope generally deliver imperfect AO corrections with Strehl ratios <0.5, and often the Strehl ratios are significantly smaller than 0.5. Thus, the PSFs delivered by the AO system can have a wide variety of shapes. The DAOPHOT package has proven to be quite successful in modeling oddly shaped PSFs delivered by the Keck AO system (Bennett et al. 2010b). An alternative method has also been presented by Vandorou et al. (2020), which is probably competitive with DAOPHOT. DAOPHOT's sophisticated semi-empirical PSF is important for our observations of MOA-2009-BLG-319 since the PSF has a prominent wing to the north that has a similar amplitude to the flux ratio of the companion star to the MOA-2009-BLG-319S source star that we interpret to be the lens star (MOA-2009-BLG-319L).

The first pass of DAOPHOT does not detect the lens, but instead produces a clear feature to the east in the residual image that can be seen in the bottom right panel (labeled "1-star res.") of Figure 2. The target is the only stellar image that has an extension in this direction, and this feature represents the position of the fainter lens star. The cyan and purple panels in Figure 2 show reference stars in the frame with similar brightness to the target that also exhibits the PSF extension to the north. This extension is accurately modeled by the DAOPHOT single-star PSF model as can be seen by the featureless residuals to the right of each reference star. The color bar on the right represents the pixel counts for the residual images only. The lens also has a separation consistent with that predicted by Miyake et al. (2011); this separation is described further in Section 5.

Fitting a two-star PSF to the target and rerunning DAOPHOT produces a nearly featureless residual, shown in the bottom right panel (labeled "2-star res.") of Figure 2. Table 2 shows the calibrated magnitudes for the two stars of K_S = 18.12 ± 0.05 and K_L = 19.98 ± 0.09. The uncertainties are derived from the "jackknife method" described in Section 4.1.2. Using the VVV extinction calculator (Gonzalez et al. 2011) and the Nishiyama et al. (2009) extinction law, we find a K-band extinction of A_K = 0.13 ± 0.05. From our reanalysis of the light-curve modeling (Section 3), we find a source color of V_L − I_L = 1.72, which leads to an extinction-corrected color of V_S0 − I_S0 = 0.80. We use the color–color relations of Kenyon & Hartmann (1995) and the I-band magnitude, I_S = 19.994, to predict a source K-band magnitude of K_S = 18.15. The fit source brightness is fainter than our measured source brightness by less than 1σ; thus, we conclude that there is virtually no evidence of additional flux from a companion to the source.

Table 2.  2019 Dual-star PSF Photometry

Star	Passband	Magnitude
Lens	Keck K	19.98 ± 0.09
Source	Keck K	18.12 ± 0.05
Source + Lens	Keck K	17.94 ± 0.06

Note. Magnitudes are calibrated to the VVV scale, as described in Section 4.

Download table as:  ASCIITypeset image

The standard version of DAOPHOT has some drawbacks for our problem of studying the closely blended images of microlens source and lens stars. First, we want to be able to study cases where the detection of the lens star may be marginal, as well as cases where we can only obtain an upper limit on the lens brightness as a function of the lens–source separation. Thus, it would be useful to have a method that will produce a probability distribution of all possible source plus lens configurations that are consistent with the data. The standard version of DAOPHOT, on the other hand, is programmed to avoid including false detections in its output star list, so it may reject some of the more marginal lens detections. Of course, because of the microlensing event, we know that another star is there, although it might be quite faint (Blackman et al. 2021). Also, as Bennett et al. (2007) have shown, constraints on the source brightness and/or lens–source separation from the light-curve models can often significantly reduce the uncertainties on parameters, such as the lens brightness, that are not significantly constrained by the light-curve data. Thus, it will be useful to be able to apply these constraints inside of DAOPHOT in order to get the most precise possible measurement of the lens star properties.

Finally, DAOPHOT does not report error bars on the star positions, which are critical for our science. King (1983) did publish a formula that can be used to estimate position error bars based on the photometry error bars, but this formula is problematic for our situation of highly blended stellar images. We can address these issues by modifying the standard version of DAOPHOT and adding a routine that uses the MCMC method to determine the distribution of source and lens star magnitudes and positions that are consistent with the data, as we explain in the next subsection.

4.1.1. Development of an MCMC Routine for DAOPHOT

We begin the MCMC routine by using the DAOPHOT empirical PSF that was described in the previous section. This PSF model is then permitted to step across the fitting box encompassing the blended targets, with a fitting box radius of ∼1.5 FWHM. For the dual-star version of the MCMC routine, there are six total parameters that are simultaneously fit: the x and y pixel location for each star (x₁, y₁, x₂, y₂), the total flux (f_T), and the flux ratio between the stars (f_R). For each step in the chain, a χ² value for the fit is measured and recorded. The routine then takes a random step in any direction (and flux), makes the same measurements, and compares the new χ² to the previous one. If the new value is smaller than the previous one, the six-parameter fit is recorded and the routine continues. However, if the new value is larger than the previous one, a weighted proposal probability distribution is calculated. If this weighted probability is less than a randomly generated probability (between 0 and 1), then the decision is reversed and the original candidate value is accepted. If the weighted probability is greater still, the candidate is rejected and the iteration moves forward with a new candidate. This procedure follows the standard Metropolis–Hastings method. Once the routine has converged, the best-fit parameters are recorded.

For a dual-star model, we calculate the flux distribution following Bhattacharya et al. (2017):

$$\begin{eqnarray}&&{f}_{T}={f}_{1}\psi (i-{x}_{1},j-{y}_{1})+(1-{f}_{1})\psi (i-{x}_{2},j-{y}_{2}),\end{eqnarray} \tag{ 2 }$$

where f₁ is the source flux contribution to the total flux, 1 − f₁ is the lens flux contribution, and ψ is the two-dimensional PSF model. The values x₁, y₁, x₂, y₂ are the initial pixel positions for the source and lens as described earlier in this section, and the indices i and j are the trial pixel positions for a given iteration. The χ² minimization routine described above computes the minimum value of the six-parameter fit as follows:

$$\begin{eqnarray}\begin{array}{rcl}{\chi }^{2} & = & \displaystyle \sum _{i,j}\left[\displaystyle \frac{1}{\sigma }\left\{{P}_{i,j}-{s}_{* }-{f}_{1}\psi (i-{x}_{1},j-{y}_{1})\right.\right.\\ & & {\left.\left.-(1-{f}_{1})\psi (i-{x}_{2},j-{y}_{2})\right\}\right]}^{2},\end{array}\end{eqnarray} \tag{ 3 }$$

where P_i,j is the intensity at pixel location i, j, σ is the uncertainty in pixel intensity, and s_* is the background flux. The MCMC chains are used as a probability distribution that we use to determine the normalized errors on the best-fit MCMC results in Tables 3 and 4.

Table 3.  DAOPHOT MCMC and Jackknife Best-fit Results

	2015 K Band		2018 K Band		2019 K Band
Parameter	MCMC	Jackknife	MCMC	Jackknife	MCMC	Jackknife
μrel,HE (mas yr−1)	6.134 ± 1.281	6.970 ± 2.187	7.172 ± 0.472	6.669 ± 0.311	6.482 ± 0.167	6.405 ± 0.072
μrel,HN (mas yr−1)	−1.351 ± 0.775	−0.555 ± 2.034	0.656 ± 0.290	0.568 ± 0.309	1.684 ± 0.158	1.788 ± 0.145
Lens flux/source flux	0.129 ± 0.069	0.158 ± 0.053	0.176 ± 0.008	0.176 ± 0.047	0.176 ± 0.007	0.180 ± 0.014

Download table as:  ASCIITypeset image

The standard version of DAOPHOT employs the Newton–Raphson method (Press et al. 1986) for fitting the positions of the two blended stars. The two-star routines were run with both the Newton–Raphson and MCMC methods, producing nearly identical results. The residual images for the reference stars shown in Figure 2 are the residuals from the Newton–Raphson analysis of standard DAOPHOT. The residual images for the target shown in the same figure are from the MCMC analysis. Figure 3 shows the 1σ, 2σ, and 3σ contour intervals for the best-fit MCMC source and lens positions, overplotted on the stellar image. The best-fit parameters from the MCMC routine and their respective error bars are listed in Table 3, along with the error bars from the jackknife method as discussed in Section 4.1.2. The lens–source separation measurement with our MCMC routine is within 1σ of the result from standard DAOPHOT.

Zoom In Zoom Out Reset image size

The routine also has the functionality to fit the simpler case of a single star. This single-star MCMC fitting was performed on the source+lens blend and produced the residual seen in Figure 2 ("1-star res."). The two-star MCMC run produces a better fit as expected, with a χ² improvement of Δχ² = 1313.0 over the single-star fit. The residual image that was created using the MCMC best-fit two-star values is nearly featureless and produced the residual shown in the bottom right panel ("2-star res.") of Figure 2.

4.1.2. Error Bars with the Jackknife Method

While the MCMC method is a powerful tool for studying the range of model parameters that are consistent with an image, there is another source of uncertainty that we must consider for our analysis of Keck AO images. It is standard practice to analyze combinations of multiple dithered infrared images in order to remove some of the instrumental artifacts from these images. However, the AO images have imperfect corrections to the optical effects of the atmosphere. The quality of the AO correction is often characterized by the Strehl ratio, which is the ratio of the brightness at the peak of a stellar PSF to the peak that would be obtained owing only to diffraction. In moderately good observing conditions, like the conditions for our 2019 K-band observations of MOA-2009-BLG-319, we typically have Strehl ratios in the range of 0.2–0.4. In H band, the Strehl ratios are worse, typically 0.1–0.2, although these images can have PSF FWHM values as good as or better than the K-band images with better Strehl ratios. Thus, greatly improved angular resolution given by these AO systems yields images that are far from perfect. Significant PSF distortions remain in the Keck AO images, and these distortions vary from image to image, and it is also likely that there is some variation across each image. Because of this, we measure the PSF with stars close to the target in our analysis, but we must also consider the effect of the variations between images.

The uncertainty due to the variations between images can be addressed by the jackknife method (Quenouille 1949, 1956; Tukey 1958). Our implementation of this method is discussed in more detail by Bhattacharya et al. (2021). To analyze a collection of N dithered images, we create N different combinations of N − 1 images, with each image missing from only one of these combinations. The error bars for each parameter are then given by $\sqrt{N-1}$ times the rms of the best-fit parameters from each of these N combinations of N − 1 images. Table 3 compares the error bars computed by the MCMC method to the error bars computed by the jackknife method. We chose to use the jackknife uncertainties because they include the uncertainties due to the PSF variations in the individual images.

4.2. 2018 K-band Analysis

4.3. 2018 H-band Analysis with Lens−Source Separation Constraint

4.4. 2015 K-band Analysis

We performed a DAOPHOT analysis of the 2015 K-band data, similar to that of the previous reductions. The PSF FWHM for this data is approximately 75 mas, which means that the lens–source separation is ∼0.53× FWHM at the time of the 2015 data approximately 6.09 yr after t₀. The lens–source relative proper motion, μ_rel, and flux ratio for the 2015 data are given in Table 4. The east and north components of the heliocentric relative proper motion from the jackknife method are consistent with both the 2018 and 2019 K-band data. Figure 4 shows the best-fit MCMC contours for the source and lens positions for each epoch, with the K-band image overplotted. The color bar refers to the pixel intensities in each frame. It is clear from these results that we are in fact measuring the lens and source moving away from one another.

Zoom In Zoom Out Reset image size

The main contribution of these 2015 images is not to increase the precision of our μ_rel,H measurements. Instead, it serves to confirm our identification of the lens star. As can be seen in Table 3, the μ_rel,H measurements from the 2015 images are consistent with the much more precise 2019 measurements. In particular, the μ_rel,HE value is within 0.25σ of the 2019 value, and the μ_rel,HN value is within 1.2σ of the 2019 value (using the jackknife error bars).

The observed motion between 2015 and 2019 rules out a possible companion to the source star as the source of the flux that we attribute to the lens star. The implied velocity is much too large for the star to be bound to the source. An unrelated star in the bulge would have to mimic the proper motion of the lens star, and the probability of this is ≲10⁻⁴ according to an analysis using the method of Koshimoto et al. (2020). There is also the possibility that we have detected the combination of the flux of the planetary host and a binary companion to the host star. The Koshimoto et al. (2020) analysis predicts a probability of 1.9% for this possibility, but this does not include a complete analysis of the triple-lens modeling for this event. There is a weak signal that could be due to an additional planet (Shin et al. 2015) or an additional star, but this will be investigated in detail in a subsequent paper.

The 2019 Keck II follow-up observations were taken 9.94 yr after peak magnification in 2009. The motion of the lens and source on the sky frame is the primary cause for their apparent separation; however, there is also a small component that can be attributed to the orbital motion of Earth. As this effect is of order ≤0.1 mas for a lens at a distance of D_L ≥ 7 kpc, we are safe to ignore this contribution in our analysis, as it is much smaller than the error bars on the stellar position measurements. The lens–source relative proper motion is measured to be μ_rel,H = (μ_rel,H,E, μ_rel,H,N) = (6.404 ± 0.072, 1.788 ± 0.145) mas yr⁻¹, where "H" indicates that these measurements were made in the heliocentric reference frame, and the "E" and "N" subscripts represent the east and north directions, respectively. Converting to Galactic coordinates, these proper motions are μ_rel,H,l = 4.670 ± 0.132 mas yr⁻¹ and μ_rel,H,b = −4.734 ± 0.095 mas yr⁻¹.

Light-curve modeling (Section 3) is most conveniently performed in the geocentric reference frame that moves with Earth at the time of the event peak. Thus, we must convert between the geocentric and heliocentric frames by using the relation given by Dong et al. (2009b):

$$\begin{eqnarray}&&{\mu }_{\mathrm{rel},{\rm{H}}}={\mu }_{\mathrm{rel},{\rm{G}}}+\displaystyle \frac{{\nu }_{\oplus }{\pi }_{\mathrm{rel}}}{\mathrm{au}},\end{eqnarray} \tag{ 4 }$$

where ν_⊕ is Earth's projected velocity relative to the Sun at the time of peak magnification. For MOA-2009-BLG-319 this value is ν_⊕E,N = (29.289, 0.347) km s⁻¹ = (6.175, 0.073) au yr⁻¹ at HJD' = 5006.99. With this information and the relative parallax relation π_rel ≡ 1/D_L − 1/D_S, we can rewrite Equation (4) in a more convenient form:

$$\begin{eqnarray}&&{\mu }_{\mathrm{rel},{\rm{G}}}={\mu }_{\mathrm{rel},{\rm{H}}}-(6.175,0.073)\times (1/{D}_{{\rm{L}}}-1/{D}_{{\rm{S}}}),\end{eqnarray} \tag{ 5 }$$

since we have directly calculated μ_rel,H from Keck. We use this relation in our Bayesian analysis of the light curve, with Galactic model and Keck constraints to determine the relative proper motion in the geocentric frame of μ_rel,G = 6.47 ± 0.12 mas. This can be compared to the value determined from the light-curve MCMC without the Keck constraints of μ_rel,G = 6.51 ± 0.59 mas, so the light-curve prediction is confirmed.

The measurement of the angular Einstein radius allows us to use a mass–distance relation if we assume that the distance to the source is known (Bennett 2008; Gaudi 2012):

$$\begin{eqnarray}&&{M}_{{\rm{L}}}=\displaystyle \frac{{c}^{2}}{4G}{\theta }_{{\rm{E}}}\displaystyle \frac{{D}_{{\rm{S}}}{D}_{{\rm{L}}}}{{D}_{{\rm{S}}}-{D}_{{\rm{L}}}},\end{eqnarray} \tag{ 6 }$$

where M_L is the lens mass, G is the gravitational constant, and c is the speed of light. D_L and D_S are the distances to the lens and source, respectively. Figure 5 shows the mass–distance plane with our new direct calculation for the lens mass and distance (black). The red curve represents the constraint from the mass–luminosity relation, with dashed lines representing the error from the Keck lens flux measurement. Additionally, the θ_E constraint is shown in green with errors dominated by the source distance uncertainty.

Zoom In Zoom Out Reset image size

As discussed in Section 3, our improved photometry and improved parameterization of Galactic bulge red clump stars yield smaller θ_*, θ_E, and μ_rel,G values. Our results from the reanalyzed light curve with detrended MOA data show a slightly fainter source star compared to Miyake et al. (2011). This yields a smaller angular Einstein radius and μ_rel that match the measured value better than the Miyake et al. (2011) value.

Table 5 shows the final planetary system results of our Bayesian analysis of the MCMC light-curve distribution constraints from our Keck observations, as well as a Galactic model. We find that the M-dwarf lens star has a mass M_L = 0.52 ± 0.05M_⊙, with a sub-Saturn planetary companion of mass m_P = 67.3 ± 6.2M_⊕. We can calculate this planet's semimajor axis using

$$\begin{eqnarray}&&{r}_{\perp }={{sD}}_{{\rm{L}}}{\theta }_{{\rm{E}}},\end{eqnarray} \tag{ 7 }$$

where s is the projected separation from the light-curve modeling; thus, we find a separation of r_⊥ = 2.03 ± 0.21 au. Additionally, the lens system is determined to be at a distance of 7.05 ± 0.71 kpc, very likely located in the Galactic bulge. Figure 6 shows the results for the physical parameters of the lens system with (red) and without (blue) the Keck constraints. The host mass and planetary mass results show very significant improvement over the unconstrained analysis, the projected separation shows marginal improvement, and the uncertainty in the lens distance is clearly still dominated by the uncertainty in the source distance, as they are highly correlated.

Zoom In Zoom Out Reset image size

Table 5.  Planetary System Properties from Lens Flux Constraints

Parameter	Units	Values and rms	2σ Range
Angular Einstein radius (θE)	mas	0.296 ± 0.006	0.283−0.309
Geocentric lens–source relative proper motion (μrel,G)	mas yr−1	6.472 ± 0.121	6.230−6.714
Host mass (Mhost)	M⊙	0.524 ± 0.048	0.428−0.621
Planet mass (Mp)	M⊕	67.3 ± 6.2	49.8−82.2
2D separation (a⊥)	au	2.03 ± 0.21	1.60−2.46
3D separation (a3D)	au	${2.90}_{-0.50}^{+1.44}$	1.88−5.78
Lens distance (DL)	kpc	7.05 ± 0.71	5.60−8.45
Source distance (DS)	kpc	8.25 ± 0.86	6.53−9.97

Download table as:  ASCIITypeset image

Our follow-up high-resolution observations of the microlensing target MOA-2009-BLG-319 have allowed us to make a direct measurement of lens flux from the host star, as well as a precise determination of the direction and amplitude of the lens–source relative proper motion. Further analysis enabled us to calculate a direct mass for the star and its planetary companion. We added a novel MCMC routine to DAOPHOT¹⁰ in order to retrieve precise astrometric and flux fits for the blended source and lens stars. It also allows constraints from the microlensing light-curve modeling to be imposed on the analysis of high angular resolution follow-up images. Following Bhattacharya et al. (2021), we performed a jackknife analysis of the Keck follow-up observations because it is able to estimate uncertainties due to variations in the Keck PSF shape in multiple images. We used these jackknife error bars for our final analysis. These methods provide more accurate results than previously used techniques for crowded field photometry in AO imaging. These routines can be used in future analyses of highly blended microlensing follow-up targets and, eventually, can form the basis for the Roman mass measurement method.

The MOA-2009-BLG-319 microlensing event has a planet-to-star mass ratio of q = (3.856 ± 0.029) × 10⁻⁴, which puts it in the range of the mass ratio desert originally predicted by Ida & Lin (2004) and confirmed more recently by Suzuki et al. (2018). This prediction was based on the runaway gas accretion scenario that has been considered a standard part of the core accretion theory (Pollack et al. 1996), but it is based on a one-dimensional calculation. The Suzuki et al. (2018) analysis found a discrepancy between the planet mass ratio distribution found by microlensing and this predicted mass ratio gap, at 10⁻⁴ < q < 4 × 10⁻⁴, thought to be caused by the rapid "runaway" growth. It was thought to be unlikely that planet growth would terminate during this predicted very rapid growth phase. However, the microlensing results of Suzuki et al. (2016) show no evidence of this predicted gap.

One possible explanation for this contradiction might be that the runaway gas accretion phase only occurs for stars of approximately solar type, which was the original focus of the core accretion theory, while microlensing probes not only solar-type stars but also lower-mass stars and even stellar remnants. Our high angular resolution follow-up observations can test this possibility by measuring host star masses for the 7 out of the 30 events in the Suzuki et al. (2016) sample that fall in the mass ratio range 10⁻⁴ < q < 4 × 10⁻⁴. Mass measurements have previously been made for two of the seven Suzuki et al. (2016) host stars with planets in this range. Bhattacharya et al. (2018) have measured a host mass of M_host = 0.58 ± 0.05M_⊙ and a planet mass of m_p = 39 ± 9M_⊕ for planetary microlensing event OGLE-2012-BLG-0950, and Bennett et al. (2016) have measured host and planet masses of M_hosts = 0.71 ± 0.12M_⊙ and m_p = 80 ± 13M_⊕ for the OGLE-2007-BLG-349L lens system, although in this case the host is a close binary pair of 0.41 ± 0.07M_⊙ and 0.31 ± 0.07M_⊙ in a ∼10-day orbit. Our group has also measured the mass of a more massive host star, OGLE-2012-BLG-0026L (Beaulieu et al. 2016), with M_host = 1.06 ± 0.05M_⊙, with one planet in the mass ratio range of the predicted gap, 10⁻⁴ < q < 4 × 10⁻⁴. The sub-Saturn planet has a mass of 46 ± 2M_⊕, and it is accompanied by a more massive planet with a mass of 265 ± 20M_⊕. However, this event is not in the Suzuki et al. (2016) statistical sample.

The addition of the MOA-2009-BLG-319L system to this collection with host and planet masses of M_host = 0.52 ± 0.05M_⊙ and m_p = 66 ± 8M_⊕ continues the trend of finding host masses within a factor of two of a solar mass, and this suggests that the lack of this mass ratio gap at 10⁻⁴ < q < 4 × 10⁻⁴ is not caused by some dramatic change in the mass ratio for host stars with very low masses. Such a conclusion would be supported by the theoretical work of J. Szulágyi et al. (2020, in preparation), who show that the runaway gas accretion phase is likely to be terminated very quickly by the formation of a circumplanetary disk, which can result in many planets in the predicted gap. Further results from our high angular resolution follow-up imaging program will provide a stronger test of these core accretion processes, with additional mass measurements for the host stars of sub-Saturn mass planets orbiting beyond the snow line. A more definitive answer to this and other questions regarding the demographics of planets in wider orbits will come from the RGES, which will have high enough angular resolution so that follow-up observations will not be needed for the majority of exoplanets discovered.

The authors thank Dr. Peter Stetson for advice on modifications to the DAOPHOT-II software. We also thank the anonymous referee for constructive comments that led to a stronger manuscript. This work was performed in part under contract with the Center for Research and Exploration in Space Sciences and Technologies (CRESST-II). The Keck observations were supported by a NASA Keck PI Data Award, 80NSSC18K0793, 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. D.P.B., A.B., and C.R. were supported by NASA through grant NASA-80NSSC18K0274. This work was supported by the University of Tasmania through the UTAS Foundation and the endowed Warren Chair in Astronomy and the ANR COLD-WORLDS (ANR-18-CE31-0002). Work by N.K. is supported by JSPS KAKENHI grant No. JP18J00897. This research was also supported in part by the Australian government through the Australian Research Council Discovery Program (project No. 200101909) grant awarded to Cole and Beaulieu. This work made use of data from the Astro Data Lab at NSF's OIR Lab, which is operated by the Association of Universities for Research in Astronomy (AURA), Inc. under a cooperative agreement with the National Science Foundation. Some of this research has made use of the NASA Exoplanet Archive, which is operated by the California Institute of Technology, under the Exoplanet Exploration Program.

Software: DAOPHOT-II (Stetson 1987), DAOPHOT-MCMC (this work), gnuplot, Matplotlib (Hunter 2007), Numpy (Oliphant 2006).