MOA-2019-BLG-008Lb: a new microlensing detection of an object at the planet/brown dwarf boundary

We report on the observations, analysis and interpretation of the microlensing event MOA-2019- BLG-008. The observed anomaly in the photometric light curve is best described through a binary lens model. In this model, the source did not cross caustics and no finite source effects were observed. Therefore the angular Einstein ring radius cannot be measured from the light curve alone. However, the large event duration, t E about 80 days, allows a precise measurement of the microlensing parallax. In addition to the constraints on the angular radius and the apparent brightness I s of the source, we employ the Besancon and GalMod galactic models to estimate the physical properties of the lens. We find excellent agreement between the predictions of the two Galactic models: the companion is likely a resident of the brown dwarf desert with a mass Mp about 30 MJup and the host is a main sequence dwarf star. The lens lies along the line of sight to the Galactic Bulge, at a distance of less then4 kpc. We estimate that in about 10 years, the lens and source will be separated by 55 mas, and it will be possible to confirm the exact nature of the lensing system by using high-resolution imaging from ground or space-based observatories.


INTRODUCTION
During the last 20 years, thousands of planets 1 have been detected and it is now clear that planets are abundant in the Milky Way (Cassan et al. 2012;Bonfils et al. 2013;Clanton & Gaudi 2016;Fulton et al. 2021). Conversely, the various methods of detection agree that brown dwarf companions (with a mass ∼13-80 M Jup ) seem much rarer (Grether & Lineweaver 2006;Lafrenière et al. 2007;Kraus et al. 2008;Metchev & Hillenbrand 2009;Kiefer et al. 2019;Nielsen et al. 2019;Carmichael et al. 2020), inspiring the idea of a "brown dwarf desert" (Marcy & Butler 2000), and such disparity raises questions about formation scenarios. Core accretion, disc instability, migration and disc evolution mechanisms are capable of producing planets up to 40M Jup (Pollack et al. 1996;Boss 1997;Alibert et al. 2005;Mordasini et al. 2009), explaining the formation of some brown dwarf companions. Brown dwarfs can also form via gas-collapse (Béjar et al. 2001;Bate et al. 2002) and several processes have been proposed to explain the cessation of gas accretion, such as ejection (see Luhman (2012) and references therein for a more complete review). However, the formation of low-mass binaries remains difficult to explain (Bate et al. 2002;Marks et al. 2017) and more detections are needed to place meaningful constraints on formation models, especially around the brown dwarf desert.
Several objects at the planet/brown dwarf mass boundary have been discovered with the microlensing technique, both in binary and single lens events (Bachelet et al. 2012a;Bozza et al. 2012;Ranc et al. 2015;Han et al. 2016;Zhu et al. 2016;Poleski et al. 2017;Shvartzvald et al. 2019;Bachelet et al. 2019;Miyazaki et al. 2020). Microlensing is particularly sensitive to exoplanets and brown dwarfs at or beyond the snow-line of their host stars, which is the region beyond which it is cold enough for water to turn to ice. Planets in this region typically have orbital periods of many years and, as such, are mostly inaccessible to other planet detection methods (Gould et al. 2010;Tsapras et al. 2016).
The location of the snow-line plays an important role during planet formation, as the prevalence of ice grains beyond that point is believed to facilitate the formation of sufficiently large planetary cores, able to trigger runaway growth and form giant planets (Ida & Lin 2004;Kley & Nelson 2012).
The lensing geometry is typically expressed in terms of the angular Einstein radius of the lens (Einstein 1936) where D L , D S are the distances from the observer to the lens and source respectively, D LS is the lens-source distance, and M L is the mass of the lens. The key observable in microlensing events that provides any connection to the physical properties of the lens is the event timescale where µ rel is the relative proper motion between lens and source, π rel is the lens-source relative parallax and κ ≡ 4G/c 2 au 8.1mas/M (Gould 2000). These two equations reveal a well-known degeneracy in microlensing event parameters. Indeed, the mass of the lens, its distance and the distance to the source are degenerate parameters when only t E is measured and at least two extra pieces of information are required to disentangle them. For binary lenses, θ E = θ * /ρ is often measured from the detection of finite source effects in the event light curve, typically parametrized with ρ. This occurs when an extended source of angular radius θ * closely approaches regions of strong magnification gradients, i.e. around caustics (Witt 1990;Tsapras 2018). Using a color-radius relation (Boyajian et al. 2014), it is then possible to estimate θ E . For sufficiently long events (i.e., t E ≥ 30 days), the microlensing parallax π E = π rel /θ E due to Earth's revolution around the Sun, can be measured. This is referred to as the annual parallax (Gould 1992(Gould , 2004). In addition, if simultaneous observations can be performed from space, as well as from the ground, it is possible to measure the space-based parallax (Refsdal 1966;Calchi Novati et al. 2015;Yee et al. 2015a). Ultimately, by obtaining high-resolution imaging several years after the event has expired, additional constraints on the relative lens-source proper motion µ rel and the lens brightness may be obtained, provided that the lens and source can be resolved (Alcock et al. 2001). See for example Beaulieu (2018) for a complete review of this technique.
It is not rare however that only t E and a single other parameter (θ E or π E ) are measured, leaving the physical parameters of the lens system only loosely constrained. As underlined by Penny et al. (2016), this is the case for about 50 % of all published microlensing planetary events. To obtain stronger constraints on these events, Galactic models may be employed to derive the probability densities of lens mass and distance along the line of sight that reproduce the fitted microlensing event parameters. Originally used by Han & Gould (1995), Galactic models are now commonly relied upon to estimate the properties of microlensing planets when no additional information is available to constrain the parameter space (Penny et al. 2016). While they generally come with large uncertainties (10% is a lower limit), Galactic model predictions have proven to be in excellent agreement with results obtained from follow-up studies using high-resolution imaging, especially for OGLE-2005-BLG-169Lb (Gould et al. 2006;Batista et al. 2015;Bennett et al. 2015), MOA-2011-BLG-293Lb (Yee et al. 2012Batista et al. 2014) and OGLE-2014-BLG-0124Lb Beaulieu et al. 2018). Galactic models developed for microlensing analysis are employed to generate distributions of stellar densities and velocities across the Galactic Disk and Bulge (Han & Gould 1995, 2003Dominik 2006;Bennett et al. 2014), and use them to reproduce microlensing observables (i.e. t E , θ E and π E ). These are then compared to the fitted event parameters in order to estimate the probability densities of the lens distance D l and mass M l . In addition to these models, there exist synthetic stellar population models for the Milky Way that have been explicitly developed to reproduce observable Galactic properties with great accuracy. Specifically, the Besançon (Robin et al. 2003) and GalMod (Pasetto et al. 2018) models have been used in many different studies, to explore the structure, kinematics and formation history of the the Milky Way (Czekaj et al. 2014;Robin et al. 2017). In addition, they have also been used to simulate astronomical sky-surveys (Rauer et al. 2014;Penny et al. 2013Penny et al. , 2019Kauffmann et al. 2020), and their predictions have been tested against real observations Bochanski et al. 2007;Pietrukowicz et al. 2012;Schmidt et al. 2020;Terry et al. 2020).
For the first time, in this study we employ both the Besançon and GalMod Galactic models to estimate the properties of a binary lens, with a companion likely located in the brown dwarf desert. The microlensing event MOA-2019-BLG-008 was observed by several microlensing teams independently, and we present the different data sets, as well as the data reduction procedures, in Section 2. The modeling of the photometric light curve and the model selection are discussed in Section 3. Section 4 presents the analysis of the properties of the source and of the blend contaminant. The methodology used to derive the physical properties of the lens system is detailed in Section 5, and we conclude in Section 6.

Survey and follow-up observations
The microlensing event MOA-2019-BLG-008 was first announced on 4 Feb 2019 by the MOA collaboration (Sumi et al. 2003), which operates the 1.8-m MOA survey telescope at Mount John observatory in New Zealand, at equatorial coordinates α = 17 h 51 m 55.89 s , δ = −29 • 59 23.03 (J2000) (l, b = 359.8049 • , −1.7203 • ). The event was also independently identified by the Early Warning System (EWS) 2 of the Optical Gravitational Lensing Experiment (OGLE) survey (Udalski 2003;Udalski et al. 2015) as OGLE-2019-BLG-0011. OGLE observations were carried out with the 1.3-m Warsaw telescope at Las Campanas Observatory in Chile, with the 32-chip mosaic CCD camera. The event occurred in OGLE bulge field BLG501, which was imaged about once per hour when not interrupted by weather or the full Moon, providing good coverage of the light curve when the bulge was visible from Chile.
Additional observations were obtained by the ROME/REA survey ) using 6×1m telescopes from the southern ring of the global robotic telescope network of the Las Cumbres Observatory (LCO) (Brown et al. 2013). The LCO telescopes are located at the Cerro Tololo International Observatory (CTIO) in Chile, South African Astronomical Observatory (SAAO) in South Africa and Siding Spring Observatory (SSO) in Australia, and they provided good coverage of the light curve, although the event occurred early in the 2019 ROME/REA microlensing season (i.e. ∼ March to September of each year, when the Galactic Bulge is observable). Obsevrations were acquired in the survey mode.
The event lies in fields BLG02 and BLG42 of the Korea Microlensing Telescopes Network (KMTNet) (Kim et al. 2016) and so, was intensely monitored by that survey, although KMTNet did not independently discover the event.
Observations were also obtained from the Spitzer satellite as part of an effort to constrain the parallax (Yee et al. 2015b). Spitzer observations will be presented in a companion paper (Han et al. 2022 in prep.).

Data reduction procedure
This analysis uses all available ground-based observations of MOA-2019-BLG-008. The list of contributing telescopes is given in Table 1. Most data were obtained in the I band (or SDSS-i) but we note that MOA observations were performed with the MOA wide-band red filter, which is specific to that survey (Sako et al. 2008). ROME/REA obtained observations in three different bands (SDSS-i , SDSS-r and SDSS-g ). The KMTNet survey observations were carried out in the I band, with a complementary V band observation every ten I exposures.
The photometric analysis of crowded-field observations is a challenging task. Images of the Galactic bulge contain hundred of thousands of stars whose point-spread functions (PSFs) often overlap, therefore aperture and PSF-fitting photometry offer very limited sensitivity to photometric deviations generated by the presence of low-mass planetary companions. For this reason, observers of microlensing events routinely perform difference image analysis (DIA) (Tomaney & Crotts 1996;Alard & Lupton 1998;Bramich 2008a;Bramich et al. 2013a), which offers superior photometric precision under such crowded conditions. Most microlensing teams have developed custom DIA pipelines to reduce their observations. OGLE, MOA and KMT images were reduced using the photometric pipelines described in Udalski (2003), Bond et al. (2001) and Albrow et al. (2009), respectively. The LCO observations were processed using the pyDANDIA pipeline (ROME/REA in prep), a customized re-implementation of the DanDIA pipeline (Bramich 2008b;Bramich et al. 2013b) in Python. The data sets presented in this paper have been carefully reprocessed to achieve greater photometric accuracy, and it is these data that we used as input when modelling the microlensing event. They are available for download from the online version of the paper.
We note the presence of a very long-term baseline trend spanning several observing seasons in the OGLE and MOA photometry that can be seen in Figure 1. As described later in this work, we determined that the source star is a red giant. Many red giants exhibit variability at the ∼ 10% level (Wray et al. 2004;Percy et al. 2008;Wyrzykowski et al. 2006;Soszyński et al. 2013;Arnold et al. 2020), and it is possible that this is also the case for this source, despite the apparently very long period P ≥ 1000 days. Because this trend manifests over very long time scales (several years), much longer than the duration of the actual microlensing event (weeks), it does not have any noticeable effect on the determination of the parameters of this event, which are primarily derived from the detailed morphology of the microlensing light curve. The baseline over the duration of the microlensing event is effectively flat. Therefore, to increase the speed of the modeling process, we only used observations with JD ≥ 2457800 and included data sets with more than 10 measurements in total during the course of the microlensing event. The latter constraint applies only to the LCO data and is limited to the observations acquired by the reactive REA mode on a different target in the same field . We verified that our data selection does not impact the overall results, by exploring models with the full-baseline.

MODELING THE EVENT LIGHT CURVE
This event displays a clear anomaly around HJD ∼ 2458580, implying that it is most likely due to a binary lens (2L1S) or a binary source (1L2S) (Dominik et al. 2019). It is morpholigicaly similar to the event MACHO 99-BLG-47 (Albrow et al. 2002), despite a different lensing geometry, as detailed below. In addition, because the event lasts for ∼ 300 days, the effect of the motion of Earth around the Sun, referred to as the annual parallax (Gould 1992;Alcock et al. 1995), needs to be taken into account. The classical approach to modeling is to first search for static binary models and subsequently gradually introduce additional second-order effects, such as parallax or the orbital motion of the lens (Dominik 1999). To model the event we use the pyLIMA software (Bachelet et al. 2017), which employs the VBBinaryLensing code (Bozza 2010;Bozza et al. 2018) to estimate the binary lens model magnification, and we search for a general solution including the annual parallax, but we also explore the static case for completeness. The first step of modeling involves identifying potential multiple minima in the parameter space, and for this we employ the differential evolution algorithm (Storn & Price 1997;Bachelet et al. 2017). During the modeling process, we rescale the data uncertainties using the same method presented in Bachelet et al. (2019), which introduces the parameters k and e min for each datasets: where σ and σ are the original and rescaled uncertainties (in magnitude units), respectively. The coefficients for each dataset are given in Table 1. Finally, we explore the posterior distribution of the parameters of each minimum that we identify using the emcee software (Foreman-Mackey et al. 2013).

(No) finite-source effects
In principle, the normalized angular source radius ρ = θ * /θ E has to be considered (Witt & Mao 1994), but preliminary models indicated that this parameter is loosely constrained. This is because the source trajectory does not pass close to caustics, as can be inferred from Figure 2. However, there exists an upper limit ρ m where the finite-source effects would start to be significantly visible in the models. Because θ E ≥ θ * /ρ m , this limit introduces constraints on the mass and distance of the lens that can be used for the analysis presented in Section 5. Indeed, it is straightforward to Table 1. Summary of telescopes and observations used for modeling the event. The number of data points per telescope represents the points used for the modeling step., i.e. JD ≥ 2457800. Lines marked with '*' indicate that this data set was not used during the modeling process, as described in the text. In cases which the rescaling parameters were not constrained, they were fixed to k = 1.0 and emin = 0.0.   (Gould 2000): where π l and π s are the parallax of the lens and source, respectively. The constraint on the mass can be written as Therefore, we sample the distribution of ρ around the best model and found that ρ m ≤ 0.01 (with a conservative 10σ limit) and consider the source as a point for the rest of the modeling presented in this analysis.

Binary lens
A binary lens model involves seven parameters. t 0 is the time when the angular distance u 0 (scaled to θ E between the source and the center of mass of the binary lens is minimal. The event duration is characterised by the angular Einstein ring radius crossing time t E = θ E /µ rel , where µ rel is the lens/source relative proper motion (in the geocentric frame, because pyLIMA follows the geocentric formalism of Gould 2004. The binary separation projected on the plane of the lens is defined as s and the mass ratio between two component as q. The angle between the trajectory and the binary axis (fixed along the x-axis) is defined as α. We also consider a source flux f s and a blend flux f b for each of the datasets, adding 2n parameters where n is the number of datasets (i.e., 29 in this study). As discussed previously, we neglect the last parameter ρ and fit the simple point-source binary lens model.
Following Gould (2004), we define the parallax vector π E by its North (π E,N ) and East (π E,E ) components. We set the parallax reference time as t 0,par = 2458570 HJD (Skowron et al. 2011) for all models considered in this analysis. This coincides with the time of the anomaly peak and corresponds to the calendar date 28 March 2019. At this time, Earth's acceleration vector was nearly parallel to the East direction. Therefore, we expect the π EE component to be the better constrained of the two. We found that the lightcurve morphology can only be explained by a single geometry, in agreement with previous results from real time modeling conducted by V. Bozza 3 and C. Han 4 . However, a second solution exists, a consequence of the binary ecliptic degeneracy (Skowron et al. 2011) with (u 0 , α, π E,N ) ⇐⇒ −(u 0 , α, π E,N ). Because the magnification pattern is symmetric relative to the binary axis, it exists two source trajectories that produce identical lightcurves for static binaries, i.e. (u 0 , α) ⇐⇒ −(u 0 , α). Moreover, the projected Earth acceleration can be considered as almost constant during the duration of the event, leading to the π E,N degeneracy for events located towards the Galactic Bulge (Smith et al. 2003;Jiang et al. 2004;Gould 2004;Poindexter et al. 2005;Skowron et al. 2011). This degeneracy is especially severe for events occurring near the equinoxes, because the projected Earth acceleration varies slowly (Skowron et al. 2011). We therefore explore both solution in the following analysis. Note that pyLIMA uses the formalism of Gould (2004) and therefore u 0 ≥ 0 if the lens passes the source on its right. Given the relatively long timescale of the event, we also explored the possibility of orbital motion of the lens (Albrow et al. 2000;Bachelet et al. 2012b) and considered the simplest linear model, parametrized with ds/dt and dα/dt, the linear separation and angular variation over time at time t 0 . For completeness, we explored the parameter space for a static model (i.e., without the annual parallax) and found that the best model converges to a similar geometry. However, the residuals display systematic trends around the event peak that are the clear signature of annual parallax, which is reflected in the high χ 2 value presented in Table 2.

Binary source
We explored the possibility that the observed light curve was due to a binary source (Gaudi 1998). Following the approach described in Hwang et al. (2013), we added to the single lens model the extra parameters ∆t 0 and ∆u 0 , respectively the shifts in the time of peak and separation of the second source relative to the first. Finally, we also added the flux ratio of the two sources in each observed band q λ . We report our results in Table 2. We also explore models with two different source angular sizes, but did not find any significant improvements in the model likelihood.  Table 2. Models parameters are defined as 16,50 and 84 MCMC chains percentile, except for the χ 2 which is reported as the minimum value (i.e., the best model in each case). The angular source radius θ * for each model is also presented, see Section 4. The two static models are denoted 2L1S, the two models with parallax are denoted 2L1SP and models with orbital motion of the lens are 2L1SPOM. + and -indicated positive or negative u0.

ANALYSIS OF THE SOURCE AND THE BLEND
In the analysis of microlensing events, the color-magnitude diagram (CMD) is used to estimate the angular source radius θ * , and ultimately the angular Einstein ring radius θ E = θ * /ρ (Yoo et al. 2004). Unfortunately, ρ is not measurable in the present case. Nonetheless, the CMD analysis provides useful information about the source and the lens that can be used to place additional constraints to the analysis presented in Section 5. The CMD constraints can also be used to inform observing decisions in the future with complementary high-resolution imaging. We conducted our CMD analysis using different and independently obtained sets of observations from our pool of available data sets. The estimation of θ * below is for the model 2L1S with parallax and u 0 < 0. The source and blend magnitudes for all models are presented in Table 3 and the angular source radius θ * derived for all models is presented in Table 2.

ROME/REA Color-Magnitude Diagrams Analysis
The ROME strategy consists of regular monitoring of 20 fields in the Galactic bulge in SDSS-g , SDSS-r and SDSS-i , as described in Tsapras et al. (2019), and is designed to improve our understanding of the source and blend properties. The photometry is obtained using the pyDANDIA algorithm (Street et al. in prep 5 ) and calibrated to the VPHAS+ catalog (Drew et al. 2016). For this event, we investigated all combinations of filters and telescope sites and selected LSC A (i.e., LCO dome A in Chile) for the ROME CMD analysis, as it provided the deepest catalog. Figure 4 presents the CMD for stars in a 2'x2' square centered on the target, while Figure 3 presents a composite image of the ROME observations from LSC A. The latter presents the variable extinction in the field of view, as well as the two clusters NGC 6451 and Basel 5.
The first step is to estimate the centroid of the Red Giant Clump (RGC). In Figure 4, stars located within 2' from the event location from the ROME/REA and VPHAS catalogs are displayed in (r-i,i) and (g-i,i) CMDs. We note that the location of the RGC is quite uncertain in the g-band for the ROME/REA data. This is due to the high extinction along this line of sight and leads to an accurate g-band calibration of the ROME data. We use the VPHAS magnitudes of these stars to estimate the centroid positions of the RGC in the three bands. We found the magnitudes of the RGC to be ( mag, E(r − i) = 1.41 ± 0.1 mag and extinction A i = 2.3 ± 0.1 mag. The best model returns a source magnitude of (g, r, i) s = (24.2 ± 0.8, 19.47 ± 0.05, 17.31 ± 0.03) mag. Because the event occured at the begining of the season, the event was poorly covered by the ROME data in the g-band and the source brightness is not well constrained. However, we can use the (r-i) color and i magnitude to estimate θ * . As described in Appendix A), we used the same catalog as Boyajian et al. (2014) to construct a new color-radius relation: This relation returns θ * = 9.8 ± 1.2 µas while the second relation using the g-band returns θ * = 22.9 ± 8.4 µas. The latter value is inaccurate due to large uncertainties in the source brightness in the g-band.

MOA color-Magnitude Diagram
The MOA magnitude system can be transformed to the OGLE-III magnitude system (i.e., the Johnson-Cousins system) by using the relation presented in Appendix B. Using the intrinsic color and magnitude of the RGC ((V − I) 0 , I 0 ) = (1.06, 14.32) mag (Nataf et al. 2013;Bensby et al. 2013) and subtracted to the measured position the RGC centroid in the Figure 5, we could estimate (E(V − I), A I ) = (2.3 ± 0.1, 2.4 ± 0.1) mag, in good agreement with the previous estimation. Knowing that the transformed magnitudes of the source are (V, I) s = (20.8 ± 0.1, 16.94 ± 0.08) mag, we found (V, I) 0,s = (16.1 ± 0.1, 14.54 ± 0.08) and we ultimately estimate θ * = 8.9 ± 1.3 µas (Adams et al. 2018), in relative agreement with the previous estimation. In Figure 5, we also display the source and blend position using the measurements from OGLE-IV in the I band and the transformed MOA V band. We estimate the source to be (V, I) s = ((20.8 ± 0.1, 16.88 ± 0.01) mag and derive θ * = 10.4 ± 1.6 µas. This estimate is likely more accurate than the previous, because it relies on a single color transformation (with the highest color term in Equation B4).

Gaia EDR3
The Gaia mission (Gaia Collaboration et al. 2016) recently released their "Early Data Release 3" data set (EDR3,(Gaia Collaboration et al. 2020)), which significantly increases the volume and precision of the Gaia catalog. We queried the Gaia catalog requesting all stars within a 3' radius 6 around the coordinates of the event to generate a Gaia CMD, which we present in Figure 6. We limit our study to stars with a Re-normalised Unit Weight Figure 3. Color composite of g,r,i reference images of the ROME survey. The inset is 2'x2' zoom around MOA-2019-BLG-008. The NGC6451 cluster center is visible, while some of the stars of the Basel5 cluster are also visible (its centre is on the left, outside of the image, see Kharchenko et al. (2013)). As indicated by the white cross, North is up and East is left.
Error (RUWE, a statistical criterion of the data quality) better than 1.4 7 . MOA-2019-BLG-008 is in the catalog (at 63 mas, Gaia EDR3 4056394717636682112) and appears in a sparse location of the CMD. The reported parallax is p = 0.39±0.12 mas, which corresponds to a distance of D = 2.56±0.79 kpc. This object is also significantly redder and brighter than the blend discussed in the previous section. Indeed, by using the magnitude transformation from Gaia to the Johnson-Cousins system 8 (Bachelet et al. 2019), we found this object to be ((V − I), I) = (2.25 ± 0.07, 16.17 ± 0.05) mag, and this is very likely the sum of the source and the blend previously discussed ((V − I), I) tot = (2.30, 16.19) mag. This is confirmed by several useful metrics available in the catalog. First, we compute the corrected BP and RP excess flux factor (Evans et al. 2018;Riello et al. 2020) and find 0.28 9 , which corresponds to a blend probability of ∼ 0.3 (see Figure 19 of Riello et al. (2020). Secondly, we note that the fraction R of visits that indicated a significant blend (defined by 'phot bp n blended transits' and 'phot rp n blended transits' for the BP and RP bands respectively) divided by the number of visits used for the astrometric solution ('astrometric matched transits') is very high for both bands R = 36/39 ∼ 90%.
Following Mróz et al. (2020), we plot in Figure 6 the distribution of galactic proper motions of the Disk (gray) and Bulge (red) populations. Note that Gaia EDR3 provides proper motion in equatorial coordinates which we transformed using the same method as in Bachelet et al. (2019). The Disk population, approximated by the main sequence population, is estimated from the Gaia CMD by using all stars with G BP −G RP ≤ 2.5. The Bulge population is estimated from the RGC population of the CMD (i.e., 2.8 < G BP − G RP ≤ 3.5 and 17.8 ≤ G ≤ 19). Because this object is blended, it is difficult to extract meaningful constraints from the proper motion distribution. However, it will be possible to do so when the source and lens will be sufficiently separated, as discussed later.

Analysis of the blend
As can be seen in the different CMD's, there is a significant blend flux in the datasets, which likely belongs to the population of foreground stars of the Galactic Disk. The measurements from MOA/OGLE are (V, I) b = (18.5 ± Figure 4. Color-Magnitude diagrams from the ROME/REA survey. The blue dots represents all stars within 2' around the event location from the VPHAS catalog (Drew et al. 2016), while the aligned ROME/REA stars are in orange. The source and blend are represented in magenta and blue, respectively. The dashed square on the left plot represents the stars that have been used to estimate the RGC centroid for both CMDs. 0.1, 16.88 ± 0.01) mag and are consistent with a late F dwarf located at ∼ 2.5 kpc (Bessell & Brett 1988;Pecaut & Mamajek 2013) (assuming half extinction). Measurements from the ROME survey indicate that the blend brightnesses are (g, r, i) b = (19.43 ± 0.01, 17.91 ± 0.01, 16.97 ± 0.02) mag, consistent with a G dwarf at ∼ 2.0 kpc (Finlator et al. 2000;Schlafly et al. 2010). These results are in agreement with the lens properties derived in the next section.
The source and blend have similar brightness in Cousins I band, but the former is much redder. Therefore, the object identified at this location by Gaia is dominated by light from the blend. At the epoch J2016, Gaia measures a total offset of 60 ± 15 mas with respect to the event location measured in 2019 (i.e., during peak magnification). The reported error on the distance has been computed from the North and East components error from ground surveys (of the order of ∼ 15 mas) and neglecting Gaia errors (of the order of ∼ 0.1 mas). Similarly, the magnified source and the baseline object in the KMTNet images are separated by ∼0.068 pixels, which is equivalent to 27 mas. Because the blend and source have similar brightness, this indicates a separation between the blend and the source of ∼ 60 mas. Therefore, the hypothesis that the blend, or a potential companion, is the lens is probable. Assuming that the blend is the lens, π E ∼ 0.2, a source distance D s ∼ 8 kpc and that the light detected by Gaia is solely due to the blend, we can estimate θ E ∼ (0.39 − 0.125)/0.2 ∼ 1.3 mas and M l ∼ 0.8M . So the astrometric solution is also compatible with a G dwarf at ∼ 2.5 kpc, with the notable exception of the relative proper motion. Indeed, we can estimate the geocentric relative proper motion to be µ geo = θ E /t E ∼ 1.3/80 ∼ 6 mas/yr. Because this event peaked in early March, the heliocentric correction v ⊕ π rel /au ∼ (0.2, −0.4) is small and we can therefore assume µ geo = µ hel (Dong et al. 2009). The separation between the lens and the source at the Gaia epoch (J2016) is therefore expected to be ∼ 18 mas, significantly smaller than the previously measured ∼60 mas. However, this argument can not, by itself, rule out the hypothesis that the blend is the lens due to the relatively large errors.
Therefore, both photometric and astrometric arguments provide sufficient evidence that the lens represents a significant fraction of the blended light, but only high-resolution imaging in the near future will provide a conclusive answer to this puzzle. Figure 5. Color-Magnitude Diagram from the MOA survey, calibrated to the OGLE III system, for stars located in the square 2'x2' around the event. We also display the position of the source and the blend using the I measurements from the OGLE lightcurve. Table 3. Source and blend magnitudes for the three 2L1S models (results are almost identical for u0 < 0 and u0 > 0). Numbers in brackets represent the 1σ errors. MOA magnitudes have been converted to the OGLE-III system using the transformation in the Appendix B. Because the normalised radius of the source ρ can not be estimated from the fit, inferring the Einstein radius is not possible without extra measurements, such as the microlensing astrometric signal (Dominik & Sahu 2000), the lens flux or the lens and source separation measurements after several years (Alcock et al. 2001;Beaulieu 2019). In order to estimate the physical properties of the lens, prior information from Galactic models can be used. By drawing random source-lens pairs from distributions of stellar physical parameters derived from the galactic models along the line of sight, and calculating the respective microlensing model parameters, the lens mass and distance probability densities can be estimated. This has been done many times in the past with parameterized models specifically designed to study microlensing events (Han & Gould 1995, 2003Dominik 2006;Bennett et al. 2014;Koshimoto et al. 2021). But there are also modern galactic models that have been extensively tested and are publicly accessible. From a theoretical point of view, these elaborate models are of great interest because they are including more relevant quantities such as color, extinction and stellar type, for instance. These quantities can be used to constrain physical parameters, but also to predict properties for follow-up observations in the more distant future. In this work, we performed a parallel analysis using the parametric model of Dominik (2006), the Besançon model and the GalMod model described thereafter.

The Besançon Model
The first galactic model we use to generate a stellar population is the Besançon Model 10 (Robin et al. 2003), version M1612. This version consists of an ellipsoidal Bulge titled by ∼ 10 • from the Sun-Galactic center direction, and populated with stellar masses drawn from a broken power law initial mass-function (IMF) dN/dM ∝ M α , with α = −1 and α = −2.35 for 0.15 ≤ 0.7M and M > 0.7M respectively (Robin et al. 2012;Penny et al. 2019). The Disk is modelled by a thin disk component with a two-slope power law IMF, with α = 1.6 and α = 3.0 for M ≤ 1M and M > 1M (Robin et al. 2012), while the density distribution is derived from Einasto (1979). The outer part of the disk model has recently been updated and is described in Amôres et al. (2017). The thick disk and halo population are fully described in Robin et al. (2014), while the kinematics of the population are described in (Bienaymé et al. 2015). We select the Marshall et al. (2006) 3D map to estimate the extinction for the simulation. Finally, we note that the Besançon Model has been used in several studies for microlensing predictions. Based on the original work of Kerins et al. (2009), Awiphan et al. (2016 and Specht et al. (2020) developed the MaBµls − 2 software that computes theoretical maps of the distribution of optical depth, event rate and timescales of microlensing events, that are in good agreement with observations. In particular, MaBµls − 2 predictions of event rate and optical depth are excellent agreement with the 8 years of observations from the OGLE survey (Mróz et al. 2019). The Besançon Model has also been used by Penny et al. (2013) to simulate the potential yields of a microlensing exoplanet survey with the Euclid space telescope. More recently, Penny et al. (2019) and Johnson et al. (2020) used an updated version of the Besançon Galactic model to estimate the expected number of detections of bound and unbound planets from the Roman (formerly known as WFIRST) microlensing survey (Spergel et al. 2015).

GalMod
The second simulation was made using the "Galaxy Model" (GalMod, version 18.21), which is a theoretical stellar population synthesis model (Pasetto et al. 2018) 11 simulating a mock catalog for a given field of view and photometric system. Similarly as for the Besançon model, the parameter range in magnitude and color permits the simulation of faint lens stars down to the dwarf and brown dwarf regime. Briefly, GalMod consists of the sum of several stellar populations including a thin and a thick disk, a stellar halo, and a bulge immersed in a halo of dark matter. Stars are generated using the multiple-stellar population consistency theorem described in Pasetto et al. (2019) with a kinematics model from Pasetto et al. (2016). For our simulation, we used the Rosin-Rammler star formation rate (SFR) (Chiosi 1980) for the Bulge and the tilted bar. The thin disk is a combination of five different stellar populations with various ages and kinematics, with a constant SFR, while the thick disk is drawn from a single population. We used the same IMF for all the different components of the model (Kroupa 2001).

Methodology
We first requested samples from the two models within a 2 cone along the line of sight to the event, and set the maximum distance to 10 kpc. We then draw samples of lens and source star combinations and apply a sequence of rules. First, the source has to be more distant than the lens. Then, we consider an event only if the angular separation between the source and the lens is below 10". Following the approach described in Shin et al. (2019), we proceed to compute the associated event parameters (i.e., t E , π E , θ * and I s in this case) and compare them with our measured observables derived from modeling. Each such combination contributes to the final derived distribution with a weight i /2), with δ 2 i being the Mahalanobis distance: are the differences between the best fit model parameters and the simulated parameters and C is the covariance matrix. Note that we also reject models with ρ ≥ 0.01, following the discussion presented in the Section 3.1. Because the galactic models return a finite number of stars (168134 for the Besançon model and 64679 for the GalMod model) and the event parameters are slightly unusual (with t E ∼ 80 days and θ * ∼ 10 µas), a large fraction of lens and source combination have a null weight. For instance, the Besançon model contains only ∼ 0.4 % of stars with |θ * − 11.5| < 3σ and about 1% of events are expected to have t E ∼ 80 days (Mróz et al. 2019). Therefore, the vast majority of trials (≥ 99.9 %) have null weights w i and it would therefore require several thousands of billions of trials to obtain meaningful parameter distributions. In light of this, we adjusted our strategy and adopted an MCMC approach. Using the Mahalanobis as the log-likelihood, we adapt the galactic models to define priors on the modeling parameters: the source and lens distances D s and D l , the proper motion of the source and lens, the mass of the lens M l , the angular radius of the source θ * and the magnitude of the source I s . We use a Kernel Density Estimation (KDE) algorithm to derive continuous distributions from the galactic-model samples. This allows a prior estimation across the entire parameter space, but at the cost of a somewhat smoother distribution and the use of extrapolation.

Results
We present the posterior distribution for the model 2L1S+POM in Figure 7 and the derived results for all models in Table 4. Results from the Galactic model of Dominik (2006) is also presented for comparison. As a supplementary control, we also derived posterior distributions from the recent Galactic model of Koshimoto et al. (2021), especially designed for microlensing studies, and found consistent results.
Despite the relatively broad distributions, we find that galactic models are in good agreement for all models. The main differences are seen in the distribution of the source and lens proper motions. The GalMod model has a much narrower distribution of stellar proper motions, as can be seen in the first line of Figure 7, which directly propagates to the source and lens proper motions. However, the relative proper motion of the two galactic models are in 1-σ agreement, at ∼5.5 mas/yr. This is because the relative proper motion is strongly constrained by t E , which is well determined from the models in this case. For all microlensing models, the derived mass and distance of the host is compatible with the measured blend light, with the exception of the 2L1S-P model, which is much fainter with V l ∼ 20.5. The companion is an object at the planet/brown dwarf boundary.
While the binary source and static binary models can be safely discarded due to their very high ∆χ 2 values relative to the best-fit model, the selection of the best overall model with/without orbital motion (∆χ 2 ∼ 200) and the sign of u 0 (∆χ 2 ∼ 50) is less trivial. In principle, the ∆χ 2 between the various models is statistically significant. However, despite error-bar rescaling, data set residuals can be affected by low level systematics, leading to errors that are not normally distributed (Bachelet et al. 2017). Because they modify the source trajectory in a similar way, the orbital motion and parallax parameters are often correlated (Batista et al. 2011) which is also the case for this event. The North component π E,N of the parallax-vector is in agreement between all non-static models, suggesting that the parallax signal is strong and real, which is expected for an event duration of t E ∼ 80 days. For models including orbital motion, we can use the results from galactic models to verify if the system is bound or not (Dong et al. 2009;Udalski et al. 2018). The condition for bound systems is (K E + P E ) < 0, where K E and P E are the kinetic and potential energies, and this can be rewritten in terms of (projected) escape velocities ratio (Dong et al. 2009 ∼ 6 for the 2L1S-POM and 2L1S+POM models, respectively. Taken at face values, the ratio of projected velocities indicates that the companion is not bound and that these models are unlikely. However, the relative errors are large, i.e. ≥ 100%, and therefore the models with orbital motion can not be completely ruled out. But given the relatively low improvement in the χ 2 of the orbital motion models, we decided to not explore more sophisticated models, such as the full Keplerian parametrization (Skowron et al. 2011).

CONCLUSIONS
We presented the analysis of the microlensing event MOA-2019-BLG-008. The modeling of this event supports a binary-lens interpretation with a mass ratio q 0.04 between the two components of the lens. Because the source a The source distance and errors have been fixed for Dominik (2006) trajectory did not approach the caustics of the system, finite-source effects were not detected, so the lens mass and distance could only be weakly constrained. We used the Besançon and GalMod synthetic stellar population models of the Milky Way to estimate the most likely physical parameters of the lens. By using samples generated by these models, in combination with available constraints on the event timescale t E , the microlensing parallax π E , the source magnitude I s and angular radius θ * , we were able to place constraints on the lens mass and distance. We found that all galactic models, including the one from Dominik (2006), converge to similar solutions for the lens mass and distance, despite different hypotheses (especially for stellar proper motions). We explore several microlensing binary lens models and they are all consistent with a main sequence star lens located at ≤ 4 kpc from Earth. The microlensing models also indicate the presence of a bright blend, separated by ∆ ∼ 60 mas from the source, with ∼ (V, I) b = (18.5, 17.0) mag. Assuming that the blend suffers half of the total extinction towards the source, this object is compatible with a late F-dwarf at ∼ 2.5 kpc (Bessell & Brett 1988;Pecaut & Mamajek 2013), consistent with the lens properties derived from the galactic models analysis. The astrometric measurement made by Gaia at this position returns D = 2.56±0.79 kpc. Assuming this object to be the lens, we derived θ E ∼ 1.3 mas and 0.8M , also consistent with the previous estimations. Depending on the exact nature of the host, the lens companion is either a massive Jupiter or a low mass brown dwarf. Given their relative proper motion, µ rel = 5.5 mas/yr, the lens and source should be sufficiently separated to be observed via high-resolution imaging in about 10 years with 10-m class telescopes. This would provide the necessary additional information needed to confirm the exact nature of the lens, including the companion. Even though the physical nature of the host star cannot yet be firmly established, it is almost certain that the companion is located at the brown dwarf/planet mass boundary. The increasing number of reported discoveries of such objects, especially by microlensing surveys (see for example Bachelet et al. (2019) and references therein), provides important observational data which can be used to improve the theoretical framework underpinning planet formation. Indeed, there is more and more evidence that the critical mass to ignite deuterium (i.e., ∼ 13 M Jup ) does not represent a clear-cut limit (Chabrier et al. 2014). While there is compelling evidence that the two classes of objects are produced by different formation processes (Reggiani et al. 2016;Bowler et al. 2020), more observational constraints will be necessary in order to better appreciate the differences between them.
Similarly to MOA-2019-BLG-008, it can be expected that a fraction of events detected by the Roman microlensing survey will miss at least one mass-distance relation, i.e. θ E or π E . In this context, the Besançon and GalMod models can be particularly helpful in estimating the most likely parameters for the lens. Indeed, while Penny et al. (2019) and Terry et al. (2020) report some discrepancies between observations and their catalogs, these models are constantly upgraded to refine their predictions. In particular, the high accuracy astrometric measurements from Gaia will offer unique constraints on the proper motions and distances of stars up to the Galactic Bulge population at ∼ 8 kpc.
Software RAS and EB gratefully acknowledge support from NASA grant 80NSSC19K0291. YT and JW acknowledge the support of DFG priority program SPP 1992 "Exploring the Diversity of Extrasolar Planets" (WA 1047/11-1). KH acknowledges support from STFC grant ST/R000824/1. J.C.Y. acknowledges support from N.S.F Grant No. AST-2108414. Work by C.H. was supported by the grants of National Research Foundation of Korea (2019R1A2C2085965 and2020R1A4A2002885). This research has made use of NASA's Astrophysics Data System, and the NASA Exoplanet Archive. The work was partly based on data products from observations made with ESO Telescopes at the La Silla Paranal Observatory under programme ID 177.D-3023, as part of the VST Photometric Halpha Survey of the Southern Galactic Plane and Bulge (VPHAS+, www.vphas.eu). This work also made use of data from the European Space Agency (ESA) mission Gaia (https://www.cosmos.esa.int/gaia), processed by the Gaia Data Processing and Analysis Consortium (DPAC, https://www.cosmos.esa.int/web/gaia/dpac/consortium). Funding for the DPAC has been provided by national institutions, in particular the institutions participating in the Gaia Multilateral Agreement. CITEUC is funded by National Funds through FCT -Foundation for Science and Technology (project: UID/Multi/00611/2013) and FEDER -European Regional Development Fund through COMPETE 2020 -Operational Programme Competitiveness and Internationalization (project: POCI-01-0145-FEDER-006922). DMB acknowledges the support of the NYU Abu Dhabi Research Enhancement Fund under grant RE124. This research uses data obtained through the Telescope Access Program (TAP), which has been funded by the National Astronomical Observatories of China, the Chinese Academy of Sciences, and the Special Fund for Astronomy from the Ministry of Finance. This work was partly supported by the National Science Foundation of China (Grant No. 11333003, 11390372 and 11761131004 to SM). This research has made use of the KMTNet system operated by the Korea Astronomy and Space Science Institute (KASI) and the data were obtained at three host sites of CTIO in Chile, SAAO in South Africa, and SSO in Australia.The MOA project is supported by JSPS KAKENHI Grant Number JSPS24253004, JSPS26247023, JSPS23340064, JSPS15H00781, JP16H06287, and JP17H02871.

APPENDIX
A. NEW SDSS COLOR-RADIUS RELATION As discussed in the main text, the source magnitude in the g-band from the ROME survey is not well known, due to the low sampling of the lightcurve. But the source brightness in the r and i bands are well measured. Because Boyajian et al. (2014) does not provide a relation for these bands, we decided to collect the data and estimate these relations. As described in Boyajian et al. (2014) and avalaible on Simbad 12 , we used the magnitudes measruements from Boyajian et al. (2013) and the angular diameter measurements from various interferometers: the CHARA Array (Di Folco et al. 2004;Bigot et al. 2006;Baines et al. 2008Baines et al. , 2012Boyajian et al. 2012;Ligi et al. 2012;Bigot et al. 2011;?;Crepp et al. 2012;Bazot et al. 2011;Huber et al. 2012), the Palomar Testbed Interferometer (van Belle & von Braun 2009), the Very Large Telescope Interferometer (Kervella et al. 2003a(Kervella et al. ,b, 2004Di Folco et al. 2004;Thévenin et al. 2005;Chiavassa et al. 2012), the Sydney university Stellar Interformeter (?), the Narrabri Intensity Interferometer (Hanbury Brown et al. 1974), Mark III (Mozurkewich et al. 2003) and the Navy Prototype Oprical Interfereometer (Nordgren et al. 1999(Nordgren et al. , 2001. Then, we fitted the color-radius relation as: where X i is the considered color and i 0 is the de-reddened magnitude in the i-band. For stars that display several brightness measurements in the g, r and i bands, we used the mean as our final values, and the error was estimated from the sample variance with the (quadratic) addition of a 0.005 mag minimum error. We also used the error on the measured radius if available, and added quadratically the error on the observed i o magnitude. We explored several solutions using polynomials of different degrees and stopped as soon as the relative error on a i reached 1. Ultimately, we obtained: The data and best fit relations can be seen in Figure 8. As expected, the (r − i) o relation is less accurate (rms ∼ 0.05) than the (g − i) o relation (rms ∼ 0.04), especially for the coolest stars with (r − i) o ≥ 1 mag. But the accuracy is similar for a MOA-2019-BLG-009 source with (r − i) o ∼ 0.75 mag.