MOA-2020-BLG-208Lb: Cool Sub-Saturn Planet Within Predicted Desert

We analyze the MOA-2020-BLG-208 gravitational microlensing event and present the discovery and characterization of a new planet, MOA-2020-BLG-208Lb, with an estimated sub-Saturn mass. With a mass ratio $q = 3.17^{+0.28}_{-0.26} \times 10^{-4}$ and a separation $s = 1.3807^{+0.0018}_{-0.0018}$, the planet lies near the peak of the mass-ratio function derived by the MOA collaboration (Suzuki et al. 2016), near the edge of expected sample sensitivity. For these estimates we provide results using two mass law priors: one assuming that all stars have an equal planet-hosting probability, and the other assuming that planets are more likely to orbit around more massive stars. In the first scenario, we estimate that the lens system is likely to be a planet of mass $m_\mathrm{planet} = 46^{+42}_{-24} \; M_\oplus$ and a host star of mass $M_\mathrm{host} = 0.43^{+0.39}_{-0.23} \; M_\odot$, located at a distance $D_L = 7.49^{+0.99}_{-1.13} \; \mathrm{kpc}$. For the second scenario, we estimate $m_\mathrm{planet} = 69^{+37}_{-34} \; M_\oplus$, $M_\mathrm{host} = 0.66^{+0.35}_{-0.32} \; M_\odot$, and $D_L = 7.81^{+0.93}_{-0.93} \; \mathrm{kpc}$. As a cool sub-Saturn-mass planet, this planet adds to a growing collection of evidence for revised planetary formation models and qualifies for inclusion in the extended MOA-II exoplanet microlensing sample.


INTRODUCTION
Gravitational microlensing (Mao & Paczynski 1991) provides a means for detecting planets that is sensitive to low-mass planets orbiting at moderate to large distances from their host star (Bennett & Rhie 1996;Gould & Loeb 1992), typically from 0.5-10 AU. Such planets may be challenging to detect via other common exoplanet detection methods (e.g., photometric transits), hence gravitational microlensing helps provide a more complete understanding of planet statistics by providing access to another population of planets (Bennett 2009;Gaudi 2012).
The first planetary microlensing event was discovered by Bond et al. (2004). Expected to launch in 2026, the Nancy Grace Roman Space Telescope (Roman), a NASA flagship mission, will survey 10 8 stars for microlensing events (Penny et al. 2019). With less than 200 planets discovered via microlensing thus far, each new planetary microlensing analysis facilitates the calibration of theory and influences the science goals and operational plan of large-scale missions such as Roman.
A recent statistical analysis of planetary signals discovered using gravitational microlensing implied that cold, Neptune-mass planets are likely to be the most common type of planets beyond the snow line (Suzuki et al. 2016;Jung et al. 2019). This was inferred from a break in the planet-to-host-star mass-ratio function for a mass-ratio q ∼ 10 −4 , with the break resulting in a peak at Neptune-mass planets. Although the Suzuki et al. (2016) sample generally supports the planet distribution predictions from core accretion theory population synthesis models for planets beyond the snow line (Ida & Lin 2004;Mordasini et al. 2009), the existence of the peak in the planet-to-host-star mass-ratio distribution at Neptune-mass planets in the sample distribution presents issues. Specifically, these models of the planet distribution predict a dearth of sub-Saturn-mass planets, which conflicts with the microlensing observations (Suzuki et al. 2018). However, the Suzuki et al. (2016) sample consisted of only 30 exoplanets. This small sample size combined with the apparent contradiction emphasizes the importance of additional analyses, such as the one presented in this work. Formation models that propose a shortage of cold sub-Saturnmass planets (Ida & Lin 2004;Mordasini et al. 2009;Ida et al. 2013) would be contested by such a population of planets (Suzuki et al. 2018;Zang et al. 2020;Terry et al. 2021), and revised formation models would be required (Ali-Dib et al. 2022).

OBSERVATIONS AND DATA REDUCTION
The microlensing event MOA-2020-BLG-208 was discovered by the Microlensing Observations in Astrophysics collaboration (MOA) and first alerted on 2020 August 11.
The event was located at  (Bond et al. 2001) based on the difference imaging method of Tomaney & Crotts (1996). A re-reduction was carried out using the method detailed in Bond et al. (2017) resulting in photometry calibrated to phase-III of the Optical Gravitational Lensing Experiment (OGLE-III, Szymański et al. 2011).
Observations of the event, particularly of the primary lens event, were made by several other collaborations.
At UT 16:45 on 6 September 2020 (HJD ′ = 9099.2), the MAP Follow-Up Collaboration and the Microlensing Follow Up Network (µFUN) found that this event could become a high magnification within two days based on the real-time MOA data and thus conducted follow-up observations. Their follow-up observations were taken using the 1.0m telescope of the Las Cumbres Observatory (LCO) global network (Brown et al. 2013) located at the South African Astronomical Observatory, South Africa (LCOS), the 0.4m telescopes at Auckland Observatory (Auck) and Possum Observatory (POS), the 0.36m telescopes at Kumeu Observatory (Kumeu), Klein Karoo Observatory (KKO), Turitea Observatory (Turitea) and Farm Cove Observatory (FCO), and the 0.3m Perth Exoplanet Survey Telescope (PEST) at Australia. The LCOS data were reduced using a custom pipeline based on the ISIS package (Alard & Lupton 1998;Alard 2000;Zang et al. 2018), and the µFUN data were reduced using DoPHOT (Schechter et al. 1993).
In addition, this event was observed by the Korea Microlensing Telescope Network (Kim et al. 2016) with two 1.6m telescopes located at the Siding Spring Observatory, Australia and the South African Astronomical Observatory, South Africa (the site at CTIO in Chile was closed due to covid).
The Observing Microlensing Events of the Galaxy Automatically project (OMEGA) observed the event with 1.0m telescopes located at the South African Astronomical Observatory, South Africa, using the Las Cumbres Observatory network of robotic telescopes (Brown et al. 2013). The OMEGA data includes SDSS-i' and SDSS-g' Figure 1. Our best-fit planetary model for a wide-orbit solution for the MOA-2020-BLG-208 event light curve, shown in black. The plot shows magnification over time (HJD'). The lower panel shows the residual of the observations from the fit model. Unfilled observation points are observations that were excluded from the fitting due to poor seeing or a χ 2 cut. See text for instrument details. Zoomed inset figures show the peak of the primary magnification, the anomaly, and a portion of light curve with a significant fitting difference compared to the close-orbit solution (see Figure 2). The best-fit single-source single-lens model is shown in gray.
bands, with data reduced on a filter basis and uses the pyDanDIA photometry pipeline (pyDanDIA Contributors 2017).

LIGHT CURVE MODEL
The primary lens peak of the MOA-2020-BLG-208 event light curve (see Figure 1) generally resembles a Paczýnski curve (Paczynski 1986), which is the expected shape for a microlensing event with a single lens. The deviation from the Paczýnski curve occurs in the form of a secondary anomaly near HJD ′ = 9113.6 1 . The secondary anomaly in concert with the primary event is suggestive of a binary lens system comprised of a star and a companion. To model the distribution of likely properties (e.g., mass-ratio, orbital separation) that de-1 HJD' = HJD − 2, 450, 000 fine this binary lens system, we apply the method described by Bennett (2010). To perform this modeling and analysis process, we include photometric measurements from the 15 instrument data sets that are listed in Section 2.

Model fitting parameters
The approach presented by Bennett (2010) uses an image-centered, ray-shooting method (Bennett & Rhie 1996) combined with the Metropolis-Hastings algorithm (Metropolis et al. 1953;Hastings 1970), which convergences to find χ 2 minima. The parameters of our modeling consist of: the Einstein crossing time (t E ); the time that the minimum separation of lens and source occurs (t 0 ); the minimum separation between source and lens as seen by the observer (u 0 ); the separation of the two masses of the binary lens system during the event (s); the counter-clockwise angle between the lens-source relative motion projected onto the sky plane and the binary lens axis (θ); the mass-ratio between the secondary lens and the primary lens (q); the source radius crossing time (t * ); two values to model parallax in polar coordinates (r πE and θ πE ); the source flux for each instrument i (f s,i ); and the blend flux per instrument i (f b,i ).
The parameters t E , t 0 and u 0 are the common parameters for the single-lens model, while s, θ and q are the additional parameters for a binary lens system model. Both length parameters, u 0 and s, are normalized by the angular Einstein radius θ E , defined by where G is the gravitational constant, M L is the total mass of the lens system, c is the speed of light, D S is the observer-source distance, and D L is the observerlens distance. The source radius crossing time, t * , is a parameter used to take into account finite source effects, where ρ is the source angular radius in Einstein units, and θ * is the source angular radius. The microlens parallax is denoted π E , with where π rel is the lens-source relative parallax (Gould et al. 2004;Gould 1992). Our model fits the parallax parameters in polar coordinates. The equivalent Cartesian coordinates are given by x π E = r πE * cos(θ πE ) and y π E = r πE * sin(θ πE ).
To account for blending with nearby unlensed stars and different photometric systems, we normalize the flux data from each data source independently to minimize the χ 2 . As the observed brightness is linearly dependent on the blend and source fluxes, all fluxes from a single data source are normalized together to determine the normalization that produces the minimal χ 2 . The total flux F i (t) for time t for passband i is given by: where A(t, x) is the magnification of the event at time t for a set of lens model parameters x = (t E , t 0 , u 0 , s, q, t * ), f s,i is the unlensed source flux for passband i, and f b,i is the blended flux for passband i. Each instrument's passbands are independently normalized, as the relative scales of the source and blend fluxes are dependent on the instrument and the method used for difference imaging. For example, the method presented by Bond et al. (2017) is used to process the MOA data, which normalizes the target flux to match the flux of the nearest star-like object in the reference frame.

Fitting procedure
Our fitting procedure begins with a manual rough estimate selection of t E , t 0 , u 0 , and t * . With t E , t 0 , u 0 , and t * fixed at these selected values, we run a grid search of model fits varying s, q, and θ as described in Bennett (2010). We then select the best-fit set of parameters from this grid search for local minima (which include models from both wide-orbit and close-orbit solutions). We repeat the remaining steps for each minima model explored.
We use an MCMC fitting method (Bennett 2010) based on the Metropolis-Hastings algorithm. Our process includes an initial simulated annealing fitting, a renormalization, a second simulated annealing fitting on the renormalized values, an initial MCMC run to determine the covariance of the parameters, and finally an MCMC run with the covariance. We use the previously determined best fit grid search model parameters for s, q, and θ, the estimated values for t E , t 0 , u 0 , and t * , and a value of 0 for r πE and θ πE as the initial state of the MCMC algorithm. The MCMC varies all parameters during the fitting process. The acceptance of a step in the Metropolis-Hastings algorithm in our case is based on the χ 2 of the fit of the model to the data. Final reported χ 2 statistics use the renormalization of the overall best-fit model. This model produced the best χ 2 result regardless of which renormalization was applied.

LIGHT CURVE MODELING RESULTS
In Figure 1, we show the best-fit model of our light curve fitting. This model is a single source binary lens wide-orbit model with the lens parameters shown in the "Wide model best-fit" column of Table 1. As a comparison, the best-fit single source single lens model is also shown in Figure 1. As expected, this single source single lens fit does not explain the anomalous data well.
In Table 1, we also show the best-fit single source binary lens close-orbit model for comparison. The light curve of the best-fit model for the close-orbit is shown in Figure 2. Here we find the wide-orbit model is favored over the close-solution model with a difference in χ 2 of -478.8. Figure 3 shows the source trajectory and magnification pattern of the wide-orbit best-fit solution. Figure 4 shows the equivalent for the close-orbit best-fit solution.
The distribution of MCMC states for the run with covariance can is shown in Figure 5 for the wide-orbit solution, and Figure 6 for the close-orbit solution.
Next, we obtain the dereddened color and corrected magnitude of the source star from the instrumental MOA-II magnitudes, which are calibrated via the procedure described by Bond et al. (2017). The colormagnitude diagram of the stars near MOA-2020-BLG-208 is shown in Figure 7. Using the unextincted red clump centroid as predicted by Nataf et al. (2013), we estimate the color and magnitudes shown in Table 2. Using the Boyajian et al. (2014) restricted analysis of stars with 3900 < T eff < 7000, we use the color and magnitude to estimate the angular size of the source star. Using this, we estimate the Einstein radius of the system. These and other estimated source and source-lens properties are shown in Table 2.

GALACTIC MODEL ANALYSIS RESULTS
Various properties of the lens system, such as lens mass and lens distance, cannot be directly inferred with-  Table 2. Source and lens-source system properties. out microlensing parallax and lens brightness measurements. These measurements do not exist for the MOA-2020-BLG-208 event, as parallax is not reliably detected and the angular source size cannot be directly measured. Indeed, the wide-solution microlensing parallax is not well constrained by our light curve fitting procedure, as is seen by the high uncertainty in the parallax parameters r πE and θ πE in Table 1. Hence, we estimate additional lens system properties using the Bayesian analysis galactic model provided by Bennett et al. (2014). This model allows its posterior distributions to be influenced by a prior based on the host mass. Specifically, it allows the posterior distributions to rely on either a mass function that assumes that all stars have an equal planet-hosting probability or one that assumes planets are more likely to orbit around more massive stars. We infer the lens properties posterior distributions using the MCMC states described in Section 3 as the input, excluding the parallax parameters. Along with the MCMC states as input, we provide the model with the estimated angular source radius and I-band extinction that we calculate as described in Section 2. To facilitate high-angular resolution follow-up observations, we additionally run the galactic model inference for the K-band, using an extinction that we calculate via the method provided by Gonzalez et al. (2012); Nishiyama   Table 3. This table shows the median values of our galactic model distribution as well as the upper and lower bounds of a 68.27% confidence interval (1 σ). These results are for the wide-orbit model.  Table 1).
et al. (2009). For MOA-2020-BLG-208, we calculate this extinction to be A K = 0.2062. For both bands, we infer the lens parameters using two different input mass priors. The difference between these two priors comes in the form of varying α in the power-law stellar mass function defined by P ∝ M α . In one case, we set α = 0, which provides a mass prior that assumes all stars have an equal probability of hosting planets. This is a common prior assumption in many existing planetary microlensing analyses and related statistical population studies (e.g. Cassan et al. 2012). Figure 8 shows the Galactic model posterior distributions with α = 0. The median values of the distributions as well as the upper and lower bounds of a 68.27% confidence interval are shown in Table 3.
However, several works have suggested that the probability of hosting a planet scales in proportion to the host star mass, including cases in radial velocity samples (Johnson et al. 2007(Johnson et al. , 2010, in direct imaging samples (Nielsen et al. 2019), and in revised analyses of microlensing samples with additional follow-up observations , suggesting mass prior with α = 0 may be incorrect. Thus, following the analysis of Silva et al. (2022), we also present an analysis using the power-law stellar mass function where α = 1, which provides a mass prior that assumes more massive stars have a greater probability of hosting planets. Fig-ure 9 shows the Galactic model posterior distributions with α = 1.
In this Bayesian analysis, we assume that the planet hosting probability does not depend on Galactocentric distance ).
The galactic model by Bennett et al. (2014) does not consider the possibility of a nearby source. Therefore, we performed a brief analysis using the galactic model by ; ; Koshimoto & Ranc (2022), which does include the possibility of nearby sources. From this, we found that probability of a source distance D S < 4 kpc is less than 1.31 × 10 −6 , suggesting a nearby source is unlikely.
6. ALTERNATIVE MODELS 6.1. Evidence against a binary source model We also attempted to fit the observed data to a binary source model. The results of this modeling suggest the observed data does not fit well to a binary source model.
For the binary source modeling, we restricted the data sources to observatories that obtained data near the anomalous event, namely, MOA and KMT data. Each of these instruments have relatively narrow pass bands. With these data, we fit both a binary lens and binary source model. The resulting binary lens model fit was similar to the wide-solution model shown in Section 4. A comparison of these two models is shown in Table 4.
The χ 2 value of the binary source model is 154.1 more than that of the binary lens model, suggesting a strong preference toward the binary lens model. Furthermore, this best fit binary source model attributes an improbable blue color to the second source (see Figure 10). Attempts to restrict the second source to a more likely color result in significant increases to the χ 2 value of the model fit.
With both the significant χ 2 preference toward the binary lens model and the best fit binary source model having unlikely physical parameters, we infer this event is unlikely to have been caused by a binary source.

Close-orbit model orbital motion
The parallax parameters of our close-orbit model are relatively constrained. This is likely due to the closesolution model being the incorrect model resulting in spurious fit values, as the wide-solution model presents a significantly improved χ 2 value. However, the constraints on the parallax parameters may also have resulted from a degeneracy with orbital motion, which was not included in the above modeling. To test this possibility, we performed preliminary modeling of the closesolution model with orbital motion being fitted in addition to the other parameters above. During this prelimi- Binary source single lens model best-fit Single source binary lens model best-fit  Table 4. This table shows a comparison of binary-source-single-lens and single-source-binary-lens models.  nary modeling, the model fit tended toward orbital periods which are unlikely according to our galactic model. Based on our orbital period distribution from our galactic modeling results, this lens parameter modeling preferred orbital periods longer than the upper bound of our 95% confidence interval at 2089 days. For practical reasons, our analysis placed an artificial limit at the 99% confidence interval, which prevented the modeling from moving toward even longer orbital periods. Notably, these preliminary fits did not improve the χ 2 compared to the close-solution without orbital motion. This can be seen as further evidence that the close-solution model is unlikely to be the correct model and that the parallax constraints are likely spurious results due to it being the incorrect solution.
6.3. Degenerate close-orbit model trajectory While our fitting procedure is able to obtain multiple χ 2 local minima models, to ensure that it did not exclude any of importance, we explicitly fit common possible degenerate models. Notably, for the close-orbit case, a trajectory that passes on opposite side of the minor caustic compared from the best-fit model is often competitive (e.g., Han et al. 2022;Wang et al. 2022). We performed our fitting process restricting parameter space to keep the resulting model within this local minima, and found that this solution performed significantly worse than our best-fit close model, with a difference in χ 2 of 354.8.

DISCUSSION
Recent exoplanet surveys have provided a plethora of planetary systems, many of which challenge established formation models (Gaudi et al. 2021). Among these surveys, gravitational microlensing surveys are apt to probe planets on significantly wider orbits (Kane 2011). Of particular note is the statistical inferences of the population of cold sub-Saturn-mass planets. Based on a limited number of detections, survey sensitivity analyses (Suzuki et al. 2016;Jung et al. 2019) suggest planets with a host-star mass-ratio, q, larger than 10 −4 are common. Formation models that propose a shortage of cold sub-Saturn-mass planets (Ida & Lin 2004;Mordasini et al. 2009;Ida et al. 2013) would be contested by such a population of planets (Suzuki et al. 2018;Zang et al. 2020;Terry et al. 2021), and revised formation models would be required (Ali-Dib et al. 2022). Notably, Suzuki et al. (2018) presents a mass-ratio gap in models from Ida & Lin (2004); Mordasini et al. (2009) which is not identifiable in microlensing observations from Suzuki et al. (2016). Our analysis suggests a planet mass which contributes to the planet population in this gap with q = 3.17 +0.28 −0.26 × 10 −4 (m planet = 69 +37 −34 M ⊕ when the mass prior uses α = 1 and m planet = 46 +42 −24 M ⊕ when α = 0). Note, the planet mass is sub-Saturn even though the mass-ratio q = 3.17 +0.28 −0.26 ×10 −4 is slightly above that of Saturn. However, this still falls within of the massratio gap described by Suzuki et al. (2018).
As the event was discovered via the MOA alert system, the planet model has a ∆χ 2 < 100 compared to the single lens fit, and data from other observatories support the planet model, the event qualifies for inclusion in the extended MOA-II exoplanet microlensing sample. The extended MOA-II sample is an upcoming statistical analysis of cold exoplanets detected by the MOA-II survey and is the expansion of the Suzuki et al. (2016) sample analysis. Future high-resolution angular followup of the lens and source may help contribute to the population of events that can be used to clarify which of the two mass priors used in this work are more appropriate.

CONCLUSION
In this work, we presented the analysis of the MOA-2020-BLG-208 gravitational microlensing event including the discovery and characterization of a new planet with an estimated mass similar to Saturn. As a cold Saturn-mass planet, this planet contributes to the evidence supporting the need for revised planetary formation models. The planet also qualifies for inclusion in the extended MOA-II exoplanet microlensing sample. We have provided evidence that the anomaly in the event is best described by a planet in a wide-solution orbit as opposed to other potential models, such as a binary source model. Our characterization was derived both from light curve modeling analysis and galactic model analysis. Notably, our galactic modeling included results from two potential mass law priors.

ACKNOWLEDGMENTS
We note that observations used in this analysis were obtained during early phases of the COVID-19 pandemic, which introduced logistical challenges for many observatories. Several collaborating observatories were unable to take measurements during this period, and those that did endured additional obstacles to do so. The MOA project is supported by JSPS KAK-ENHI Grant Number JSPS24253004, JSPS26247023, JSPS23340064, JSPS15H00781, JP16H06287, and JP17H02871. Work by C.R. was supported by a Research fellowship of the Alexander von Humboldt Foundation.