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

Sean Terry, Aparna Bhattacharya, David Bennett, Jean-Philippe Beaulieu, Naoki Koshimoto, Joshua Blackman, Ian Bond, Andrew Cole, Calen Henderson, Jessica Lu, et al.

► To cite this version:

Sean Terry, Aparna Bhattacharya, David Bennett, Jean-Philippe Beaulieu, Naoki Koshimoto, et al. MOA-2009-BLG-319Lb: A Sub-Saturn Planet inside the Predicted Mass Desert. The Astronomical Journal, 2021, 161 (2), pp.54. 10.3847/1538-3881/abcc60 . hal-03430352

HAL Id: hal-03430352

https://hal.science/hal-03430352v1

Submitted on 2 Aug 2025

**HAL** is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>
<figure><img src="image_2.png" /></figure>

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

Sean K. Terry<fnref n="1" />,<fnref n="2" /> ID, Aparna Bhattacharya<fnref n="2" />,<fnref n="3" />, David P. Bennett<fnref n="2" />,<fnref n="3" /> ID, Jean-Philippe Beaulieu<fnref n="4" />, Naoki Koshimoto<fnref n="2" />,<fnref n="3" />,<fnref n="5" /> ID,
Joshua W. Blackman<fnref n="4" /> ID, Ian A. Bond<fnref n="6" />, Andrew A. Cole<fnref n="4" /> ID, Calen B. Henderson<fnref n="7" /> ID, Jessica R. Lu<fnref n="8" /> ID, Jean Baptiste Marquette<fnref n="9" />,
Clément Ranc<fnref n="2" /> ID, and Aikaterini Vandorou<fnref n="4" /> ID

Received 2020 September 15; revised 2020 November 10; accepted 2020 November 16; published 2021 January 6

## Abstract

We present an adaptive optics (AO) analysis of images from the Keck II telescope NIRC2 instrument of the planetary microlensing event MOA-2009-BLG-319. The ~10 yr baseline between the event and the Keck observations allows the planetary host star to be detected at a separation of 66.5 ± 1.7 mas from the source star, consistent with the light-curve model prediction. The combination of the host star brightness and light-curve parameters yields host star and planet masses of $M_{\text{host}} = 0.524 \pm 0.048 M_{\odot}$ and $m_p = 67.3 \pm 6.2 M_{\oplus}$ at a distance of $D_{\text{L}} = 7.1 \pm 0.7$ kpc. The star–planet projected separation is $2.03 \pm 0.21$ au. The planet-to-star mass ratio of this system, $q = (3.857 \pm 0.029) \times 10^{-4}$, places it in the predicted “planet desert” at $10^{-4} < q < 4 \times 10^{-4}$ according to the runaway gas accretion scenario of the core accretion theory. Seven of the 30 planets in the Suzuki et al. sample fall in this mass ratio range, and this is the third with a measured host mass. All three of these host stars have masses of $0.5 \le M_{\text{host}}/M_{\odot} \le 0.7$, which implies that this predicted mass ratio gap is filled with planets that have host stars within a factor of two of $1M_{\odot}$. This suggests that runaway gas accretion does not play a major role in determining giant planet masses for stars somewhat less massive than the Sun. Our analysis has been accomplished with a modified DAOPHOT code that has been designed to measure the brightness and positions of closely blended stars. This will aid in the development of the primary method that the Nancy Grace Roman Space Telescope mission will use to determine the masses of microlens planets and their hosts.

*Unified Astronomy Thesaurus concepts: Exoplanets (498); Gravitational microlensing (672); Computational methods (1965)*

*Supporting material: data behind figure*

## 1. Introduction

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 \pm 0.02 M_{\odot}$, with a planetary companion of mass $14.1 \pm 0.9 M_{\oplus}$.

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

<footnote n="1">¹Department of Physics, Catholic University of America, 620 Michigan Ave. NE, Washington, DC 20064, USA; 41terry@cua.edu</footnote>

<footnote n="2">²Code 667, NASA Goddard Space Flight Center, Greenbelt, MD 20771, USA</footnote>

<footnote n="3">³Department of Astronomy, University of Maryland, College Park, MD 20742, USA</footnote>

<footnote n="4">⁴School of Natural Sciences, University of Tasmania, Private Bag 37 Hobart, Tasmania, 70001, Australia</footnote>

<footnote n="5">⁵Department of Astronomy, Graduate School of Science, The University of Tokyo, 7-3-1 Hongo, Bunkyo-ku, Tokyo 113-0033, Japan</footnote>

<footnote n="6">⁶Institute of Mathematical and Natural Sciences, Massey University, Auckland 0745, New Zealand</footnote>

<footnote n="7">⁷NASA Exoplanet Science Institute, IPAC/Caltech, Pasadena, CA 91125, USA</footnote>

<footnote n="8">⁸Department of Astronomy, University of California Berkeley, Berkeley, CA 94701, USA</footnote>

<footnote n="9">⁹Laboratoire d'astrophysique de Bordeaux, Univ. Bordeaux, CNRS, B18N, alle Geoffroy Saint-Hilaire, F-33615 Pessac, France</footnote>
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.

## 2. Event MOA-2009-BLG-319 and New Photometry

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.

## 3. New Light-curve Model

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 $\chi^2$ improvement, $\Delta\chi^2 > 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
<figure><img src="image_3.png" /><figcaption>**Figure 1.** Best-fit planetary light-curve model for MOA-2009-BLG-319 with the data used for the analysis in this paper. Only the sparsely sampled V-band data are not shown. The CTIO-H0 through CTIO-H4 data are treated as independent data sets, shown in different shades of green. The data behind the figure are available in machine-readable format. The data provided include the Dan, Pal, and WISE I-band measurements. All the data are presented in magnitudes units. (The data used to create this figure are available.)</figcaption></figure>

<table><caption>**Table 1**
Best-fit MOA-2009-BLG-319L Model Parameters</caption><thead><tr><th>Parameter</th><th>Units</th><th>Value</th><th>MCMC Averages</th><th>Miyake+2011</th></tr></thead><tbody><tr><td>t<sub>E</sub></td><td>days</td><td>16.762</td><td>16.72 ± 0.10</td><td>16.56 ± 0.08</td></tr><tr><td>t<sub>0</sub></td><td>HJD'</td><td>5006.9951</td><td>5006.9952 ± 0.0008</td><td>5006.995 ± 0.001</td></tr><tr><td>u<sub>0</sub></td><td></td><td>-0.006103</td><td>-0.0061 ± 0.0004</td><td>-0.0062 ± 0.0003</td></tr><tr><td>s</td><td></td><td>0.97564</td><td>0.9756 ± 0.0001</td><td>0.975 ± 0.001</td></tr><tr><td>α</td><td>rad</td><td>-2.62995</td><td>-2.6299 ± 0.0007</td><td>-2.629 ± 0.001</td></tr><tr><td>q × 10<sup>4</sup></td><td></td><td>3.8463</td><td>3.856 ± 0.029</td><td>3.95 ± 0.02</td></tr><tr><td>t<sub>*</sub></td><td>days</td><td>0.03186</td><td>0.0319 ± 0.0006</td><td>0.0320 ± 0.0033</td></tr><tr><td>I<sub>s</sub></td><td></td><td>19.994</td><td>19.992 ± 0.007</td><td>19.78 ± 0.07</td></tr><tr><td>V<sub>s</sub></td><td></td><td>21.714</td><td>21.712 ± 0.007</td><td>21.52 ± 0.09</td></tr><tr><td>χ<sup>2</sup>/dof</td><td></td><td>10746.24/10805</td><td></td><td></td></tr></tbody></table>

**Note.** HJD' = HJD - 2,450,000. Miyake et al. (2011) values are for their best-fit $u_0 < 0$ solution without parallax. We have performed a change of coordinate for $\alpha$ reported in Miyake et al. (2011) by $\pi \to \pi - \alpha$, based on the choice of “mass one” for the planet.

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 \pm 0.069$ and $V_{s0} = 19.678 \pm 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:

$$ \log(2\theta_*) = 0.5014 + 0.4197(V_{s0} - I_{s0}) - 0.2I_{s0}. \quad (1) $$

This yields $\theta_* = 0.576 \pm 0.077$ µas, which is smaller than the Miyake et al. (2011) value of $\theta_* = 0.66 \pm 0.06$ µas. Our
measurement is consistent with the $\mu_{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 $\mu_{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 \pm 0.86$ kpc. These new light-curve modeling results produce smaller best-fit values for the mass ratio, $q$, and angular Einstein radius, $\theta_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 $\pi_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 \pm 0.02$ kpc.

# 4. Keck Follow-up and Analysis

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$^{-1}$, 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.

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, $\mu_{rel,H}$, and lens brightness measurements will be less precise than the later images. Even with these lower-signal data, a careful DAO-PHOT 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 \pm 0.05$ and $K_L = 19.98 \pm 0.09$. The uncertainties are derived from the "jackknife method" described in Section 4.1.2. Using the VVV extinction calculator (Gonzalez
<figure><img src="image_4.png" /><figcaption>**Figure 2.** Left panel: co-added sum of nine 60 s NIRC2 K-band narrow camera images from 2019. Cyan and purple panels: close-up of single stars in the frame, with one-star PSF residuals plotted next to each. Red panel: close-up of MOA-2009-BLG-319 showing center position of the source (red point) and lens (yellow point), with one-star and two-star PSF residuals, respectively. The color bar refers to the PSF residual images only.</figcaption></figure>

<table><caption>**Table 2**
2019 Dual-star PSF Photometry</caption><thead><tr><td>Star</td><td>Passband</td><td>Magnitude</td></tr></thead><tbody><tr><td>Lens</td><td>Keck K</td><td>19.98 ± 0.09</td></tr><tr><td>Source</td><td>Keck K</td><td>18.12 ± 0.05</td></tr><tr><td>Source + Lens</td><td>Keck K</td><td>17.94 ± 0.06</td></tr></tbody></table>

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

et al. 2011) and the Nishiyama et al. (2009) extinction law, we find a K-band extinction of $A_K = 0.13 \pm 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\sigma$; thus, we conclude that there is virtually no evidence of additional flux from a companion to the source.

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_1, y_1, x_2, y_2$), the total flux ($f_T$), and the flux ratio between the stars ($f_R$). For each step in the chain, a $\chi^2$ 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 $\chi^2$ 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
<table><caption>**Table 3**
DAOPHOT MCMC and Jackknife Best-fit Results</caption><thead><tr><th rowspan="2">Parameter</th><th colspan="2">2015 K Band</th><th colspan="2">2018 K Band</th><th colspan="2">2019 K Band</th></tr><tr><th>MCMC</th><th>Jackknife</th><th>MCMC</th><th>Jackknife</th><th>MCMC</th><th>Jackknife</th></tr></thead><tbody><tr><td>μ<sub>rel,HE</sub> (mas yr<sup>-1</sup>)</td><td>6.134 ± 1.281</td><td>6.970 ± 2.187</td><td>7.172 ± 0.472</td><td>6.669 ± 0.311</td><td>6.482 ± 0.167</td><td>6.405 ± 0.072</td></tr><tr><td>μ<sub>rel,HN</sub> (mas yr<sup>-1</sup>)</td><td>-1.351 ± 0.775</td><td>-0.555 ± 2.034</td><td>0.656 ± 0.290</td><td>0.568 ± 0.309</td><td>1.684 ± 0.158</td><td>1.788 ± 0.145</td></tr><tr><td>Lens flux/source flux</td><td>0.129 ± 0.069</td><td>0.158 ± 0.053</td><td>0.176 ± 0.008</td><td>0.176 ± 0.047</td><td>0.176 ± 0.007</td><td>0.180 ± 0.014</td></tr></tbody></table>

<table><caption>**Table 4**
Best-fit MCMC Results for Relative Proper Motion and Flux Ratio</caption><thead><tr><th>Parameter</th><th>2015 K</th><th>2018 H</th><th>2018 K</th><th>2019 K</th></tr></thead><tbody><tr><td>μ<sub>rel,HE</sub> (mas yr<sup>-1</sup>)</td><td>6.134 ± 1.281</td><td>6.183 ± 0.449</td><td>7.172 ± 0.472</td><td>6.482 ± 0.167</td></tr><tr><td>μ<sub>rel,HN</sub> (mas yr<sup>-1</sup>)</td><td>-1.351 ± 0.775</td><td>1.823 ± 0.889</td><td>0.656 ± 0.290</td><td>1.684 ± 0.158</td></tr><tr><td>Lens flux/source flux</td><td>0.129 ± 0.069</td><td>0.034 ± 0.009</td><td>0.176 ± 0.008</td><td>0.175 ± 0.007</td></tr></tbody></table>

**Note.** The 2018 *H*-band lens–source flux ratio is unreliable, as described in Section 4.3, and we regard the small flux ratio MCMC error as significantly underestimated.

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):

$$f_T = f_1 \psi(i - x_1, j - y_1) + (1 - f_1) \psi(i - x_2, j - y_2), \quad (2)$$

where $f_1$ is the source flux contribution to the total flux, $1 - f_1$ is the lens flux contribution, and $\psi$ is the two-dimensional PSF model. The values $x_1, y_1, x_2, y_2$ 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 $\chi^2$ minimization routine described above computes the minimum value of the six-parameter fit as follows:

$$\chi^2 = \sum_{i,j} \left[ \frac{1}{\sigma} \left\{ P_{i,j} - s_* - f_1 \psi(i - x_1, j - y_1) - (1 - f_1) \psi(i - x_2, j - y_2) \right\} \right]^2, \quad (3)$$

where $P_{i,j}$ is the intensity at pixel location *i*, *j*, $\sigma$ 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.

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\sigma$, $2\sigma$, and $3\sigma$ 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\sigma$ of the result from standard DAOPHOT.

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 $\chi^2$ improvement of $\Delta\chi^2 = 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
<figure><img src="image_5.png" /><figcaption>**Figure 3.** Best-fit MCMC contours (68.3%, 99.5%, 99.7%) for the source and lens positions, respectively, overplotted on the *K*-band image of the target. The lens contributes ~15% of flux to the total blend.</figcaption></figure>

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

In addition to the 2019 *K*-band data discussed in detail above, we also obtained a set of 13 30 s exposure NIRC2 narrow camera images on 2018 May 25. A total of 20 calibration frames were used for flat-fielding, dark subtraction, and sky subtraction.

The 2018 *K*-band images have a PSF FWHM similar to the 2019 *K*-band images, although the PSF appears to be slightly elongated in the east-west direction instead of having the extended wing to the north, like the 2019 *K*-band images. This is a complication because the lens star is located toward the east, but the more serious issue is that these images are much noisier. They have been taken through ~0.7 mag of extinction due to cirrus clouds, and there appears to have been a substantial amount of moonlight reflected off the clouds. This generated a much higher background and probably prevented the sky subtraction from removing some systematic errors.

We reduced these data with the same procedures used for the 2019 *K*-band data described above, and the results were very similar to the 2019 *K*-band results. However, as shown in Table 3, the error bars from the jackknife method were significantly larger than for the 2019 *K*-band data, particularly for the lens/source flux ratio and the $\mu_{rel,HE}$ component of the relative proper motion. Therefore, we use the 2019 *K*-band data for our constraints on the properties of the lens system, although the results with the weighted sum of the 2018 and 2019 *K*-band data are indistinguishable.

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

The clear lens detection from the 2019 *K*-band data allows us to carefully test the capabilities of our observations and analysis on a data set with a marginal detection (i.e., the *H*-band data for MOA-2009-BLG-319). Our initial reduction of this data with standard DAOPHOT did not detect the lens or show an obvious feature in the best-fit single-star residual to indicate the presence of the lens star. In addition, our first attempts at two-star fits with the MCMC version of DAOPHOT also did not successfully converge on a lens location. Following the methods described in Bhattacharya et al. (2017), we implemented a separation constraint to our MCMC analysis based on the known $\mu_{rel}$ from our light-curve reanalysis. While we could also constrain the 2018 separation based on our 2019 *K*-band lens-source separation measurement, our goal is to show the reliability of a marginal detection with MCMC on future targets that do not have any such better data. With this lens-source separation constraint, along with a renormalization of the pixel errors such that the best-fit $\chi^2$/dof $\simeq 1$, the MCMC converged on a solution for the lens location of 57.5 $\pm$ 2.4 mas to the NE of the source, consistent with the 2019 data. The renormalization factor for our *H*-band analysis was 0.256, and the total number of fitted pixels was 2304. Finally, we test the stability of the PSF model by calculating the total $\chi^2$ of the pixels from a radius of 1 pixel from the center of the bright source to a radius of the size of the fitting box. We find a relatively smooth distribution in $\chi^2$/pixel space, which indicates a stable PSF model.

We subsequently reran the MCMC routine with the separation constraint and renormalized errors, and our best-fit results show that the lens is detected, albeit with less confidence than the *K*-band result. The best-fit results for the *H* band are shown in Table 4. One drawback we find during this marginal detection test is that the best-fit lens-source flux ratio is not consistent with the 2019 result. The contrast should be somewhat lower in *H* band since the lens is redder than the source; however, the *H*-band results are more than 10$\sigma$ lower than the *K*-band results.

## 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 $\sim$0.53× FWHM at the time of the 2015 data approximately 6.09 yr after $t_0$. The lens-source relative proper motion, $\mu_{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.

The main contribution of these 2015 images is not to increase the precision of our $\mu_{rel,H}$ measurements. Instead, it serves to confirm our identification of the lens star. As can be seen in Table 3, the $\mu_{rel,H}$ measurements from the 2015 images are consistent with the much more precise 2019 measurements. In particular, the $\mu_{rel,HE}$ value is within 0.25$\sigma$ of the 2019
<figure><img src="image_6.png" /><figcaption>Figure 4. Best-fit MCMC contours (68.3%, 99.5%, 99.7%) for the source and lens positions overplotted on the 0.3" × 0.3" K-band images from 2015 (left), 2018 (middle), and 2019 (right). The color bar refers to the pixel intensity. North is up and east is left in all panels. This series of data clearly show that the lens and source are separating from each other. While the MCMC calculations provide enough resolution to calculate contours, they can often be underestimated because they exclude any effects of PSF variations between images.</figcaption></figure>

value, and the $\mu_{rel,HN}$ value is within $1.2\sigma$ 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 $\lesssim 10^{-4}$ 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.

## 5. Lens–Source Relative Proper Motion

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 $\le 0.1$ mas for a lens at a distance of $D_L \ge 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 $\mu_{rel,H} = (\mu_{rel,H,E}, \mu_{rel,H,N}) = (6.404 \pm 0.072, 1.788 \pm 0.145)$ mas yr$^{-1}$, 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 $\mu_{rel,H,1} = 4.670 \pm 0.132$ mas yr$^{-1}$ and $\mu_{rel,H,2} = -4.734 \pm 0.095$ mas yr$^{-1}$.

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):

$$ \mu_{rel,H} = \mu_{rel,G} + \frac{\nu_{\oplus}\pi_{rel}}{au}, \quad (4) $$

where $\nu_{\oplus}$ is Earth’s projected velocity relative to the Sun at the time of peak magnification. For MOA-2009-BLG-319 this value is $\nu_{\oplus,E,N} = (29.289, 0.347)$ km s$^{-1} = (6.175, 0.073)$ au yr$^{-1}$ at HJD' = 5006.99. With this information and the relative parallax relation $\pi_{rel} = 1/D_L - 1/D_S$, we can rewrite Equation (4) in a more convenient form:

$$ \mu_{rel,G} = \mu_{rel,H} - (6.175, 0.073) \times (1/D_L - 1/D_S), \quad (5) $$

since we have directly calculated $\mu_{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 $\mu_{rel,G} = 6.47 \pm 0.12$ mas. This can be compared to the value determined from the light-curve MCMC without the Keck constraints of $\mu_{rel,G} = 6.51 \pm 0.59$ mas, so the light-curve prediction is confirmed.

## 6. Lens System Properties

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; Gaudí 2012):

$$ M_L = \frac{c^2}{4G} \theta_E \frac{D_S D_L}{D_S - D_L}, \quad (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 $\theta_E$ constraint is shown in green with errors dominated by the source distance uncertainty.

As discussed in Section 3, our improved photometry and improved parameterization of Galactic bulge red clump stars yield smaller $\theta_*$, $\theta_E$, and $\mu_{rel,G}$ values. Our results from the
<figure><img src="image_7.png" /><figcaption>**Figure 5.** Mass–distance relation for MOA-2009-BLG-319 with constraints from the K-band lens flux measurement (red curve) and angular Einstein radius measurement (green curve).</figcaption></figure>

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 $\mu_{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 \pm 0.05 M_\odot$, with a sub-Saturn planetary companion of mass $m_p = 67.3 \pm 6.2 M_\odot$. We can calculate this planet's semimajor axis using

$$r_{\perp} = sD_L\theta_E, \quad (7)$$

where *s* is the projected separation from the light-curve modeling; thus, we find a separation of $r_{\perp} = 2.03 \pm 0.21$ au. Additionally, the lens system is determined to be at a distance of $7.05 \pm 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.

## 7. Discussion and Conclusion

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 \pm 0.029) \times 10^{-4}$, 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^{-4} < q < 4 \times 10^{-4}$, 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^{-4} < q < 4 \times 10^{-4}$. 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 \pm 0.05 M_\odot$ and a planet mass of $m_p = 39 \pm 9 M_\odot$ for planetary microlensing event OGLE-2012-BLG-0950, and Bennett et al. (2016) have measured host and planet masses of $M_{hosts} = 0.71 \pm 0.12 M_\odot$ and

<table><caption>**Table 5**
Planetary System Properties from Lens Flux Constraints</caption><thead><tr><th>Parameter</th><th>Units</th><th>Values and rms</th><th>2σ Range</th></tr></thead><tbody><tr><td>Angular Einstein radius (θ<sub>E</sub>)</td><td>mas</td><td>0.296 ± 0.006</td><td>0.283–0.309</td></tr><tr><td>Geocentric lens–source relative proper motion (μ<sub>rel,G</sub>)</td><td>mas yr<sup>-1</sup></td><td>6.472 ± 0.121</td><td>6.230–6.714</td></tr><tr><td>Host mass (M<sub>host</sub>)</td><td>M<sub>☉</sub></td><td>0.524 ± 0.048</td><td>0.428–0.621</td></tr><tr><td>Planet mass (M<sub>p</sub>)</td><td>M<sub>⊕</sub></td><td>67.3 ± 6.2</td><td>49.8–82.2</td></tr><tr><td>2D separation (a<sub>⊥</sub>)</td><td>au</td><td>2.03 ± 0.21</td><td>1.60–2.46</td></tr><tr><td>3D separation (a<sub>3D</sub>)</td><td>au</td><td>2.90<sup>+1.44</sup><sub>-0.50</sub></td><td>1.88–5.78</td></tr><tr><td>Lens distance (D<sub>L</sub>)</td><td>kpc</td><td>7.05 ± 0.71</td><td>5.60–8.45</td></tr><tr><td>Source distance (D<sub>S</sub>)</td><td>kpc</td><td>8.25 ± 0.86</td><td>6.53–9.97</td></tr></tbody></table>
<figure><img src="image_8.png" /><figcaption>**Figure 6.** Bayesian posterior probability distributions for the planetary companion mass, host mass, their separation, and the distance to the lens system are shown with only light-curve constraints in blue and with the additional constraints from our Keck follow-up observations in red. The central 68.3% of the distributions are shaded in darker colors (dark red and dark blue), and the remaining central 95.4% of the distributions are shaded in lighter colors. The vertical black line marks the median of the probability distribution for the respective parameters.</figcaption></figure>

$m_p = 80 \pm 13M_⊕$ for the OGLE-2007-BLG-349L lens system, although in this case the host is a close binary pair of $0.41 \pm 0.07M_⊕$ and $0.31 \pm 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_{\text{host}} = 1.06 \pm 0.05M_⊕$, with one planet in the mass ratio range of the predicted gap, $10^{-4} < q < 4 \times 10^{-4}$. The sub-Saturn planet has a mass of $46 \pm 2M_⊕$, and it is accompanied by a more massive planet with a mass of $265 \pm 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_{\text{host}} = 0.52 \pm 0.05M_⊕$ and $m_p = 66 \pm 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^{-4} < q < 4 \times 10^{-4}$ 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

<footnote n="10">¹⁰ https://github.com/skterry/DAOPHOT-MCMC</footnote>
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).

## ORCID iDs

Sean K. Terry https://orcid.org/0000-0002-5029-3257

David P. Bennett https://orcid.org/0000-0001-8043-8413

Naoki Koshimoto https://orcid.org/0000-0003-2302-9562

Joshua W. Blackman https://orcid.org/0000-0001-5860-1157

Andrew A. Cole https://orcid.org/0000-0003-0303-3855

Calen B. Henderson https://orcid.org/0000-0001-8877-9060

Jessica R. Lu https://orcid.org/0000-0001-9611-0009

Clément Ranc https://orcid.org/0000-0003-2388-4534

Aikaterini Vandorou https://orcid.org/0000-0002-9881-4760

## References

Batista, V., Beaulieu, J. P., Bennett, D. P., et al. 2015, ApJ, 808, 170

Beaulieu, J. P., Batista, V., Bennett, D. P., et al. 2018, AJ, 155, 78

Beaulieu, J. P., Bennett, D. P., Batista, V., et al. 2016, ApJ, 824, 83

Bennett, D. 2018, Development of the WFIRST Exoplanet Mass Measurement Method, Keck Observatory Archive N021

Bennett, D. P. 2008, in Exoplanets, ed. J. W. Mason (Berlin: Springer), 47

Bennett, D. P. 2010, ApJ, 716, 1408

Bennett, D. P., Anderson, J., Beaulieu, J. P., et al. 2010a, arXiv:1012.4486

Bennett, D. P., Anderson, J., Bond, I. A., Udalski, A., & Gould, A. 2006, ApJL, 647, L171

Bennett, D. P., Anderson, J., & Gaudi, B. S. 2007, ApJ, 660, 781

Bennett, D. P., Batista, V., Bond, I. A., et al. 2014, ApJ, 785, 155

Bennett, D. P., Bhattacharya, A., Anderson, J., et al. 2015, ApJ, 808, 169

Bennett, D. P., Bhattacharya, A., Beaulieu, J.-P., et al. 2020, AJ, 159, 68

Bennett, D. P., Bond, I. A., Abe, F., et al. 2017, AJ, 154, 68

Bennett, D. P., & Rhie, S. H. 1996, ApJ, 472, 660

Bennett, D. P., & Rhie, S. H. 2002, ApJ, 574, 985

Bennett, D. P., Rhie, S. H., Nikolaev, S., et al. 2010b, ApJ, 713, 837

Bennett, D. P., Rhie, S. H., Udalski, A., et al. 2016, AJ, 152, 125

Bennett, D. P., Sumi, T., Bond, I. A., et al. 2012, ApJ, 757, 119

Bertin, E. 2010, SWarp: Resampling and Co-adding FITS Images Together, Astrophysics Source Code Library, ascl:1010.068

Bertin, E., & Amouts, S. 1996, A&AS, 117, 393

Bhattacharya, A., Beaulieu, J. P., Bennett, D. P., et al. 2018, AJ, 156, 289

Bhattacharya, A., Bennett, D. P., Anderson, J., et al. 2017, AJ, 154, 59

Bhattacharya, A., Bennett, D. P., Beaulieu, J. P., et al. 2021, AJ, submitted (arXiv:2009.02329)

Bhattacharya, A., Bennett, D. P., Bond, I. A., et al. 2016, AJ, 152, 140

Blackman, J. W., Beaulieu, J. P., & Bennett, D. P. 2021, NatAs, submitted

Bond, I. A., Abe, F., Dodd, R. J., et al. 2001, MNRAS, 327, 868

Bond, I. A., Bennett, D. P., Sumi, T., et al. 2017, MNRAS, 469, 2434

Borucki, W. J., Koch, D. G., Basri, G., et al. 2011, ApJ, 736, 19

Boyajian, T. S., van Belle, G., & von Braun, K. 2014, AJ, 147, 47

Delfosse, X., Forveille, T., Ségransan, D., et al. 2000, A&A, 364, 217

Dong, S., Bond, I. A., Gould, A., et al. 2009a, ApJ, 698, 1826

Dong, S., Gould, A., Udalski, A., et al. 2009b, ApJ, 695, 970

Gaudi, B. S. 2012, ARA&A, 50, 411

Gonzalez, O. A., Rejkuba, M., Minniti, D., et al. 2011, A&A, 534, L14

Gould, A., & Loeb, A. 1992, ApJ, 396, 104

Gould, A., Udalski, A., Monard, B., et al. 2009, ApJL, 698, L147

Gould, A., Udalski, A., Shin, I.-G., et al. 2014, Sci, 345, 46

Hardy, S. J., & Walker, M. A. 1995, MNRAS, 276, L79

Henry, T. J., Franz, O. G., Wasserman, L. H., et al. 1999, ApJ, 512, 864

Henry, T. J., & McCarthy, D. W. J. 1993, AJ, 106, 773

Holz, D. E., & Wald, R. M. 1996, ApJ, 471, 64

Hunter, J. D. 2007, CSE, 9, 90

Ida, S., & Lin, D. N. C. 2004, ApJ, 604, 388

Johnson, S. A., Penny, M., Gaudi, B. S., et al. 2020, ApJ, 160, 123

Kenyon, S. J., & Hartmann, L. 1995, ApJS, 101, 117

King, I. R. 1983, PASP, 95, 163

Koshimoto, N., Bennett, D. P., & Suzuki, D. 2020, AJ, 159, 268

Lissauer, J. J. 1993, ARA&A, 31, 129

Mao, S., & Paczynski, B. 1991, ApJL, 374, L37

Minniti, D., Lucas, P. W., Emerson, J. P., et al. 2010, NewA, 15, 433

Miyake, N., Sumi, T., Dong, S., et al. 2011, ApJ, 728, 120

Nataf, D. M., Gould, A., Fouqué, P., et al. 2013, ApJ, 769, 88

Nishiyama, S., Tamura, M., Hatano, H., et al. 2009, ApJ, 696, 1407

Oliphant, T. E. 2006, A Guide to NumPy, Vol. 1 (Spanish Fork, UT: Trelgol Publishing)

Penny, M. T., Gaudi, B. S., Kerins, E., et al. 2019, ApJS, 241, 3

Pollack, J. B., Hubickyj, O., Bodenheimer, P., et al. 1996, Icar, 124, 62

Press, W. H., Flannery, B. P., Teukolsky, S. A., & Vetterling, W. T. 1986, Numerical Recipes. The Art of Scientific Computing (Cambridge: Cambridge Univ. Press), 839

Quenouille, M. H. 1949, Ann. Math. Stat., 20, 355

Quenouille, M. H. 1956, Biometrika, 43, 353

Schechter, P. L., Mateo, M., & Saha, A. 1993, PASP, 105, 1342

Service, M., Lu, J. R., Campbell, R., et al. 2016, PASP, 128, 095004

Shin, I. G., Han, C., Choi, J. Y., et al. 2015, ApJ, 802, 108

Spergel, D., Gehrels, N., Baltay, C., et al. 2015, arXiv:1503.03757

Stetson, P. B. 1987, PASP, 99, 191

Sumi, T., Abe, F., Bond, I. A., et al. 2003, ApJ, 591, 204

Suzuki, D., Bennett, D. P., Ida, S., et al. 2018, ApJL, 869, L34

Suzuki, D., Bennett, D. P., Sumi, T., et al. 2016, ApJ, 833, 145

Szymański, M. K., Udalski, A., Soszyński, I., et al. 2011, AcA, 61, 83

Tukey, J. W. 1958, Ann. Math. Stat., 29, 614

Vandorou, A., Bennett, D. P., Beaulieu, J.-P., et al. 2020, AJ, 160, 121

Yelda, S., Lu, J. R., Ghez, A. M., et al. 2010, ApJ, 725, 331