# OGLE-2018-BLG-0532Lb: Cold Neptune with Possible Jovian Sibling

Yoon-Hyun Ryu, Andrzej Udalski, Jennifer C. Yee, Matthew T. Penny, Weicheng Zang, Michael D. Albrow, Sun-Ju Chung, Andrew Gould, Cheongho Han, Kyu-Ha Hwang, et al.

► To cite this version:

Yoon-Hyun Ryu, Andrzej Udalski, Jennifer C. Yee, Matthew T. Penny, Weicheng Zang, et al. OGLE-2018-BLG-0532Lb: Cold Neptune with Possible Jovian Sibling. The Astronomical Journal, 2020, 160, 10.3847/1538-3881/abaa3f. insu-04836347

HAL Id: insu-04836347

https://insu.hal.science/insu-04836347v1

Submitted on 4 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>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_2.png" /></figure>

# OGLE-2018-BLG-0532Lb: Cold Neptune with Possible Jovian Sibling

Yoon-Hyun Ryu<fnref n="1" />, Andrzej Udalski<fnref n="2" />, Jennifer C. Yee<fnref n="3" />, Matthew T. Penny<fnref n="4" />, Weicheng Zang<fnref n="5" />, Michael D. Albrow<fnref n="6" />, Sun-Ju Chung<fnref n="1" />,<fnref n="7" />, Andrew Gould<fnref n="4" />,<fnref n="8" />, Cheongho Han<fnref n="9" />, Kyu-Ha Hwang<fnref n="1" />, Youn Kil Jung<fnref n="1" />, In-Gu Shin<fnref n="1" />, Dong-Joo Lee<fnref n="1" />, Yongseok Lee<fnref n="1" />,<fnref n="11" />, Byeong-Gon Park<fnref n="1" />,<fnref n="7" />, Richard W. Pogge<fnref n="4" />

Yossi Shvartzvald<fnref n="10" />, Sang-Mok Cha<fnref n="1" />,<fnref n="11" />, Dong-Jin Kim<fnref n="1" />, Hyoun-Woo Kim<fnref n="1" />,<fnref n="12" />, Seung-Lee Kim<fnref n="1" />,<fnref n="7" />, Chung-Uk Lee<fnref n="1" />,<fnref n="7" />, Dong-Joo Lee<fnref n="1" />, Yongseok Lee<fnref n="1" />,<fnref n="11" />, Byeong-Gon Park<fnref n="1" />,<fnref n="7" />, Richard W. Pogge<fnref n="4" />

(KMTNet Collaboration),

Przemek Mróz<fnref n="2" />, Michał K. Szymański<fnref n="2" />, Jan Skowron<fnref n="2" />, Radek Poleski<fnref n="4" />, Igor Soszyński<fnref n="2" />, Paweł Pietrukowicz<fnref n="2" />, Szymon Kozłowski<fnref n="2" />, Krzysztof Ulaczyk<fnref n="13" />, Krzysztof A. Rybicki<fnref n="2" />, Patryk Iwanek<fnref n="2" />, Marcin Wrona<fnref n="2" />

(OGLE Collaboration),

and
Shude Mao<fnref n="5" />,<fnref n="14" />, Pascal Fouque<fnref n="15" />,<fnref n="16" />, Wei Zhu<fnref n="17" />, and Tianshu Wang<fnref n="5" />

(CFHT microlensing collaboration)

Received 2019 April 15; revised 2020 July 24; accepted 2020 July 27; published 2020 September 25

## Abstract

We report the discovery of the planet OGLE-2018-BLG-0532Lb, with very obvious signatures in the light curve that lead to an estimate of the planet-host mass ratio $q = M_{\text{planet}}/M_{\text{host}} \simeq 1 \times 10^{-4}$. Although there are no obvious systematic residuals to this double-lens/single-source (2L1S) fit, we find that $\chi^2$ can be significantly improved by adding either a third lens (3L1S, $\Delta\chi^2 = 81$) or second source (2L2S, $\Delta\chi^2 = 77$) to the lens-source geometry. After thorough investigation, we conclude that we cannot decisively distinguish between these two scenarios and therefore focus on the robustly detected planet. However, given the possible presence of a second planet, we investigate to what degree and with what probability such additional planets may affect seemingly single-planet light curves. Our best estimates for the properties of the lens star and the secure planet are a host mass $M \sim 0.25 M_\odot$, system distance $D_L \sim 1$ kpc, and planet mass $m_{p,1} = 8 M_\oplus$ with projected separation $a_{1,\perp} = 1.4$ au. However, there is a relatively bright $I = 18.6$ (and also relatively blue) star projected within $<50$ mas of the lens, and if future high-resolution images show that this is coincident with the lens, then it is possible that it is the lens, in which case, the lens would be both more massive and more distant than the best-estimated values above.

Unified Astronomy Thesaurus concepts: Gravitational microlensing (672)

## 1. Introduction

Based on the second detection of a Neptune-class planet beyond the snow line, Gould et al. (2006) had already suggested that such “cold Neptunes” are the most common type of planet. Sumi et al. (2010) then showed that because there were roughly equal numbers of Neptunes and Jupiters, despite the decreasing sensitivity with mass ratio, Neptunes must be more common than Jupiters. As microlensing-planet discoveries continued to accumulate and populate the ($\log q, \log s$) diagram (Figure 7 from Mróz et al. 2017), it became manifest that cold Neptunes are at least more numerous than cold planets of greater mass.

Here, $q$ is the planet-host mass ratio and $s$ is the planet-host projected separation normalized to the Einstein radius $\theta_E$. These ($\log q, \log s$) diagrams show that planet detections are roughly uniform over $-4.3 < \log q < -2$. Because higher-mass planets have larger caustics, and so larger cross sections for detection, a uniform rate of detection implies more planets at lower mass.

The same ($\log q, \log s$) diagrams show a sharp cutoff in detections at $\log q \simeq -4.3$, i.e., $q \simeq 5 \times 10^{-5}$. See Figure 3 of Udalski et al. (2018). This could in principle reflect a sharp cutoff in the existence of cold planets at lower masses, but it also might simply reflect a cutoff in sensitivity of present-day microlensing experiments.

<footnote n="1">¹Korea Astronomy and Space Science Institute, Daejon 34055, Republic of Korea</footnote>

<footnote n="2">²Warsaw University Observatory, Al. Ujazdowskie 4, 00-478 Warszawa, Poland</footnote>

<footnote n="3">³Center for Astrophysics|Harvard & Smithsonian, 60 Garden St., Cambridge, MA 02138, USA</footnote>

<footnote n="4">⁴Department of Astronomy, Ohio State University, 140 W. 18th Ave., Columbus, OH 43210, USA</footnote>

<footnote n="5">⁵Physics Department and Tsinghua Centre for Astrophysics, Tsinghua University, Beijing 100084, People's Republic of China</footnote>

<footnote n="6">⁶University of Canterbury, Department of Physics and Astronomy, Private Bag 4800, Christchurch 8020, New Zealand</footnote>

<footnote n="7">⁷Korea University of Science and Technology, Korea, (UST), 217 Gajeong-ro, Yuseong-gu, Daejeon, 34113, Republic of Korea</footnote>

<footnote n="8">⁸Max-Planck-Institute for Astronomy, Königstuhl 17, D-69117 Heidelberg, Germany</footnote>

<footnote n="9">⁹Department of Physics, Chungbuk National University, Cheongju 28644, Republic of Korea</footnote>

<footnote n="10">¹⁰IPAC, Mail Code 100-22, Caltech, 1200 E. California Blvd., Pasadena, CA 91125, USA</footnote>

<footnote n="11">¹¹School of Space Research, Kyung Hee University, Yongin, Kyeonggi 17104, Republic of Korea</footnote>

<footnote n="12">¹²Department of Astronomy and Space Science, Chungbuk National University, Cheongju 28644, Republic of Korea</footnote>

<footnote n="13">¹³Department of Physics, University of Warwick, Gibbet Hill Rd., Coventry, CV4 7AL, UK</footnote>

<footnote n="14">¹⁴National Astronomical Observatories, Chinese Academy of Sciences, A20 Datun Rd., Chaoyang District, Beijing 100012, People's Republic of China</footnote>

<footnote n="15">¹⁵CFHT Corporation, 65-1238 Mamalahoa Hwy., Kamuela, HI 96743, USA</footnote>

<footnote n="16">¹⁶Université de Toulouse, UPS-OMP, IRAP, Toulouse, France</footnote>

<footnote n="17">¹⁷Canadian Institute for Theoretical Astrophysics, University of Toronto, 60 St George St., Toronto, ON M5S 3H8, Canada</footnote>
However, two studies concluded that the first explanation is correct: cold-planet frequency reaches a peak at cold Neptunes and then declines toward lower masses. First, Suzuki et al. (2016) studied planets detected by MOA and compared these to the planet sensitivity of the MOA sample, as judged by a $\Delta\chi^2$ criterion. They concluded that the cold-planet frequency (as a function of log *q*) peaks around log *q* ~ -4. Suzuki et al. (2018a) found the mass-ratio distribution from the sample of Suzuki et al. (2016) was inconsistent with the predicted mass-ratio distribution from the core accretion theories. Second, Udalski et al. (2018) studied the complete sample of seven microlensing planets with well-measured mass ratios log *q* < -4. They developed a new "V/Vmax" method that is logically independent of the Suzuki et al. (2016) method. The data samples were also largely independent. Udalski et al. (2018) concluded that if the mass-ratio function in this low-mass-ratio regime is modeled as a single power law, then it was falling toward lower log *q* in this range, thus confirming the results of Suzuki et al. (2016). In particular, Udalski et al. (2018) found that OGLE-2017-BLG-1434Lb would have been detected (and been well characterized) even if it had been 30 times lighter than it actually is, i.e., log *q* ≈ -5.7.

Pascucci et al. (2018) derived the planet mass-ratio function of Kepler transiting planets by combining the planet-radius distribution and the empirical mass-radius relation. They found a break at log *q* ≈ -4.5. A similar break (or pileup) is also found by Wu (2019), who inferred the Kepler planet masses from the photoevaporation gap. Because microlensing planets are typically outside the snow line by a factor of 2 or more, while Kepler transiting planets are inside the snow line by a factor of 10, these results enable one to study whether (and how) the break in the planet mass-ratio function is strongly influenced by the presence (or absence) of icy material in the protoplanetary disk.

Jung et al. (2019) subsequently analyzed the 15 planets with well-determined mass ratios $q < 3 \times 10^{-4}$ and concluded from their approximately uniform distribution in log *q* that the break in the microlensing mass-ratio function is both quite low, $q_{break} \approx 5.6 \times 10^{-5}$ (i.e., log $q_{break} \approx -4.25$), and quite severe, to the point of approximating a cutoff. They also pointed toward possible evidence for a “pileup,” at or just above the break. This break point would be more consistent with the Kepler break found by Pascucci et al. (2018), although probably much more severe.

Additional detections of planets in this low-mass regime will be crucial for resolving the position and severity of the break in the microlensing mass-ratio function. Here we analyze OGLE-2018-BLG-0532, whose lensing system contains at least one planet, a cold Neptune with mass ratio $q = (1.0 \pm 0.2) \times 10^{-4}$, i.e., right at the somewhat arbitrary “low mass-ratio” boundary of Udalski et al. (2018) and near the upper end of the possible “pileup” noted by Jung et al. (2019). We note that Ranc et al. (2019) have reported a planet that also straddles this boundary, with $q = (1.00 \pm 0.17) \times 10^{-4}$. That is, both of these planets would meet the sample conditions of the Jung et al. (2019) study.

Moreover, OGLE-2018-BLG-0532 also contains evidence for a second planet, a Jovian-class body that would lie at a projected distance that is either 2.65 times greater or 2.65 times smaller than the Neptune mass-ratio planet. Fitting for the additional planet results in an improvement of $\Delta\chi^2 = 81$, which is formally very significant. At the same time, the

candidate planet does not result in any recognizable signatures, either in the original light curve or in the residuals to the single-planet fit.

Shin et al. (2015) studied the problem of such weak third-body signatures in apparently single-planet microlensing events. At the time of their study (and still today), there were only two cases of two-planet microlensing systems with unambiguous signatures for both planets: OGLE-2006-BLG-109 (Gaudi et al. 2008; Bennett et al. 2010) and OGLE-2012-BLG-0026 (Han et al. 2013). In both cases, the underlying events had high magnifications (normalized impact parameters $u_0 = 0.0035$ and $u_0 = 0.0095$, respectively), which made them especially sensitive to planets via their “central caustics” (Griest & Safizadeh 1998), and thus particularly sensitive to multiple planets in the same system (Gaudi et al. 1998). In fact, both systems followed the “factoring” (or “superposition”) of caustic signatures already predicted by Han et al. (2001).

However, Shin et al. (2015) considered that if there were additional planets in a given system with one detected planet, they would be at least as likely to give rise to weak signatures that may have escaped notice as they would be to give rise to strong signatures that were obvious, at least in the residuals to the single-planet fit. This led them to search for additional bodies in eight apparently single-planet systems in events of high or moderately-high magnification, OGLE-2005-BLG-071 (Udalski et al. 2005), OGLE-2005-BLG-169 (Gould et al. 2006), MOA-2007-BLG-400 (Dong et al. 2009a), MOA-2008-BLG-310 (Janczak et al. 2010), MOA-2009-BLG-319 (Miyake et al. 2011), MOA-2009-BLG-387 (Batista et al. 2011), MOA-2010-BLG-477 (Bachelet et al. 2012), and MOA-2011-BLG-293 (Yee et al. 2012). These have, respectively, normalized impact parameters $u_0 = (0.023, 0.00012, 0.0005, 0.0030, 0.0062, 0.08, 0.0034, 0.0035)$.

Shin et al. (2015) found that allowing for a third body led to $\Delta\chi^2$ improvements of 143, 78, and 50 for events MOA-2009-BLG-319, MOA-2008-BLG-310, and MOA-2010-BLG-477, respectively. All other events had $\Delta\chi^2 < 30$. Given that none of these events met their adopted threshold of $\Delta\chi^2 > 500$, they only set upper limits on additional planets in all cases. They showed only one two-planet fit, i.e., for MOA-2009-BLG-319, which had the highest $\Delta\chi^2$. As we will show is also the case for OGLE-2018-BLG-0532, the single-planet fit did not display any particularly noticeable residuals, and the two-planet fit did not result in any obvious improvements.

Subsequently, Suzuki et al. (2018b) showed that the OGLE-2014-BLG-1722 light curve has two unambiguous deviations from a Paczyński (1986) fit that are well explained by a 3L1S model with two planets. One deviation is a dip and so can only be due to a planet, but the other is a smooth bump, which could also in principle be caused by a second source. This 2L2S model is disfavored by only $\Delta\chi^2 = 5.7$, but after taking consideration of auxiliary information, Suzuki et al. (2018b) concluded that the 3L1S model is favored by $3.1\sigma$ and so consider that this is the “likely” interpretation. Hence, this is currently the best case for a relatively weak detection of a second planet, although the 3L1S interpretation remains less than absolutely secure.

For OGLE-2018-BLG-0532, the smooth character of the residuals from the 2L1S fit also implies that one should test for 2L2S (in addition to 3L1S) solutions. Although we find that 3L1S is preferred by $\Delta\chi^2 = 4$ relative to 2L2S, we consider that the evidence is too weak to claim the existence of a second
<figure><img src="image_3.png" /><figcaption>**Figure 1.** Light curve and best-fit 3L1S (star plus two planets) model of OGLE-2018-BLG-0532, with a zoom of the perturbation region in the upper set of three panels, which also shows the best-fit 2L2S and 2L1S models. The other degenerate models from Tables 2, 3, 4, and 1 are indistinguishable from these. The broad flat trough lasting ~1.5 days, together with the adjacent peaks is due to a mass ratio $q = 1 \times 10^{-4}$ planet. The second planet, which is about 25 times more massive, does not give rise to obvious signatures. The panels immediately below the light-curve panels show the residuals to this fit, while the bottom panels show residuals to the 2L2S and 2L1S fits. The differences between these panels are not obvious to the eye. The 3L1S fit is favored by $\Delta\chi^2 = 81$ relative to the 2L1S fit and by $\Delta\chi^2 = 16$ relative to the 2L2S fit (without orbital motion).</figcaption></figure>

planet. Nevertheless, the possibility of such a planet prompts us to investigate how such “weak detections” are likely to appear in the data.

## 2. Observations

OGLE-2018-BLG-0532 is at (R.A., decl.) = (17:59:56.02, -28:59:51.9) corresponding to (l, b) = (1.54, -2.73). It was announced to the community by the Optical Gravitational Lensing Experiment (OGLE; Udalski et al. 2015) Early Warning System (Udalski et al. 1994; Udalski 2003) at UT 13:54 on 2018 April 9. OGLE observed this field (BLG505) in 2018 at a cadence $\Gamma = 1 \text{ hr}^{-1}$ from its 1.3 m telescope with 1.4 deg$^2$ camera at Las Campanas Observatory in Chile. These

data show a clear anomaly, but are adversely impacted by long weather gaps during the anomaly.

The Korea Microlensing Telescope Network (KMTNet; Kim et al. 2016) also observed this event from its three identical 1.6 m telescopes at Cerro Tololo Inter-American Observatory (CTIO; Chile, KMTC), South African Astronomical Observatory (SAAO; South Africa, KMTS), and Siding Spring Observatory (SSO; Australia, KMTA), each equipped with identical 4 deg$^2$ cameras. The event lies in three overlapping KMTNet fields, BLG03, and BLG43, with a combined cadence of $\Gamma = 4 \text{ hr}^{-1}$ from KMTC and $\Gamma = 6 \text{ hr}^{-1}$ from KMTS and KMTA. In fact, there are some observations from KMTC BLG04 as well, with a cadence of $\Gamma = 1 \text{ hr}^{-1}$. The event was rediscovered at the end of the season by the KMTNet eventfinder (Kim et al. 2018) as KMT-2018-BLG-1161.
The great majority of observations were carried out in the I band with occasional V-band observations made solely to determine source colors. All reductions for the light-curve analysis were conducted using variants of image subtraction (Tomaney & Crotts 1996; Alard & Lupton 1998), either DIA (Woźniak 2000) or pySIS (Albrow et al. 2009). The errors were renormalized so that the $\chi^2$ per degree of freedom (dof) at each observatory is approximately equal to unity.

### 3. Analysis

#### 3.1. Search for 2L1S Solutions

Figure 1 shows the OGLE and KMTNet data, together with the final best-fit triple-lens/single-source (3L1S) model. The duration of noticeable magnification is quite long, $\gtrsim 100$ days, and the peak of the main event is relatively sharp. Together, these characteristics would normally indicate a high-magnification event, $A_{max} \gg 1$. However, the flux at the peak is only a factor of $F_{peak}/F_{base} \sim 3.6$ above baseline. These two indicators can only be reconciled if the event is heavily blended, $f_b \gg f_s$, where $f_s$ and $f_b$ are the source flux and the blended flux, respectively, and where $F_{base} = f_s + f_b$.

The major anomaly, which is shown in the upper (zoomed) panel of Figure 1, is comprised of a trough that lasts $\sim 1.5$ days, flanked by two spikes, with the post-trough spike being particularly sharp. Such troughs can only be produced by a minor-image perturbation. The unperturbed light curve is produced by two images that, according to Fermat's principle, occur at extrema of the time-delay surface. The larger (major) image, on one side of the primary lens, is at a local minimum of this surface and is extremely stable. The smaller (minor) image, on the opposite side, is at a saddle point and is highly unstable. Hence, it can be virtually annihilated by a planet sitting at or near the position of the minor image. For high-magnification events $A_+ + A_- = A \gg 1$, the ratio of these magnifications $A_+/A_- = (A+1)/(A-1) \to 1$. Hence, up to half the light can be briefly eliminated as a planet passes near this image. On the flanking sides of such troughs are two caustics (lines of formally infinite magnification), which continue as extended ridges in magnification beyond the finite extent of the caustics themselves. Thus, the form of the observed perturbation is exactly as expected for such a geometry, and there are no other known geometries that can generate this light-curve morphology.

Notwithstanding this qualitative analysis, we begin by conducting a systematic grid search for 2L1S models over a range of $(s, q)$ geometries. We hold these two parameters fixed at the grid points, while five others $(t_0, u_0, t_E, \rho, \alpha)$ are allowed to vary in a Markov Chain Monte Carlo (MCMC). The first three (Paczyński 1986) parameters $(t_0, u_0, t_E)$ are seeded at values derived from a 1L1S fit to the light curve with the anomaly removed. The angle $\alpha$ between the binary axis and lens-source relative proper motion $\mu_{rel}$ is seeded at six equally spaced points around the unit circle, while $\rho = \theta_*/\theta_E$, i.e., the source radius normalized to the Einstein radius, is seeded at $\rho = 5 \times 10^{-4}$. As usual, there are two flux parameters $(f_s, f_b)_i$ for each observatory $i$, so that the predicted flux is given by $F_i(t) = f_{s,i}A(t; s, q, \alpha, t_0, u_0, t_E, \rho) + f_{b,i}$. During the brief intervals when the source passes over or very close to the caustic, it is partially resolved by the caustic, and so the source brightness profile becomes relevant. We characterize this by a

<table><caption>**Table 1**
Best-fit Solutions for 2L1S Models</caption><thead><tr><th rowspan="2">Parameters</th><th rowspan="2">Standard</th><th colspan="2">Parallax Models</th></tr><tr><th>u<sub>0</sub> > 0</th><th>u<sub>0</sub> < 0</th></tr></thead><tbody><tr><td>χ<sup>2</sup>/dof</td><td>9240.291/9124</td><td>9169.300/9122</td><td>9169.877/9122</td></tr><tr><td>t<sub>0</sub> (HJD')</td><td>8219.579 ± 0.004</td><td>8219.579 ± 0.004</td><td>8219.579 ± 0.004</td></tr><tr><td>u<sub>0</sub> (10<sup>-3</sup>)</td><td>9.341 ± 0.371</td><td>8.295 ± 0.353</td><td>-7.873 ± 0.369</td></tr><tr><td>t<sub>E</sub> (days)</td><td>118.327 ± 4.709</td><td>133.269 ± 5.849</td><td>140.226 ± 5.962</td></tr><tr><td>s</td><td>1.014 ± 0.0005</td><td>1.012 ± 0.001</td><td>1.013 ± 0.001</td></tr><tr><td>q (10<sup>-4</sup>)</td><td>1.389 ± 0.053</td><td>1.015 ± 0.057</td><td>0.962 ± 0.056</td></tr><tr><td>α (rad)</td><td>-0.431 ± 0.002</td><td>-0.465 ± 0.005</td><td>0.465 ± 0.005</td></tr><tr><td>ρ (10<sup>-4</sup>)</td><td>3.154 ± 0.099</td><td>2.734 ± 0.105</td><td>2.680 ± 0.101</td></tr><tr><td>π<sub>E,N</sub></td><td>...</td><td>-0.785 ± 0.105</td><td>0.765 ± 0.097</td></tr><tr><td>π<sub>E,E</sub></td><td>...</td><td>-0.104 ± 0.018</td><td>-0.054 ± 0.015</td></tr><tr><td>f<sub>S</sub></td><td>0.0201 ± 0.0008</td><td>0.0181 ± 0.0008</td><td>0.0171 ± 0.0008</td></tr><tr><td>f<sub>B</sub></td><td>0.7764 ± 0.0008</td><td>0.7787 ± 0.0008</td><td>0.7797 ± 0.0008</td></tr><tr><td>t<sub>*</sub> (days)</td><td>0.037 ± 0.001</td><td>0.036 ± 0.001</td><td>0.038 ± 0.002</td></tr></tbody></table>

Note. $t_* \equiv \rho t_E$ is not an independent quantity.

linear limb-darkening coefficient $\Gamma_I = 0.50$ (Claret 2000), in accordance with the source-type determination of Section 4.2.

As anticipated, we find only one 2L1S solution, the parameters for which are given in Table 1. Before continuing, we remark on the extreme level of blending $f_b/f_s \sim 40$. That is, despite the fact that the baseline appears relatively bright, $I_{base} = 18.24$ in the OGLE-IV system (see Section 4.1.2), the source is extremely faint, $I_s \sim 22.25$.

We transform the instrumental OGLE-IV system (Udalski et al. 2015) to the standard Johnson-Cousins system using the calibration constants:

$$ \begin{aligned} (V - I)_{\text{calib}} &= 0.93[(V - I)_{\text{OGLE-IV}} + 0.309]; \\ I_{\text{calib}} &= I_{\text{OGLE-IV}} - 0.005 - 0.008(V - I)_{\text{calib}}. \end{aligned} \quad (1) $$

Thus, in the standard Johnson-Cousins system, $I_{base,standard} = 18.23$ and $I_{s,standard} = 22.23$.

Figure 2 shows the caustic geometry for the 2L1S solution presented here, as well as for additional solutions presented in Sections 3.5 and 3.6. It shows that the pronounced dip near the peak is due to the source transiting the “back end” of a resonant caustic.

#### 3.2. 2L1S Solutions with Parallax

In spite of the faintness of the source (hence, high fractional errors $\sigma(F)/f_s$), the long Einstein timescale $t_E \sim 120$ days implies that the microlens parallax may well be measurable. In particular, Smith et al. (2003) showed that the strength of the parallax signal basically scales as $t_E^4$. Moreover, in the case of caustic crossing events, the parallax signal can be augmented by information from the caustic crossing (An & Gould 2001; Muraki et al. 2011; Sumi et al. 2016). This introduces two additional parameters $(\pi_{E,N}, \pi_{E,E})$, i.e., the components in equatorial coordinates of the “vector microlens parallax” $\pi_E$ (Gould 1992, 2000),

$$ \pi_E = \frac{\pi_{rel}}{\theta_E} \frac{\mu_{rel}}{\mu_{rel}}, \quad (2) $$

where $\theta_E$ is the Einstein radius,

$$ \theta_E = \sqrt{\kappa M \pi_{rel}}; \quad \kappa \equiv \frac{4G}{c^2 \text{au}} \approx 8.14 \frac{\text{mas}}{M_{\odot}}, \quad (3) $$
<figure><img src="image_4.png" /><figcaption>**Figure 2.** Caustic geometries for the 2L1S (Section 3.1), 2L2S (Section 3.6), and 3L1S (Section 3.5) solutions. The left panels show zooms of these geometries in the neighborhood of the host. In all cases, the main feature, a pronounced dip near the peak, is caused by the source transiting the magnification trough at the “back end” of a resonant caustic. In the 2L2S solution, there is a second source (red line) passing well away from the caustic that accounts for the residuals. The two sources are shown at the same time. In the bottom two panels, a second planet (3L1S) distorts the caustic near the central cusp (left panels), which then accounts somewhat better for the residuals to the 2L1S fit.</figcaption></figure>

$M$ is the total mass of the lens, and $\pi_{\text{rel}}$ is the relative lens-source parallax.

Because the parallax effect (due to Earth's orbital motion) can be degenerate with orbital motion of the lens (Batista et al. 2011; Skowron et al. 2011), one must also introduce two additional parameters representing linearized orbital motion in order to test for such correlations. We adopt $d\alpha/dt$ and $ds/dt$ for the time rates of change of the binary's orientation and separation, respectively. However, we find no significant correlation between the orbital motion and parallax parameters. Instead, the orbital parameters become “attracted” to solutions that model minor fluctuations in the 2017 (i.e., baseline) data as being due to a planetary caustic that fortuitously “orbited” to the source position during 2017. This leads to small $\chi^2$ “improvements” that are entirely spurious. Because the orbital-motion parameters are not in fact correlated with the parallax

parameters, which was the original reason for introducing them, we suppress the orbital-motion parameters.

For sources lying near the ecliptic, which includes essentially all microlensing toward the Galactic bulge, one must check for the “ecliptic degeneracy,” which approximately takes $(u_0, \pi_E, \alpha, d\alpha/dt) \rightarrow -(u_0, \pi_E, \alpha, d\alpha/dt)$. We indeed find two such solutions with very similar $\chi^2$, which are both shown in Table 1. We also search for jerk-parallax solutions (Gould 2004), which would approximately take $\pi_{E,N} \rightarrow -\pi_{E,N}$, while leaving the other parameters approximately unchanged, but this search does not lead to viable solutions.

We note that most parameters are fairly similar between the “parallax” and “standard” models, except of course the parallax, which is newly introduced. The $\Delta\chi^2 = 71$ is highly significant, and corresponding to this, the parallax amplitude $\pi_E$ is measured to about 13%. By far the largest change, apart from
<figure><img src="image_5.png" /><figcaption>**Figure 3.** $\chi^2$ distribution of xallarap model according to the orbital period.</figcaption></figure>

the parallax, is the mass ratio, which drops by roughly 25% in the parallax solutions.

One point of possible concern is that $\pi_{E,E}$ is very close to zero, with very small measurement errors, while $\pi_{E,N}$ is much larger. This might be regarded as worrisome because Earth's acceleration is essentially east, and this induces an asymmetry in the light curve, which is quite easy to measure because it is not strongly correlated with any other parameters. Hence, if the vector parallax were actually close to zero, $\pi_E \sim 0$, then this would be robustly reflected in $\pi_{E,\|} \simeq \pi_{E,E}$ but might be masked by subtle long-term systematics in the light curve that gave rise to a spurious $\pi_{E,\perp} \simeq \pi_{E,N}$.

However, first, it will be straightforward to show that the parallax must obey $\pi_E \gtrsim 0.2$ (see Section 6.3). We do not discuss these arguments in detail here to avoid repeating them later. Second, while $|\pi_{E,E}|$ is small compared to $|\pi_{E,N}|$, it is strongly inconsistent with zero, and quite consistent with the lower limit on $\pi_E$ just given. Therefore, at this point there is no reason to doubt the parallax measurement. Nevertheless, we will return in Section 6.3 to the question of whether it could be several $\sigma$ lower than the values reported in Table 1.

### 3.3. Xallarap Analysis

The parallax effect can be mimicked by the orbital motion of the source star with its companion, which is called the “xallarap” effect. Therefore, we check whether there is a real xallarap effect that is being misinterpreted as a spurious

parallax effect. The xallarap model (with assumed circular orbit) requires five additional parameters, which are the xallarap-vector components in the north and east directions ($\xi_{E,N}$ and $\xi_{E,E}$), the orbital period $P$, and the ecliptic longitude ($l_{eclip}$) and latitude ($b_{eclip}$) of the binary-source orbit.

As discussed in Poindexter et al. (2005), if the xallarap model perfectly reproduces the parallax effect, we expect the best-fit solution to have $P = 1$ yr and ($l_{eclip}, b_{eclip}$) equal to the ecliptic coordinates of the event. However, because of the three additional degrees of freedom in this model relative to the parallax models, this model also has more freedom to fit other unrelated effects, such as low-level systematics.

We perform two tests of xallarap models. First, we fit models at several fixed periods $P = (0.2, 0.4, 0.8, 1.0, 1.2, 1.4, 1.6, 1.8)$ yr. We find the local minima using MCMC and obtain an overall $\chi^2$ minimum near $P \sim 1$ yr (see Figure 3). This is what we expect if the xallarap effect matched the reflex effect of parallax due to the motion of Earth around the Sun.

As a second test, we examine the best-fit coordinates for $P = 1$ yr. These coordinates are offset from the true coordinates by $\sim 27^\circ$ and improve the fit by $\Delta\chi^2 = 19$. In principle, a xallarap signal could have any set of coordinates over $\pi$ sr, taking account of the four-fold symmetry of xallarap solutions. Thus, there is only a $\sim 20\%$ probability of lying within $27^\circ$ of the true coordinates.

Therefore, because the best-fit period is the same as Earth's and the coordinates are close to the true coordinates, we
conclude that parallax is the dominant effect. However, the
offset in the coordinates is an indication that xallarap or
systematics could be present at a modest level. As mentioned
above, we will revisit the robustness of the parallax measure-
ment in Section 6.3.

3.4. Decision to Investigate 3L1S Solutions

The residuals to the 2L1S model shown in Figure 1 do not exhibit pronounced deviations, and therefore one would not under normal circumstances search for a third body. Our decision to carry out such a 3L1S investigation was prompted by “accidental” developments in the course of the 2L1S investigation, i.e., apparent “problems” in the 2L1S solution that were all eventually resolved. For completeness, we describe this process in the Appendix, but we do not divert the reader’s attention with the details here.

3.5. *Search for 3L1S Solutions*

We begin by considering static 3L1S models with 10 geometric parameters, i.e., the seven “standard” 2L1S parameters (with ($q,s$) renamed ($q_1, s_1$)), plus three additional parameters, ($q_2, s_2, \psi$). These are, respectively, the mass ratio of the second companion to the primary, the normalized separation between these, and the angle between the primary–secondary and primary–tertiary directions. We conduct this search by simply adding a new component with $q_2 = 0.3$ and $s_2 = 0.05$ and at several different values of $\psi$ and then using MCMC to search for a local minimum. Such $s_2 \ll 1$ models are in the extreme Chang–Refsdal regime, in which the pseudo-shear<fnref n="18" /> $\gamma_2^c \equiv q_2s_2^2$ is an approximate invariant. Therefore, one expects that (if the choice of $\psi$ approximately conforms to the physical configuration of the triple lens) the $\gamma_2^c$ parameter-combination will quickly approach the true value in the MCMC, while the subsequent disambiguation between $q_2$ and $s_2$ at fixed $\gamma_2^c$ will proceed much more slowly. In fact, we found that the first expectation was confirmed but the second was not. Instead, the MCMC became “stuck” and could not proceed toward the ($q_2, s_2$) minimum, probably due to the weakness of the $\chi^2$ gradient. We therefore substituted $s_2 \rightarrow \gamma_2^c$ as an MCMC variable, after which the minimum was found relatively quickly.

This search yields solutions with $q_2 \ll 1$ and $s_2 < 0.5$
similar model light curve from the corresponding wide solution
($q_2, s_2) \rightarrow (q_2, s_2^{-1}$) (Griest & Safizadeh 1998; Han et al. 2013;
Song et al. 2014). We find that this yields a model with very
similar $\chi^2$.

We then add the two parallax parameters ($\pi_{E,N}$, $\pi_{E,E}$). After finding the best fit, we attempt to add orbital motion (to first one, then the other companion). However, as with the 2L1S solutions, the $\chi^2$ only improves when the planetary caustic becomes “attached” to noise features in the baseline. Hence, we again remove the orbital-motion parameters.

The best fit among these four solutions (close, wide) × ($u_0 < 0, u_0 > 0$) is shown in Figure 1, together with the residuals of the data from this model. By eye, the improvement relative to the 2L1S residuals looks modest.

The four solutions are shown in Table 2. Comparison of this
table with Table 1 shows that the 3L1S and 2L1S solutions are
overall very similar (except for the addition of three parameters
due to the second planet). This additional planet is about 25
times more massive than the first planet, which is qualitatively
similar to the ratio of the masses of Jupiter and Neptune. The
addition of this planet improves the fit by $\Delta\chi^2 = 81$. This is
formally very significant, but it would be by far the lowest $\Delta\chi^2$
to justify an “officially accepted” planet. Therefore, some
discussion is required to assess the reality of this putative
planet. See Section 6.

3.6. *Search for 2L2S Solutions*

As mentioned in Section 3.4, there are no obvious systematic
residuals to the 2L1S solution. This means that the evidence for
the third body is not an additional caustic crossing (unlike, e.g.,
OGLE-2006-BLG-109 and OGLE-2012-BLG-0026). When-
ever the evidence for an additional body lacks such caustic
features, one must always test whether it can be generated by
an additional source rather than an additional lens, i.e., in this
case 2L2S. We therefore begin by adding three parameters to
the 2L1S geometry, which add a second source with fixed
projected offset from the first source. That is, we rename
($t_0, u_0) \rightarrow (t_{0,1}, u_{0,1})$, introduce a second pair ($t_{0,2}, u_{0,2}$),
and then also introduce the flux ratio of the second source to the
first in I band, $q_{F,I}$. This yields a substantial improvement
$\Delta\chi^2 = 124$ relative to the standard 2L1S solution. We then
add the two parallax parameters, ($\pi_{E,N}$, $\pi_{E,E}$), which yields an
additional $\Delta\chi^2 = 12$ improvement. Normally, such a modest
improvement would not be considered sufficient for a reliable
parallax measurement. However, in the present case, we must
include these parameters in order to allow a fair comparison
with the 3L1S solutions. Comparison of Tables 3 and 2 shows
that the 3L1S solution is preferred by $\Delta\chi^2 = 16$ relative to the
2L2S solution.

3.6.1. 2L2S Solutions with Source Orbital Motion

The projected separation of the two sources in the static 2L2S models in Table 3 is ~0.02θ<sub>E</sub>, corresponding to about 0.25 au at the distance of the Galactic bulge and hence to periods possibly as short as ~45 days. Therefore, we should also allow for orbital motion of the binary source, particularly because, as just noted, the 3L1S solution is only preferred by $\Delta\chi^2 = 16$ relative to the static 2L2S model. For this purpose, we assume circular orbits, and (based on the source magnitude 7.4 mag below the clump, and source-flux ratio $q_{F,I} \sim 0.125$) we adopt source masses of 0.5 $M_\odot$ and 0.3 $M_\odot$ for the primary and secondary, respectively. We adopt a source distance of 8.0 kpc, which means that the source period is determined from the (3D) source separation by Kepler's Third Law. Then, to describe the time-varying source positions, we replace the previous four parameters ($t_{0,1}, u_{0,1}, t_{0,2}, u_{0,2}$) by six parameters ($t_{0,c}, u_{0,c}, s_c, \psi, \alpha_s, \delta_s$). Here ($t_{0,c}, u_{0,c}$) describe the Einstein ring position of ($t_0, u_0$) of the source center of mass, ($\alpha_s, \delta_s$) describe the orientation of the source-orbit angular momentum vector, $s_s$ is the source (3D) separation in units of $\theta_E$, and $\psi$ is the phase of the orbit at $t_{0,c}$.

Note that while this 2L2S-orbital-motion model has the same
number of sources and lenses as the “xallarap” model
considered in Section 3.3, the two physical systems that
underlie these models differ substantially. Most importantly,

<footnote n="18">¹⁸ Strictly speaking, the shear refers to the wide regime $s_2 \gg 1$, where it is $\gamma_2^w = q_2/s_2^2$. Based on the close-wide degeneracy, model pairs with $\gamma_2^c \sim \gamma_2^w$ should produce similar light curves. For simplicity, we therefore call $\gamma_2^c$ the “pseudo-shear.”</footnote>
<table><caption>**Table 2**
Best-fit Solutions for 3L1S Models</caption><thead><tr><th rowspan="2">Parameters</th><th colspan="2">Close</th><th colspan="2">Wide</th></tr><tr><th>u<sub>0</sub> > 0</th><th>u<sub>0</sub> < 0</th><th>u<sub>0</sub> > 0</th><th>u<sub>0</sub> < 0</th></tr></thead><tbody><tr><td>χ<sup>2</sup>/dof</td><td>9088.321/9119</td><td>9092.882/9119</td><td>9090.243/9119</td><td>9093.705/9119</td></tr><tr><td>t<sub>0</sub> (HJD')</td><td>8219.590 ± 0.009</td><td>8219.571 ± 0.008</td><td>8219.479 ± 0.028</td><td>8219.493 ± 0.028</td></tr><tr><td>u<sub>0</sub> (10<sup>-3</sup>)</td><td>8.228 ± 0.613</td><td>-8.893 ± 0.645</td><td>7.133 ± 0.503</td><td>-7.428 ± 0.709</td></tr><tr><td>t<sub>E</sub> (days)</td><td>139.317 ± 9.735</td><td>128.530 ± 9.344</td><td>151.039 ± 8.858</td><td>145.361 ± 11.011</td></tr><tr><td>s<sub>1</sub></td><td>1.013 ± 0.001</td><td>1.013 ± 0.001</td><td>1.011 ± 0.001</td><td>1.011 ± 0.001</td></tr><tr><td>q<sub>1</sub> (10<sup>-4</sup>)</td><td>0.975 ± 0.120</td><td>1.187 ± 0.133</td><td>0.927 ± 0.100</td><td>0.963 ± 0.168</td></tr><tr><td>α (rad)</td><td>-0.477 ± 0.007</td><td>0.475 ± 0.007</td><td>-0.478 ± 0.007</td><td>0.474 ± 0.009</td></tr><tr><td>ρ (10<sup>-4</sup>)</td><td>2.793 ± 0.240</td><td>3.310 ± 0.255</td><td>2.905 ± 0.198</td><td>2.884 ± 0.308</td></tr><tr><td>s<sub>2</sub></td><td>0.364 ± 0.043</td><td>0.406 ± 0.041</td><td>2.656 ± 0.375</td><td>2.634 ± 0.396</td></tr><tr><td>q<sub>2</sub> (10<sup>-3</sup>)</td><td>3.081 ± 0.871</td><td>2.231 ± 0.712</td><td>2.456 ± 1.059</td><td>2.208 ± 1.022</td></tr><tr><td>ψ (rad)</td><td>-0.036 ± 0.047</td><td>0.144 ± 0.042</td><td>-0.086 ± 0.047</td><td>0.084 ± 0.055</td></tr><tr><td>π<sub>E,N</sub></td><td>-0.675 ± 0.129</td><td>0.473 ± 0.157</td><td>-0.630 ± 0.100</td><td>0.625 ± 0.227</td></tr><tr><td>π<sub>E,E</sub></td><td>-0.105 ± 0.018</td><td>-0.041 ± 0.016</td><td>-0.087 ± 0.016</td><td>-0.038 ± 0.015</td></tr><tr><td>f<sub>S</sub></td><td>0.0175 ± 0.0013</td><td>0.0189 ± 0.0013</td><td>0.0161 ± 0.0011</td><td>0.0168 ± 0.0015</td></tr><tr><td>f<sub>B</sub></td><td>0.7790 ± 0.0013</td><td>0.7776 ± 0.0013</td><td>0.7804 ± 0.0011</td><td>0.7799 ± 0.0015</td></tr><tr><td>t<sub>*</sub> (days)</td><td>0.039 ± 0.002</td><td>0.043 ± 0.002</td><td>0.044 ± 0.002</td><td>0.042 ± 0.002</td></tr></tbody></table>

Note. $t_* \equiv \rho t_E$ is not an independent quantity.

<table><caption>**Table 3**
Best-fit Solutions for 2L2S Models</caption><thead><tr><th rowspan="2">Parameters</th><th rowspan="2">Standard</th><th colspan="2">Parallax</th></tr><tr><th>u<sub>0</sub> > 0</th><th>u<sub>0</sub> < 0</th></tr></thead><tbody><tr><td>χ<sup>2</sup>/dof</td><td>9116.315/9121</td><td>9104.215/9119</td><td>9106.616/9119</td></tr><tr><td>t<sub>0.1</sub> (HJD')</td><td>8219.598 ± 0.005</td><td>8219.594 ± 0.005</td><td>8219.594 ± 0.005</td></tr><tr><td>u<sub>0.1</sub> (10<sup>-3</sup>)</td><td>7.832 ± 0.429</td><td>5.926 ± 0.392</td><td>-5.557 ± 0.438</td></tr><tr><td>t<sub>E</sub> (days)</td><td>131.607 ± 6.148</td><td>173.809 ± 8.555</td><td>182.493 ± 11.562</td></tr><tr><td>s</td><td>1.014 ± 0.0006</td><td>1.013 ± 0.0006</td><td>1.014 ± 0.0006</td></tr><tr><td>q (10<sup>-4</sup>)</td><td>1.110 ± 0.060</td><td>0.705 ± 0.056</td><td>0.660 ± 0.064</td></tr><tr><td>α (rad)</td><td>-0.403 ± 0.004</td><td>-0.436 ± 0.007</td><td>0.431 ± 0.007</td></tr><tr><td>ρ<sub>1</sub> (10<sup>-4</sup>)</td><td>2.747 ± 0.112</td><td>2.149 ± 0.122</td><td>1.986 ± 0.135</td></tr><tr><td>t<sub>0.2</sub> (HJD')</td><td>8217.655 ± 0.598</td><td>8219.763 ± 0.575</td><td>8220.006 ± 0.342</td></tr><tr><td>u<sub>0.2</sub></td><td>0.048 ± 0.007</td><td>0.026 ± 0.006</td><td>-0.023 ± 0.004</td></tr><tr><td>qF, I</td><td>0.137 ± 0.017</td><td>0.119 ± 0.018</td><td>0.130 ± 0.017</td></tr><tr><td>π<sub>E,N</sub></td><td>...</td><td>-0.557 ± 0.095</td><td>0.527 ± 0.098</td></tr><tr><td>π<sub>E,E</sub></td><td>...</td><td>-0.072 ± 0.020</td><td>-0.041 ± 0.015</td></tr><tr><td>f<sub>S</sub></td><td>0.019 ± 0.0011</td><td>0.014 ± 0.0009</td><td>0.013 ± 0.0010</td></tr><tr><td>f<sub>B</sub></td><td>0.777 ± 0.0010</td><td>0.782 ± 0.0009</td><td>0.783 ± 0.0010</td></tr><tr><td>t<sub>*</sub> (days)</td><td>0.036 ± 0.002</td><td>0.037 ± 0.002</td><td>0.036 ± 0.002</td></tr></tbody></table>

Note. $t_* \equiv \rho t_E$ is not an independent quantity.

the xallarap model sought to explain the apparent parallax effect by source orbital motion, whereas the 2L2S-orbital-motion model explains these long-term effects by parallax, while seeking to explain shorter-term deviations by an additional xallarap effect. Second, the companion to the source is dark in the xallarap model, while it is luminous in the 2L2S-orbital-motion model.

We find that inclusion of orbital motion further improves the fit by $\Delta\chi^2 = 12$ for two additional degrees of freedom. See Table 4. Hence, comparing the best 2L2S solution with the best 3L1S solution, one sees that the former is disfavored by $\Delta\chi^2 = 3.8$ despite two fewer degrees of freedom. By formal criteria under the assumption of Gaussian errors, this would imply that the 3L1S solution was mildly favored. However, this small $\chi^2$ difference is on the order of or smaller than typical systematic errors in microlensing. Hence the two classes of

<table><caption>**Table 4**
Best-fit Solutions for 2L2S Models with Source Orbital Motion</caption><thead><tr><th>Parameters</th><th>u<sub>0</sub> > 0</th><th>u<sub>0</sub> < 0</th></tr></thead><tbody><tr><td>χ<sup>2</sup>/dof</td><td>9092.077/9117</td><td>9094.720/9117</td></tr><tr><td>t<sub>0,c</sub> (HJD')</td><td>8219.683 ± 0.065</td><td>8219.461 ± 0.076</td></tr><tr><td>u<sub>0,c</sub> (10<sup>-2</sup>)</td><td>1.353 ± 0.065</td><td>-1.704 ± 0.012</td></tr><tr><td>t<sub>E</sub> (days)</td><td>185.023 ± 3.914</td><td>155.299 ± 5.269</td></tr><tr><td>s</td><td>1.014 ± 0.001</td><td>1.014 ± 0.001</td></tr><tr><td>q (10<sup>-4</sup>)</td><td>0.553 ± 0.018</td><td>0.750 ± 0.034</td></tr><tr><td>α (rad)</td><td>-0.487 ± 0.007</td><td>0.408 ± 0.017</td></tr><tr><td>ρ<sub>1</sub> (10<sup>-4</sup>)</td><td>1.813 ± 0.060</td><td>2.234 ± 0.074</td></tr><tr><td>s<sub>s</sub></td><td>0.024 ± 0.002</td><td>0.029 ± 0.003</td></tr><tr><td>ψ (rad)</td><td>-0.362 ± 0.089</td><td>3.511 ± 0.059</td></tr><tr><td>qF, I</td><td>0.153 ± 0.013</td><td>0.154 ± 0.021</td></tr><tr><td>π<sub>E,N</sub></td><td>-0.393 ± 0.063</td><td>0.140 ± 0.146</td></tr><tr><td>π<sub>E,E</sub></td><td>-0.090 ± 0.018</td><td>-0.021 ± 0.017</td></tr><tr><td>α<sub>s</sub> (deg)</td><td>107.198 ± 5.614</td><td>84.458 ± 4.336</td></tr><tr><td>δ<sub>s</sub> (deg)</td><td>40.817 ± 12.182</td><td>-29.558 ± 8.312</td></tr><tr><td>f<sub>S</sub></td><td>0.012 ± 0.001</td><td>0.015 ± 0.001</td></tr><tr><td>f<sub>B</sub></td><td>0.784 ± 0.001</td><td>0.781 ± 0.001</td></tr><tr><td>t<sub>*</sub> (days)</td><td>0.034 ± 0.002</td><td>0.035 ± 0.001</td></tr></tbody></table>

Note. $t_* \equiv \rho t_E$ is not an independent quantity.

solutions cannot be distinguished based on goodness-of-fit criteria alone.

We note that the microlensing parameters of the 3L1S and 2L2S parallax solutions are qualitatively similar, but with some quantitative differences. The most notable differences are that the mass ratio of the robustly detected planet and the (first) source radius are smaller by a factor of 0.5–0.7 in the 2L2S solution. This mainly reflects the fact that in high-magnification events, $t_* \equiv \rho t_E$ and $q t_E$ are approximate invariants (Yee et al. 2012), and the value of $t_E$ is substantially higher in the 2L2S solution. Note in particular that the value of $\pi_E$ is qualitatively similar.

In order to assess which, if either, of the 3L1S and 2L2S solutions is preferred, we must first be able to evaluate the implied physical parameters for each.
<figure><img src="image_6.png" /><figcaption>**Figure 4.** CFHT images (14" × 6") of the field near OGLE-2018-BLG-0532 in g (upper left), r (upper right), i (lower left), and combined (lower right). The red circle shows the astrometric position of the microlensing event derived from KMTC difference images near peak, transformed to the CFHT frame. The “baseline object” shown here is nominally composed of the source and blend, but the former contributes negligibly to the total light. The blend is ≲50 mas from the lens, and may be coincident, although this is less likely. The point-spread function (PSF) of the “baseline object” (especially as seen in the i-band image) is not symmetric, but rather has a trefoil-like pattern to the north (up) and west (right) indicating that two fainter stars are partially blended with the PSF. See Sections 4.1.1 and 6.3.</figcaption></figure>

# 4. Color–Magnitude Diagram and Blended Light

The first step toward this goal is to measure $\theta_*$ from the position of the source on the color–magnitude diagram (CMD). This allows one to determine $\theta_E = \theta_*/\rho$ and $\mu_{rel} = \theta_E/t_E$, which are both important inputs into this calculation. In the present case, we also measure the astrometric position and calibrated flux of the blended light, which will also provide important constraints for interpreting the lens system.

## 4.1. Baseline Object

We first note that the various light-curve models described above all obey $F_{base}/f_s > 40$, where $F_{base} \equiv f_s + f_b$. Therefore, to an excellent approximation, the “baseline object” is essentially the same as the blended light.

We conduct our investigation by analyzing relatively high-resolution (FWHM ≃ 0"/45) images taken by the 3.6 m Canada–France–Hawaii Telescope (CFHT) from 2018 June to 2018 August. Because of higher resolution, these images might in principle resolve out stars that are blended with the “baseline object” as seen in OGLE and/or KMT images. Moreover, the higher resolution and smaller pixel size permit a more precise astrometric measurement. See Figure 4.

### 4.1.1. Astrometry

Here our goal is to find the astrometric offset between the microlensing source and the baseline object. If this is consistent with zero, then the blend (essentially the same as the baseline object) could be the lens. And if this offset is small (and possibly zero), it would be evidence that the blend was associated with the event, i.e., either the lens itself, a companion to the lens, or a companion to the source. On the

other hand, if the baseline object were well displaced from the microlensing source, then the lens light could at most comprise a fraction of the baseline light, which would lead to stronger upper limits on the lens flux (see the Appendix).

There are three steps to measuring this offset. First, we must find the position of the microlensed source relative to neighboring field stars from difference images near peak. Second, we must find the position of the CFHT baseline object relative to these same neighboring stars. Third, we must transform the microlensing-template image coordinates to the CFHT image coordinates. We discuss each of these procedures, with a focus on the error estimates of each.

In fact, we have carried out the entire procedure described below twice: once using the source position on the KMT template and once using the source position on the OGLE template. For clarity, we report the KMT analysis first, and then the OGLE analysis. Finally, we investigate the origin of the differences.

**Source Position.** As part of the normal process of image subtraction, each image is astrometrically aligned to the template image, which is then subtracted from it. The difference image then basically consists of an isolated star. We measure this position for 12 highly magnified KMTC images. This sample has standard deviations in the (x, y) (essentially west and north) directions of (0.0435, 0.0320) pixels, which are 400 mas. Hence, the standard errors of the mean of this measurement are (5, 4) mas. For completeness we note that the correlation coefficient is 0.43, although this has no practical impact.

**Baseline Object Position.** In contrast to the microlens difference image, which consists of a single star that is magnified and hence relatively bright, the baseline image is filled with ambient stars. Hence, in principle, the measurement
could be corrupted by ambient light either from gradients from relatively nearby stars, or from very nearby, semi-resolved stars. Fortunately, we find that the baseline object is isolated from relatively bright stars that could cause the first problem. We do find two semi-resolved very nearby faint stars. However, we ignore them here and treat them further below under “systematic errors.”

We construct transformations between pairs of CFHT images, one with excellent (“best”) seeing and the other with good seeing. We measure the scatter in the positions of stars with similar magnitude as the baseline object, finding 0.06 pixels, each pixel being 185 mas. We attribute these as equally due to the two images, i.e., 0.04 pixels each. Using the five best images (and conservatively using the scatter of these measurements, rather than the slightly smaller mean scatter of similar stars mentioned above), we derive standard errors of the mean (6,4) mas in the same (x, y) directions.

**Transformation.** The transformation from the KMT to CFHT images uses the FITSH GRMATCH and GRTRANS routines (Pál 2012) on stars within 1' that are substantially brighter than the baseline object. We find a scatter of these stars in the transformation of about 8 mas. While roughly 200 stars were used to derive the transformation, we conservatively consider only the 50 closest to the baseline object. Then the error in the transformation is ~8/√50 ~ 1 mas. This error is negligible compared to the other errors, and so we ignore it. Hence, finally we derive

$$ \Delta\theta_{\text{KMT}}(E, N) = \theta_{\text{base}} - \theta_{\text{source}} = (-43, +20) \pm (6, 4)\text{mas.} \quad (4) $$

**Systematic Errors.** Equation (4) seems to imply a highly significant (Δχ² ≳ 50) difference in the positions of the baseline object and source, indicating that the baseline object is displaced to the west and north of the source. However, one can see from Figure 4 that there are two semi-resolved stars, which are not registered by the DoPhot (Schechter et al. 1993; Alonso-García et al. 2012) astrometry of the CFHT image, that lie to the west and north. These are certainly biasing the tabulated position of the source in these directions, particularly in the CFHT *i* image. However, the standard configuration of DoPhot does not allow a reliable evaluation of the magnitude of these errors. We discuss their qualitative impact below.

**Comparison to OGLE.** We repeat the entire procedure beginning with the OGLE-based source position and transforming this to the CFHT image using the OGLE template and find

$$ \Delta\theta_{\text{OGLE}}(E, N) = (-61, +31) \pm (5, 5)\text{mas.} \quad (5) $$

While the difference between Equations (4) and (5) is relatively modest relative to the errors, it does have implications for the interpretation of the event. Therefore we investigate it further.

We first note that the OGLE template is from 2010 while the KMT template is from 2018. Because the magnified source is only a factor of η ≃ 1.5 times brighter than the baseline object, the apparent position of the “difference source” is displaced by +μ<sub>base</sub>Δt/η relative to the true position of the source, where μ<sub>base</sub> is the proper motion of the baseline object relative to the reference frame and Δt = 8 yr. Hence, to account for the difference between the OGLE- and KMT-based offsets in the east direction (the only one for which there is even a marginally significant discrepancy) would require a proper motion of μ<sub>base,E</sub> = +3.4 ± 1.5 mas yr⁻¹. This is a very plausible number. Hence, even the marginal discrepancy between the

OGLE and KMT values is entirely consistent with the KMT value being correct. Because the KMT reference image is from the same year as the event, the KMT-based procedure does not suffer from the added uncertainty of the unknown proper motion of the baseline object. Thus, we adopt Equation (4).

We conclude that |Δθ| < 50 mas with the balance of evidence indicating that the baseline object is at least somewhat displaced from the lens. However, due to systematics, the displacement could be substantially smaller than 50 mas, and zero is not ruled out. Hence, the blend could be the lens, but more likely a large fraction of this light is due to another star.

### 4.1.2. Photometry

We calibrate the CFHT instrumental (g, i) photometry (Zang et al. 2018) to calibrated OGLE-III (V, I) photometry (Szymański et al. 2011) to derive $I = i - 0.328 - 0.054(g - i)$ and $(V - I) = 0.326 + 0.716(g - i)$. Combining these with the instrumental magnitude measurements $(g, i)_{\text{base}} = (19.992, 18.596) \pm (0.014, 0.023)$ yields

$$ [(V - I), I]_{\text{base}} = (1.362, 18.596) \pm (0.019, 0.024). \quad (6) $$

### 4.1.3. Role of Surface-brightness Fluctuations

The measurements underlying Equation (6) are derived from point-spread function (PSF) photometry, which treats the light distribution as being the sum of a set of detected stars (convolved with the PSF) plus a uniform background. In fact, in addition to the detected stars, the field also contains stars that are “undetected,” either because they are “unresolved” (too many per PSF to be separately registered) or simply below the detection threshold. In the limit where these “undetected” stars are so numerous that they form an essentially uniform background (on scales of the PSF), then there is no effect. That is, the program measures the mean surface brightness from this “background,” together with the true background due to ambient sky light, and it subtracts this surface brightness from the image before fitting for individual fluxes from the detected stars.

In fact, however, the Poisson distribution of these “undetected” stars will give rise to surface brightness fluctuations (SBF) after they are convolved with the PSF. As first pointed out by Park et al. (2004), this can lead to a “hole” in this distribution relative to neighboring average areas where the mean surface brightness is evaluated, and so to spurious “negative blending.” The same effect can occur when a detected star is present, in which case it will cause the flux of the star to be underestimated. If, as in the present case, the detected star is relatively blue, then the “color” of the “hole” will be redder than the detected star because the fainter, undetected, stars that make up the background are generally redder than the detected star. Hence, the effect of such a “hole” will be to make the detected star appear both bluer and fainter than it actually is. Of course, the opposite case is also possible, i.e., that a relative excess of stars would make the detected star appear brighter and redder than it actually is.

We evaluate this effect quantitatively below and find that it is small in the present case. However, because the effect could be important in other cases, we document the method that we employ.

We model the background by a Holtzman et al. (1998) luminosity function (HLF), which is based on Hubble Space Telescope (HST) images of Baade’s Window (BW). We first
multiply the HLF by 1.57 to account for the higher surface density of bulge stars at this location relative to BW (Nataf et al. 2013; D. Nataf 2019, private communication). Next, we account for the extinction $[E(V - I), A_I] = (0.80, 0.97)$ and the mean bulge distance $D_{\text{bulge}} = 8.0$ kpc toward this field (Nataf et al. 2013). We assume that stars fainter than $I > 20.0$ are detected, so that the “undetected background” is of fainter stars. We create 10,000 realizations of the HLF, each time placing each star at a random position. We recover the apparent flux of the baseline object by convolving the PSF at its location with all the “undetected” stars and then subtracting the background. As foreshadowed above, we find that the error induced by these “undetected” stars is small, $\delta I = \pm 0.09$.

In order to assess the effect on the color, we adopt for these background stars, which are essentially all on the main sequence, $(V - I) = (M_I - 2.9)/2.37$. We then find an even smaller effect, $\delta(V - I) = \pm 0.03$. As a general trend, when the true baseline object is redder than measured, it is also brighter than measured, since both effects are caused by a “hole.”

We note that while a “hole” can only comprise (the absence of) stars below the detection limit, an “excess” could in principle be due, in addition, to stars that are above this limit but are just too close to the baseline object to be separately resolved. However, this effect is already covered under the idea that the baseline object could be due to one or more stars, either related or unrelated to the event.

In the present case, the main interest is in possible systematic errors in the color measurement that might account for why the baseline object appears so blue. That is, if we subtract out the field reddening, the baseline object would have a color $(V - I)_0 = 0.56$. This exercise tells us that, in addition to the effect of the measurement error in Equation (6) the true color might be redder by 0.03 mag at $1\sigma$. For reference, we note that at $2\sigma$, it could be 0.04 mag redder.

### 4.1.4. Historical Case of MOA-2013-BLG-220

There are very few cases for which predictions about the lens that take account of SBF can be tested in practice. In fact, there are only four planetary events for which the lens has been separately resolved from the source, OGLE-2005-BLG-071 (Bennett et al. 2020), OGLE-2005-BLG-169 (Batista et al. 2015; Bennett et al. 2015), OGLE-2012-BLG-0950 (Bhattacharya et al. 2019), and MOA-2023-BLG-220 (Vandorou et al. 2020). The one case among these with such predictions, MOA-2013-BLG-220 (Yee et al. 2014; Vandorou et al. 2020), provides a cautionary tale.

In their Section 4.2, Yee et al. (2014) adopted a limit on (the logarithm of) the lens-source flux ratio of $(I_L - I_S) > 1$, which allowed for lenses that were substantially brighter than their measurement of the blended light $(I_B - I_S) = 1.9$. This limit corresponds to $I_L < 20.20$. Yee et al. (2014) did not report their reasoning for this limit. However, by applying the method described just above to account for “holes” in the mottled bulge background, we find that this represents the 90% confidence limit, which could be regarded as moderately conservative.

Nevertheless, when Vandorou et al. (2020) separately imaged the source and lens, they found that the lens was relatively bright: $K_L = 17.92 \pm 0.05$. To relate this measurement in K

with the predicted limit in I, we write

$$I_L = K_L + (I - K)_{0,L} + E(I - K)_L. \quad (7)$$

In the solution presented by Vandorou et al. (2020) based on their $K_L$ measurement and their analysis of the microlensing light curve, they found that the lens lies in or near the bulge. Hence, $E(I - K)_L = E(I - K)_S = 1.005 \pm 0.05$. Vandorou et al. (2020) find that their measurement of $K_L$ combined with the light-curve analysis yields a lens mass $M_L = 0.89 M_\odot$, for which we adopt $(I - K)_0 = 0.85 \pm 0.04$, which is characteristic of a G5 star (Bessell & Brett 1988). Hence, Equation (7) yields $I_L = 19.77 \pm 0.08$, which is considerably brighter than the Yee et al. (2014) “limit” $I_L > 20.20$. A full analysis of this case to resolve this contradiction would be well beyond the scope of the present work. However, it does make clear the need for additional work on understanding blended-light limits (Gould et al. 2020).

### 4.2. CMD

We will ultimately wish to place the source color and magnitude $(V - I, I)_S$ on an instrumental CMD in order to find its offset from the centroid of the red clump and so infer its angular radius $\theta_*$ (Yoo et al. 2004). This in turn will yield the angular Einstein radius via the relation $\theta_E = \theta_*/\rho$, where $\rho$ is the normalized source radius, which is a parameter of the microlens modeling.

While $I_S$ can only be determined by conducting detailed modeling of the microlensing event, $(V - I)_S$ is independent of the model, and in fact can often be determined by the regression of the V on I fluxes without any model at all.<fnref n="19" /> In the present case, we will be considering models with somewhat different values of $I_S$. The purpose of this section is therefore to derive a scaled relation for $\theta_*$ that is valid for all of these. This is feasible because for fixed source color (and hence fixed inferred source surface brightness), $\theta_* \propto f_S^{1/2}$.

Hence, we will simply adopt a fiducial source flux $f_{s,fid} = 0.020$, which corresponds to $I_{s,fid} = 22.23$ on an $I = 18$ flux scale. This value corresponds to the source flux in the “standard” 2L1S model in Table 1, but the choice of this particular value is simply a matter of convenience. Then, other than (temporarily) treating the error in $I_S$ as zero, the analysis proceeds in the usual way.

We begin by constructing the OGLE-IV CMD (Figure 5) and measuring the position of the red clump centroid $[(V - I), I]_{cl} = (1.82, 15.28) \pm (0.04, 0.08)$. Next we measure the source color $(V - I)_S = 2.48 \pm 0.10$. In practice we do this from the models rather than regression, but we confirm that the measurement and error are the same for the different models (as expected). Together with the adopted fiducial $I_{s,fid} = 22.23$, these imply an offset from the clump of $\Delta[(V - I), I] = (0.66, 6.95) \pm (0.11, 0.08)$.

We repeat a similar procedure with KMTC43 data and obtain $\Delta[(V - I), I]_{standard,KMT} = (0.69, 6.86) \pm (0.10, 0.08)$. These results are consistent at the $1\sigma$ level. Hence, we average them and obtain $\Delta[(V - I), I] = (0.67, 6.91) \pm (0.07, 0.08)$.

We adopt $[(V - I), I]_{0,\text{clump}} = (1.06, 14.39)$ from Bensby et al. (2013) and Nataf et al. (2013), which yields $[(V - I), I]_{0,S} = (1.73, 21.30) \pm (0.07, 0.08)$, where we are ignoring for the moment any errors in the intrinsic position of

<footnote n="19">¹⁹ This is not strictly true for 1L2S or 2L2S events in which the source stars have different colors. However, we will show that this has no significant impact in the present case.</footnote>
<figure><img src="image_7.png" /><figcaption>**Figure 5.** OGLE-IV color-magnitude diagram (CMD) for OGLE-2018-BLG-0532. The source star, blended light, and clump centroid are shown as black, green, and red circles, respectively. These results are combined with a similar CMD derived from KMTC43 data to measure the offset of the source star from the clump and so to derive the angular source radius, $\theta_*$.</figcaption></figure>

The astrometry and photometry of the blended light are discussed in Section 4.1, based on the CFHT images shown in Figure 4. The clump. Then converting from $V/I$ to $V/K$ using the color-color relations of Bessell & Brett (1988), $(V-K)_{0,s} = 3.55 \pm 0.07$ and $A_K = 0.134$, and then applying the “dwarf and sub-giant” color/surface-brightness relations of Kervella et al. (2004), we obtain

$$ \theta_{*,\mathrm{fid}} = 0.388 \pm 0.019 \ \mu\mathrm{as}. \qquad (8) $$

As noted above, the error bar in Equation (8) reflects only the errors in centroiding the clump and in measuring the source color. In particular, it does not include the measurement error of $f_s$, nor does it include possible systematic errors in the color/surface-brightness relations or the intrinsic position of the clump. To take account of all of these errors, we add 5% in quadrature to the $\approx 5\%$ error reported in Equation (8), i.e., 7% in total. Because we discuss other solutions with somewhat different source fluxes $f_s$, we note that these also have somewhat different $\theta_*$. Including this additional 5%, the general formula for $\theta_*$ becomes

$$ \theta_* = \sqrt{\frac{f_s}{0.020}} (0.388 \pm 0.027 \ \mu\mathrm{as}). \qquad (9) $$

We note that the best-fit values for $\theta_E = \theta_*/\rho$ for the 12 models considered in Tables 1–4 are therefore $\theta_E = (1.23, 1.35, 1.34)\mathrm{mas}$, $\theta_E = (1.30, 1.14, 1.20, 1.23)\mathrm{mas}$, $\theta_E = (1.29, 1.43, 1.48)\mathrm{mas}$, and $\theta_E = (1.66, 1.50)\mathrm{mas}$, respectively.

## 5. Physical Parameters

We now determine the physical parameters, such as the mass and distance of the lens system,

$$ M = \frac{\theta_E}{\kappa\pi_E}, \quad D_L = \frac{au}{\theta_E\pi_E + \pi_s}, \qquad (10) $$

for the 10 solutions presented in Tables 1–4. We directly evaluate these from the output of the MCMC. For example, the projected physical separation of a given MCMC realization of the 2L1S solutions is

$$ \begin{aligned} a_{\perp} &\equiv sD_L\theta_E = \frac{s \cdot au}{\pi_E + \pi_s/\theta_E} \\ &= \frac{s \cdot au}{(\pi_{E,N}^2 + \pi_{E,E}^2)^{1/2} + \pi_s\rho/[\theta_{*,\mathrm{fid}}(f_s/0.020)^{1/2}(1+\delta)]} \end{aligned} \qquad (11) $$
<table><caption>**Table 5**
Physical Parameters for 2L1S Models</caption><thead><tr><th>Quantity</th><th>u<sub>0</sub> > 0</th><th>u<sub>0</sub> < 0</th></tr></thead><tbody><tr><td>M<sub>lens</sub> [M<sub>☉</sub>]</td><td>0.195<sup>+0.022</sup><sub>-0.020</sub></td><td>0.203<sup>+0.022</sup><sub>-0.020</sub></td></tr><tr><td>M<sub>planet</sub> [M<sub>☉</sub>]</td><td>6.228<sup>+0.905</sup><sub>-0.881</sub></td><td>6.554<sup>+0.905</sup><sub>-0.808</sub></td></tr><tr><td>a<sub>⊥</sub> [au]</td><td>1.064<sup>+0.119</sup><sub>-0.115</sub></td><td>1.103<sup>+0.118</sup><sub>-0.107</sub></td></tr><tr><td>a<sub>⊥</sub>/a<sub>snow</sub></td><td>2.015<sup>+0.051</sup><sub>-0.052</sub></td><td>2.011<sup>+0.051</sup><sub>-0.051</sub></td></tr><tr><td>D<sub>L</sub> [kpc]</td><td>0.773<sup>+0.102</sup><sub>-0.099</sub></td><td>0.804<sup>+0.100</sup><sub>-0.092</sub></td></tr><tr><td>μ<sub>geo</sub> [mas yr<sup>-1</sup>]</td><td>3.587<sup>+0.198</sup><sub>-0.170</sub></td><td>3.594<sup>+0.209</sup><sub>-0.170</sub></td></tr><tr><td>μ<sub>hel,N</sub> [mas yr<sup>-1</sup>]</td><td>-2.880<sup>+0.154</sup><sub>-0.177</sub></td><td>4.230<sup>+0.273</sup><sub>-0.221</sub></td></tr><tr><td>μ<sub>hel,E</sub> [mas yr<sup>-1</sup>]</td><td>2.043<sup>+0.412</sup><sub>-0.341</sub></td><td>2.196<sup>+0.360</sup><sub>-0.317</sub></td></tr><tr><td>˜v<sub>lsr,l</sub> [km s<sup>-1</sup>]</td><td>5.872<sup>+1.240</sup><sub>-1.814</sub></td><td>32.279<sup>+1.540</sup><sub>-1.386</sub></td></tr><tr><td>˜v<sub>lsr,b</sub> [km s<sup>-1</sup>]</td><td>-6.049<sup>+0.663</sup><sub>-0.705</sub></td><td>7.809<sup>+1.148</sup><sub>-0.984</sub></td></tr></tbody></table>

Note. The error bars for $a_\perp/a_{\text{snow}}$ are calculated under the definition $a_{\text{snow}} \equiv 2.7 \text{ au}(M/M_\odot)$.

Because $\pi_E \theta_E = \pi_{\text{rel}} \gg \pi_s$ (i.e., the lens is much closer to the Sun than to the source), we treat the source distance as a constant $D_S = \text{au}/\pi_s = 8.0$ kpc. The “$\delta$” in the denominator represents the 7% error in $\theta_*$.

It is implemented by integrating over a Gaussian for each MCMC realization. The remaining variables ($s$, $\pi_{E,N}$, $\pi_{E,E}$, $\rho$, $f_s$) come directly from the MCMC.
Tables 5–8, show the resulting physical parameters for the 2L1S, 3L1S, 2L2S-static, and 2L2S-orbital-motion solutions, respectively. To scale the projected separation to the snow line, we adopt $a_{\text{snow}} = 2.7(M/M_\odot)$au. The heliocentric and local standard of rest (LSR) quantities are derived by

$$ \mu_{\text{hel}} = \mu_{\text{geo}} + \frac{\pi_{\text{rel}}}{\text{au}} \mathbf{v}_{\perp, \oplus}: \quad \tilde{\nu}_{\text{hel}} \equiv \frac{\text{au}}{\pi_E} \frac{\mu_{\text{hel}}}{\mu_{\text{hel}}^2} \quad (12) $$

and

$$ \tilde{\nu}_{\text{lsr}} = \tilde{\nu}_{\text{hel}} + \mathbf{v}_{\perp, \odot}. \quad (13) $$

Here $\mu_{\text{geo}} = \mu_{\text{rel}}$ is the geocentric relative proper motion, $\mu_{\text{hel}}$ is the heliocentric proper motion, $\mathbf{v}_{\perp, \oplus}(N, E) = (2.7, 10.2) \text{ km s}^{-1}$ is Earth’s velocity projected onto the plane of the sky at $t_0$, $\mathbf{v}_{\perp, \odot}(l, b) = (12, 7) \text{ km s}^{-1}$ is the Sun’s peculiar velocity projected onto the plane of the sky, and $\tilde{\nu}_{\text{hel}}$ and $\tilde{\nu}_{\text{lsr}}$ are the “projected velocities” (Gould 1992) in the heliocentric and LSR frames, respectively.

The physical properties of the host-Neptune system are quantitatively similar for the six single-source (2L1S + 3L1S) solutions shown in Tables 5 and 6 except for the direction of $\mu_{\text{hel}}$. This follows from the similarity of the underlying microlensing parameters in Tables 1 and 2. It is expected that the four pairs of [(u₀ < 0), (u₀ > 0)] will be similar to each other (apart from the direction of $\pi_E$). The fact that the properties of the Neptune do not change much when the light curve is fit for an additional planet conforms to the prediction of Zhu et al. (2014), who found from simulations that failure to take account of “unrecognized” second planets generally does not significantly affect the parameter estimates of the first planet. However, Zhu et al. (2014) did not simultaneously fit for parallax, so the fact that the estimate of $\pi_E$ is robust against the presence of a planet with a $\Delta\chi^2 \sim 80$ signature could not necessarily have been anticipated.

The host mass is estimated at 0.20–0.25 $M_\odot$, at a distance $D_L \sim 1$ kpc. At this mass, its absolute magnitude should be

<table><caption>**Table 6**
Physical Parameters for 3L1S Models</caption><thead><tr><th rowspan="2">Quantity</th><th colspan="2">Close</th><th colspan="2">Wide</th></tr><tr><th>u<sub>0</sub> > 0</th><th>u<sub>0</sub> < 0</th><th>u<sub>0</sub> > 0</th><th>u<sub>0</sub> < 0</th></tr></thead><tbody><tr><td>M<sub>lens</sub> [M<sub>☉</sub>]</td><td>0.238<sup>+0.056</sup><sub>-0.029</sub></td><td>0.254<sup>+0.110</sup><sub>-0.043</sub></td><td>0.242<sup>+0.041</sup><sub>-0.026</sub></td><td>0.254<sup>+0.076</sup><sub>-0.034</sub></td></tr><tr><td>M<sub>planet,1</sub> [M<sub>☉</sub>]</td><td>8.041<sup>+2.877</sup><sub>-1.274</sub></td><td>9.505<sup>+5.645</sup><sub>-2.224</sub></td><td>7.697<sup>+1.894</sup><sub>-1.108</sub></td><td>8.309<sup>+3.995</sup><sub>-1.501</sub></td></tr><tr><td>M<sub>planet,2</sub> [M<sub>☉</sub>]</td><td>0.679<sup>+0.309</sup><sub>-0.188</sub></td><td>0.613<sup>+0.346</sup><sub>-0.197</sub></td><td>0.746<sup>+0.341</sup><sub>-0.239</sub></td><td>0.743<sup>+0.440</sup><sub>-0.241</sub></td></tr><tr><td>a<sub>⊥</sub> [au]</td><td>1.363<sup>+0.309</sup><sub>-0.156</sub></td><td>1.484<sup>+0.548</sup><sub>-0.248</sub></td><td>1.364<sup>+0.225</sup><sub>-0.136</sub></td><td>1.434<sup>+0.402</sup><sub>-0.181</sub></td></tr><tr><td>a<sub>2⊥</sub> [au]</td><td>0.506<sup>+0.146</sup><sub>-0.087</sub></td><td>0.599<sup>+0.244</sup><sub>-0.126</sub></td><td>3.819<sup>+0.817</sup><sub>-0.583</sub></td><td>3.963<sup>+1.080</sup><sub>-0.672</sub></td></tr><tr><td>a<sub>⊥</sub>/a<sub>snow</sub></td><td>2.113<sup>+0.080</sup><sub>-0.084</sub></td><td>2.130<sup>+0.086</sup><sub>-0.102</sub></td><td>2.083<sup>+0.078</sup><sub>-0.079</sub></td><td>2.075<sup>+0.088</sup><sub>-0.102</sub></td></tr><tr><td>a<sub>2⊥</sub>/a<sub>snow</sub></td><td>0.784<sup>+0.104</sup><sub>-0.102</sub></td><td>0.851<sup>+0.099</sup><sub>-0.104</sub></td><td>5.724<sup>+0.720</sup><sub>-0.560</sub></td><td>5.536<sup>+0.740</sup><sub>-0.621</sub></td></tr><tr><td>D<sub>L</sub> [kpc]</td><td>1.091<sup>+0.307</sup><sub>-0.153</sub></td><td>1.239<sup>+0.529</sup><sub>-0.255</sub></td><td>1.077<sup>+0.233</sup><sub>-0.133</sub></td><td>1.150<sup>+0.398</sup><sub>-0.184</sub></td></tr><tr><td>μ<sub>geo</sub> [mas yr<sup>-1</sup>]</td><td>3.263<sup>+0.216</sup><sub>-0.191</sub></td><td>3.314<sup>+0.217</sup><sub>-0.190</sub></td><td>3.162<sup>+0.203</sup><sub>-0.180</sub></td><td>3.178<sup>+0.205</sup><sub>-0.170</sub></td></tr><tr><td>μ<sub>hel,N</sub> [mas yr<sup>-1</sup>]</td><td>-2.785<sup>+0.181</sup><sub>-0.203</sub></td><td>3.670<sup>+0.274</sup><sub>-0.226</sub></td><td>-2.686<sup>+0.186</sup><sub>-0.186</sub></td><td>3.569<sup>+0.244</sup><sub>-0.234</sub></td></tr><tr><td>μ<sub>hel,E</sub> [mas yr<sup>-1</sup>]</td><td>1.220<sup>+0.342</sup><sub>-0.492</sub></td><td>1.228<sup>+0.474</sup><sub>-0.658</sub></td><td>1.286<sup>+0.298</sup><sub>-0.376</sub></td><td>1.375<sup>+0.375</sup><sub>-0.570</sub></td></tr><tr><td>˜v<sub>lsr,l</sub> [km s<sup>-1</sup>]</td><td>1.167<sup>+2.800</sup><sub>-6.400</sub></td><td>38.570<sup>+10.095</sup><sub>-4.299</sub></td><td>1.957<sup>+2.400</sup><sub>-2.283</sub></td><td>36.000<sup>+6.500</sup><sub>-2.838</sub></td></tr><tr><td>˜v<sub>lsr,b</sub> [km s<sup>-1</sup>]</td><td>-7.633<sup>+1.033</sup><sub>-2.033</sub></td><td>12.255<sup>+8.475</sup><sub>-3.020</sub></td><td>-7.486<sup>+0.914</sup><sub>-1.500</sub></td><td>10.375<sup>+5.250</sup><sub>-2.000</sub></td></tr></tbody></table>

Note. The error bars for $a_\perp/a_{\text{snow}}$ are calculated under the definition $a_{\text{snow}} \equiv 2.7 \text{ au}(M/M_\odot)$.

$M_I \gtrsim 9$, and so allowing for an extinction of $A_{l,l} \sim 0.5$ mag to this distance, the predicted lens flux is $I_l \gtrsim 19.5$. This is below the limit from the blended light from Equation (6), $I_l > 18.6$. Hence, the blended light appears to be dominated by some other star, which would be consistent with the “balance of evidence” in Section 4.1.1 that the baseline object is offset from the lens. Nevertheless, the fact that the predicted lens flux is of the same order as the blend flux does raise some subtle issues, the discussion of which we defer to Section 6.

The first (robustly detected) planet has a mass of 6–10 $M_⊕$ across the six solutions, which is near the mass $M_p \sim 10 M_⊕$ generally thought to be required for rapid growth by gas accretion. We also leave the implications of this mass estimate to the discussion. It lies projected at about twice the snow line distance.

The second (putative) planet has a mass midway between that of Jupiter and Saturn. It lies projected either at $\sim 0.8$ or $\sim 5.6$ times the snow line, depending on whether the close or wide solution is correct. The close solution is very slightly preferred by $\chi^2$. From general statistics of such gas giants, they are more likely to be found outside than inside the snow line, but for this individual case we cannot distinguish.

On the other hand, the physical parameters of the host and planet in the 2L2S solution with orbital motion (Table 8) differ substantially from the single-source solutions. This is particularly notable for the $u_0 < 0$ solution in Table 8. The fact that such a marginally detected second source can have such an impact is potentially of wide significance. We defer discussion of its relevance to the present case of OGLE-2018-BLG-0532 to Section 6.

## 6. Discussion

The OGLE-2018-BLG-0532L system is interesting for two major reasons. First, it shows significant evidence for a second planet, which would potentially make it only the third two-planet system discovered by microlensing. Second, its securely detected planet, OGLE-2018-BLG-0532Lb lies near the possible “pileup” of planets identified by Jung et al. (2019)
that is centered just below $q = 10^{-4}$. In addition, the blended light in this event raises some puzzles that deserve further discussion.

## 6.1. Cold Neptune

OGLE-2018-BLG-0532Lb is the second cold Neptune to be discovered since the recent analysis by Jung et al. (2019) of the mass-ratio function below $q < 3 \times 10^{-4}$. They had argued for a sharp break in the mass-ratio function at $q_{brk} \sim 0.56 \times 10^{-4}$ and/or a “pile up” of cold Neptunes just above the putative break point. Both of these planets (the other being OGLE-2015-BLG-1670Lb; Ranc et al. 2019), have mass-ratio estimates whose error bars straddle $q = 1.0 \times 10^{-4}$. This places them near the upper end of the putative “pileup.” While it would be premature to revisit the form of the mass-ratio function at this point, we note that the pace of discovery of planets in the regime $q \lesssim 1 \times 10^{-4}$ is accelerating. The discovery years for the nine planets in this regime are (2005, 2005, 2007, 2009, 2013, 2015, 2016, 2017, and 2018). That is, the discovery pace has roughly doubled starting in 2015. This can at least be partly attributed to the inauguration of KMTNet in that year in the sense that KMTNet played a critical role in two of the four discoveries since 2015, and an auxiliary role in the other two. This tends to confirm that the apparent doubling of the discovery rate for low-$q$ planets is real and not the result of a statistical fluctuation. Hence, it is likely that within a few years we will gain a significantly better picture of the mass-ratio function at the low end.

## 6.2. Two-planet System?

While the “first planet” (the Neptune) in this system is readily apparent from the light curve, even without detailed modeling, visual inspection of the light curve does not give even the slightest hint of the “second planet” (the Jupiter). This contrasts sharply with the first two cases of microlensed two-planet systems, OGLE-2006-BLG-109Lb,c (Gaudi et al. 2008; Bennett et al. 2010) and OGLE-2012-BLG-0026Lb,c (Han et al. 2013), wherein the light curves basically “factored” into perturbations due to each planet separately (Han et al. 2001). The present case is more similar to the simulations of two-planet systems in the KMTNet-like data streams that were investigated by Zhu et al. (2014), in which the “second-planet” signature was often a weak featureless deviation. Nevertheless, while these second planets were often not obvious to the eye, they did clearly stand out in residuals to the fits to the first planet. This is not true in the present case.

Rather, as discussed in the Appendix, 3L1S solutions were only investigated because some intermediate results of the analysis seemed inconsistent with the upper limits on lens flux. This investigation consumed a lot of resources (both human and computer), so such investigations are not generally undertaken in the analysis of planetary microlensing events in the absence of discernible evidence for a third body. As we noted in Section 1, however, Shin et al. (2015) did search eight archival planetary events for evidence of third bodies and found $\Delta\chi^2 \sim 50-142$ in three of these cases, which are similar to the $\Delta\chi^2 = 81$ improvement for OGLE-2018-BLG-0532.

We argue that this planetary signal is plausibly real, but we consider that the evidence for the planet is not sufficient to definitively claim its detection. In this regard, there are two distinct questions that must be examined. First, is the

$\Delta\chi^2 = 81$ difference between 3L1S and 2L1S secure enough so that this improvement cannot plausibly be attributed to systematics? Second, is evidence in favor of 3L1S over 2L2S strong enough to confidently claim that the additional “body” relative to 2L1S is an additional lens (3L1S) rather than an additional source (2L2S)? We will answer “yes” to the first question and “no” to the second.

Regarding the preference of 3L1S to 2L1S, there are three factors that all favor accepting the 3L1S solution. First, the $\Delta\chi^2 = 81$ improvement is relatively high. It is true that various authors have advocated much higher thresholds for “detectability” of planets, but the reasons for these higher thresholds must be clearly understood. For example, Gould et al. (2010) advocated for a range of possible thresholds centered on $\Delta\chi^2 = 500$. However, their argument was not that this was the minimum required to detect, or even to characterize a planet. Rather, it was the minimum $\Delta\chi^2$ using re-reduced data needed to guarantee that a planet would be sufficiently apparent in pipeline data to trigger a deeper investigation. Without this “guarantee” that all planets meeting some objective criterion have been found in a given sample, one cannot carry out a statistical analysis of the sample. But the problem of finding a second planet in a system already known to contain one is quite different from that of finding the first planet. First, the data are already re-reduced. Second, the sample (events with planets) that is being probed for second planets is much smaller than the sample (all events) that was probed for first planets. Hence, the probability of some unknown systematics introducing a “planetary signal” is likewise reduced. Finally, in the present case, we are not trying to create a statistical sample, but only to assess the evidence that the planet is real. Therefore there is no need to set a threshold that is “high enough” that the great majority of planets meeting it would be recognized “by eye.”

Second, the cumulative distribution $\Delta\chi^2(t) = \chi^2_{2\text{L1S}}(t) - \chi^2_{3\text{L1S}}(t)$ is quite consistent with a real planet and would require a remarkable set of coincidences if it were due to systematics. As shown in Figure 6, all eight observatory/field combinations contribute positively to the total $\Delta\chi^2$. Moreover, as would be expected, most of the contribution is from the regions of the caustic structure, where a weak shear (or pseudo-shear) would generate the most pronounced effects. See Figure 7. And again, within this critical region, all observatory/field combinations contribute positively. Moreover, the whole light curve does weakly contribute as well, also as one would expect.<fnref n="20" />

Third, a weak signal from a second planet in a high-magnification event should not be regarded as unexpected. As discussed in Section 1, high-magnification events are simultaneously sensitive to all planets in the system, provided that they are close enough to the Einstein ring to generate a sizable

<footnote n="20">²⁰ In principle, all eight observatory/field combinations could also contribute positively if there were a common physical cause due to unmodeled physics of the microlensing event or coherent variability of the source or blended light. As we discussed in Sections 3.1 and 3.5, orbital motion does not give rise to significant signals during the event. The blend is an upper main-sequence or turnoff star and so is not expected to be significantly variable. Nevertheless, we search for such putative variability in binned residuals to the 2L1S fit, but find no coherent deviations except near the peak. Source variability at the required level could not be detected except when highly magnified (i.e., near the peak) because it is too faint. However, the source is a late K or early M dwarf and so is also not expected to vary on few day timescales. Moreover, a minority (but still significant part) of the 3L1S signature comes from the long-term behavior of the event, which would require that the source vary on two different timescales in a “cooperating” fashion. We conclude that all such explanations by real physical effects are unlikely.</footnote>
<figure><img src="image_8.png" /><figcaption>**Figure 6.** Cumulative $\Delta\chi^2(t) = \chi^2_{2\text{LIS},(u_0>0)}(t) - \chi^2_{3\text{LIS},(u_0>0,\text{close})}(t)$ function for the addition of a second planet in OGLE-2018-BLG-0532. All eight observatory/ field combinations contribute positively, and most of the signal comes from the region of the anomaly. Both of these facts contribute to the confidence that the putative second planet is not due to systematics (however, see Section 3.6).</figcaption></figure>

central caustic. At the time of the Shin et al. (2015) study, there were a total of nine “high-magnification” ($u_0 < 0.01$) events that contained published planets. This set consists of six of the eight events analyzed by Shin et al. (2015; i.e., not including OGLE-2005-BLG-071 and MOA-2009-BLG-387, which had $u_0 > 0.01$), plus three additional events: OGLE-2006-BLG-109 (Gaudi et al. 2008; Bennett et al. 2010), OGLE-2012-BLG-0026 (Han et al. 2013), and OGLE-2007-BLG-349 (Gould et al. 2010; Bennett et al. 2016). The first two of these three additional events showed very clear evidence of two planets, while the last had very clear residuals from the 2L1S fit due to a binary companion to the host, i.e., circumbinary planet. Of the six events that did not have discernible systematic residuals, i.e., did not require a third body to achieve a satisfactory fit to the light curve, Shin et al. (2015) found that two had evidence for a third body at a similar or higher $\Delta\chi^2$ as OGLE-2018-BLG-0532. Further, one other had

evidence at a lower level ($\Delta\chi^2 = 50$) that still could be a plausible candidate. While Shin et al. (2015) did not investigate all of these marginal detections in detail, they did mention that three different observatories contributed significantly to the $\Delta\chi^2 = 142$ of their best case. Stated otherwise, only three of the nine cases in this complete sample of high-magnification planetary events showed no significant evidence ($\Delta\chi^2 < 30$) for third bodies.

To further illuminate this issue, we carry out simulations based on OGLE-2018-BLG-0532 to determine how often one would expect “weak-but-detectable” signals versus “strong-obvious” signatures, assuming that a Jovian-mass-ratio planet like the possible second planet in OGLE-2018-BLG-0532 was somewhere in a system basically defined by the robustly detected planet OGLE-2018-BLG-0532Lb. To do so, we create simulated light curves following the procedure of Udalski et al. (2018). That is, we measure the residuals from the 3L1S model
<figure><img src="image_9.png" /><figcaption>Figure 7. Zoom of the cumulative $\Delta\chi^2(t) = \chi^2_{2\text{LIS}}(t) - \chi^2_{3\text{LIS}}(t)$ function shown in Figure 6.</figcaption></figure>

and add these to models with the same $q_2$ but different $s_2$ and $\psi$ compared to the best-fit model. The seven standard parameters ($t_0, u_0, t_E, s_1, q_1, \alpha, \rho$) are kept the same. Then we fit to a standard 2L1S model (i.e., with seven parameters) and find the increase in $\Delta\chi^2$ due to the absence of the second planet in the model. Note that to enhance the computational speed, we both create 3L1S models and fit to 2L1S models with $\pi_E = 0$. Conceivably, this might alter the $\Delta\chi^2$ due to the additional degrees of freedom. However, this seems unlikely on general grounds and in the few cases that we checked, the fractional change in $\Delta\chi^2$ was a few percent, i.e., well below the level of interest for this theoretical study.

Figures 8–10 show the results of these simulations for $\ln s_2 = (1, 0.75, 0.50)$, i.e., $s = (e^1, e^{3/4}, e^{1/2})$, respectively. Note that the first of these is very similar to the best-fit value for the real data, $s_2 = 2.65$. In each case, we examine 10 values of $\psi' = \psi_{\text{best}} - 2\pi n/10$ with $n = 0, 1, \dots, 9$. For each simulation, we show the caustic geometry on the left and $\Delta\chi^2$ on the right.

We see that at $\ln s = 1$, only two of the 2L1S fits look “clearly suspicious,” i.e., $n = 5$ and $n = 6$. The first shows strong systematic residuals in KMTA and OGLE/KMTC data on 8217.xx. The second shows such residuals for OGLE/KMTC data on 8217.xx and KMTA data on 8218.xx. However, none of the 10 cases show residuals that are obviously due to a third body.

By contrast, for $\ln s_2 = 0.75$, there are two cases (again $n = 5$ and $n = 6$) with residuals that clearly look like perturbations due to a third body, although if one were sufficiently suspicious of the data one might be inclined to dismiss them at “just systematics in the data.” Only at $\ln s_2 = 0.5$ are there several cases ($n = 5, n = 6$, and $n = 8$) for which the light curve appears to be a superposition of two planetary perturbations, as predicted by Han et al. (2001) and as actually observed in the cases of OGLE-2006-BLG-109 (Gaudi et al. 2008; Bennett et al. 2010) and OGLE-2012-BLG-0026 (Han et al. 2013). However, for this $\ln s_2 = 0.5$ case, nearly all the examples have very obvious residuals that would almost certainly prompt
<figure><img src="image_10.png" /><figcaption>**Figure 8.** Simulated 3L1S light-curve data (colored points) with best-fit 2L1S models (cyan curves), for separation $ln s_2 = 1$ and various values of the angle $\psi' = \psi - 2\pi n/10$ ($n = 0, 1, ..., 9$) between the two planets, but with the second planet having the same mass ratio as in the best fit. The 3L1S caustic geometries are shown on the left, and the $\Delta\chi^2 = \chi^2(2L1S) - \chi^2(3L1S)$ values are shown on the right. The $n=5$ and $n=6$ cases show noticeable offsets of the data from the models, but in other cases the differences are difficult to discern by eye.</figcaption></figure>

investigation for an additional body or bodies in the system. We expect that for second planets even closer to the Einstein ring, $0 < \ln s_2 < 0.5$, essentially all cases would yield obvious residuals that would have to be investigated. Further, due to the $s \leftrightarrow s^{-1}$ degeneracy, we expect that detectability would be symmetric with respect to the Einstein ring.

If we then consider potential planets that are uniformly distributed in $\ln s$ (i.e., Opik's law), the instances of "obvious signatures" and "significant but non-obvious" signatures are roughly comparable for the geometry and light-curve coverage of OGLE-2018-BLG-0532. This study is meant to be only illustrative. A full investigation of the statistics of such weak signals would require systematic modeling of all high-magnification events with planets. This would be well beyond the scope of the present paper.

Thus, the answer to our first question is that based solely on comparison of the 3L1S to the 2L1S solutions, we would conclude that the apparent second (i.e., Jovian) planetary signal

is not due to systematics. Given that the signature of this planet was invisible to the eye, even in the residuals, we think that it would be prudent to systematically search for third bodies in all high-magnification planetary events.

Moreover, such a systematic search could reveal additional binary companions to planetary systems. To date, there are four microlensing planets in microlens-binary systems: OGLE-2007-BLG-349 (Bennett et al. 2016), OGLE-2013-BLG-0341 (Gould et al. 2014), OGLE-2016-BLG-0613 (Han et al. 2017), and OGLE-2008-BLG-092 (Poleski et al. 2014). The first two were in high-magnification events. For both of the middle two, the binary companion provides by far the dominant signal in the light curve, while for the first it generates very noticeable residuals to the single-planet fit. The last case has a completely different geometry, in which each of the three bodies gives rise to a nearly isolated microlensing event. Detection of additional systems would be of interest in their own right. In addition,
<figure><img src="image_11.png" /><figcaption>**Figure 9.** Simulated 3L1S light-curve data with 2L1S models, similar to Figure 8, except that $\ln s_2 = 0.75$. For $n = 5$ and $n = 6$, the light-curve deviations strongly suggest a third body, and several other panels show obvious deviations.</figcaption></figure>

inclusion of the third (binary-companion) body in the fit could change the parameters of the planet.

We turn now to the second question: can we confidently accept 3L1S over 2L2S? Based on the light-curve analysis alone, the 3L1S solution is better by only $\Delta\chi^2 = 3.8$. Given the level of systematics that is typical of microlensing light curves, this small $\chi^2$ difference would not be enough to clearly prefer (let alone decisively favor) the 3L1S solution.

Another possible method to distinguish between “extra-source” versus “extra-lens” solutions is a color test of the “extra-source” solution (i.e., 2L2S in the present case). For example, Hwang et al. (2018) showed that the color-offset between two of the three sources in the 1L3S solution of OGLE-2015-BLG-1459 clearly confirmed it over the otherwise degenerate 2L2S and 3L1S solutions. In the present case, the two sources of the 2L2S solution are both on the main sequence with the secondary being $\Delta I \simeq 1.9$ mag fainter than the primary. Hence, one would expect that it would also be redder by $\Delta(V-I) \simeq 0.8$. On the other hand, if the 2L2S solution were simply mimicking a light curve generated by a

3L1S (i.e., single-source, hence achromatic) event, then we would expect $\Delta(V-I) \to 0$. Thus, this provides a potentially clear test.

We therefore include all V data into the fit and model these by an additional parameter $q_{F,V}$ and evaluate $\Delta(V-I) = 2.5 \log(q_{F,I}/q_{F,V})$. Unfortunately, due to the redness of the sources as well as the relatively bright blue blend and the relative paucity of V-band data (1 for every 10 I-band), this test does not yield decisive results. We find that the best fit is $\Delta(V-I) \sim 0.33$, which is midway between the expected values for the 2L2S and 3L1S. Moreover, the predicted 3L1S value (zero) is disfavored by only $\Delta\chi^2 = 3.5$ (for 1 dof). Hence, this test provides no clear evidence for either solution.

Another point that prevents definitely claiming the detection of the second planet is that our search for 3L1S models focused on static models in the Chang-Refsdal regime. Thus, some classes of solutions, such as those with significant orbital motion, may have been missed as a result. For example, we have not ruled out a system in which the Neptune might orbit a
<figure><img src="image_12.png" /><figcaption>**Figure 10.** Simulated 3L1S light-curve data with 2L1S models, similar to Figure 8, except that $\ln s_2 = 0.5$. For $n = 6$ and $n = 8$, small insets show the excursions of the data beyond the standardized window-size of the panel. The cases $n = 5, n = 6$, and $n = 8$ suggest that the light curve “factors” into contributions from two planets, as described by Han et al. (2001).</figcaption></figure>

very close binary whose orbital motion had significant effects on the light curve. Because the evidence for a second planet is not decisive, we have not exhaustively searched such models, which are substantially more complex.

As we will show, the main evidence against the 2L2S solutions is that they are inconsistent (or, at least, in very strong tension) with the limits on lens flux from the blended light. However, before investigating this issue, we must first assess to what extent the 3L1S solutions are consistent with the blended light.

## 6.3. Nature of the Blended Light

The blended light in this event represents something of a puzzle. We first assess this puzzle within the context of the 3L1S models and then discuss how this assessment is altered within the context of 2L2S models. The four facts that need to be evaluated are as follows. (1) The source contributes only a tiny fraction of the light from the “baseline object,” meaning that the “baseline object” can be effectively identified with the

blend. (2) The blend is astrometrically offset from the lens by <50 mas, but we concluded that the balance of evidence was that it was not coincident with the lens. (3) The predicted flux from the lens in the best-fit models is about 1 mag fainter than the blend, so the upper limits on lens light are clearly satisfied. (4) The blend is “relatively blue” in a sense that we will make clear shortly.

The blend must be either (1) the lens itself, (2) a companion to the lens, (3) a companion to the source, (4) an unrelated ambient star, or (5) combined light from two (or more) of the above, for which both light sources contribute significantly.

First, consider the scenario that the blend is behind the same dust column as the clump. Recall from Section 4.1.2 that $[(V-I), I]_{\text{base}} = (1.36, 18.60)$. The extinction toward the clump is $[E(V-I), A_I]_{\text{cl}} = (0.80, 0.97)$ (Nataf et al. 2013). If the blend is in the bulge or more than a few kiloparsecs from us in the disk, then it will suffer the same extinction. In this case, $[(V-I), I]_{0,\text{base}} = (0.56, 17.63)$. In particular, if the blend is in the bulge, then $[(V-I)_{0}, M_I]_{\text{base}} \simeq (0.56, 3.1)$.
This is quite consistent with a metal-poor turnoff star, particularly if we account for the roughly 0.05 mag combined color error from the blend measurement and the estimate of $E(V-I)$, as well as the 0.03 mag error in the blend-color estimate due to SBF (Section 4.1.3). Hence, considering the blend light in itself, by far the simplest explanation is that it is either a companion to the source or an ambient star in the bulge. The small astrometric offset between the blend and the lensing event makes the first possibility more likely, but the two possibilities are essentially the same in their implications for the lens system.

The problem with this hypothesis, however, is that the event models predict that the lens itself contributes a significant amount of light to the blend and that this light is red. This means that the “remaining light” from the blend must be even bluer than just deduced. That is, at $M = 0.25 M_\odot$ and $D_L \simeq 1$ kpc, the lens would have $((V-I)_0, I_0)_L \sim (3.0, 19.0)$. Assuming, for example, $(E(V-I), A_I)_L \sim (0.4, 0.5)$ (i.e., half the dust column in front of the lens), this implies $((V-I), I)_L \sim (3.4, 19.5)$. That is, the lens would contribute about 44% of the I-band flux, but only 7% of the V-band light. Hence, $((V-I), I)_{\text{remain}} \sim (0.8, 19.2)$. Thus, if this remaining light were at the distance to (i.e., extinction of) the lens, it would have $[(V-I)_0, M_I]_{\text{remain}} \sim (0.4, 8.7)$, while if it were at the distance to (i.e., the extinction of) the source, it would have $[(V-I)_0, M_I]_{\text{remain}} \sim (0, 3.7)$. Neither of these positions on the CMD correspond to any normal star. Nor would the problem be resolved by placing the “remaining light” somewhere else along the line of sight between the lens and the source. Of course, we have made these calculations using one particular estimate for the lens mass and hence $V/I$ fluxes. At $2\sigma$, the lens could be $M = 0.17 M_\odot$ (so, $((V-I)_0, I_0)_L \sim (3.5, 20.1))$, and at $1\sigma$ the baseline object (assumed to be behind the dust) could have dereddened values $[(V-I), I]_{0,\text{base}} = (0.62, 17.5)$. Then, the same reasoning given above would lead to $[(V-I)_0, M_I]_{\text{remain}} \sim (0.53, 3.1)$, which is consistent with a metal-poor turnoff star.

This means that there is only $2\sigma$ tension between the light-curve-based solution and the constraints from blended light. Such tension is fairly typical of microlensing planets that have undergone deep investigation. For example, Dong et al. (2009b) found $2\sigma$ tension between various lines of evidence from higher-order light-curve effects together with HST imaging when they reanalyzed OGLE-2005-BLG-071 (Udalski et al. 2005). Then, 14 yr after the event, Bennett et al. (2020) separately resolved the source and lens using Keck adaptive optics (AO), which showed that the true value of $\pi_{E,\perp}$ (the component of $\pi_E$ perpendicular to Earth’s acceleration) differed from the measured value by $2\sigma$ according to the error bars of the original estimate. See Figure 2 of Dong et al. (2009b) and Figure 2 of Bennett et al. (2020). That is, the original analysis reported a $2\sigma$ tension, and this was resolved by a $2\sigma$ correction to the original measurement.

Another possibility is that the blend is actually the lens. This would require two “adjustments” to the results of the measurements described in the body of this paper. First, the blend would have to be coincident with the lens rather than being displaced by $\lesssim 50$ mas. As we have discussed in some detail in Section 4.1.1, this is possible, although against the “balance of evidence.” Second, the lens would have to be substantially farther from us than the microlensing models predict. This is because the blue light of the blend (if attributed

<table><caption>**Table 7**

Physical Parameters for 2L2S Models</caption><thead><tr><th rowspan="2">Quantity</th><th colspan="2">Parallax</th></tr><tr><th>u<sub>0</sub> > 0</th><th>u<sub>0</sub> < 0</th></tr></thead><tbody><tr><td>M<sub>lens</sub> [M<sub>☉</sub>]</td><td>0.324<sup>+0.053</sup><sub>-0.041</sub></td><td>0.342<sup>+0.055</sup><sub>-0.046</sub></td></tr><tr><td>M<sub>planet</sub> [M<sub>☉</sub>]</td><td>8.389<sup>+1.629</sup><sub>-1.110</sub></td><td>8.058<sup>+1.462</sup><sub>-1.102</sub></td></tr><tr><td>a<sub>⊥</sub> [au]</td><td>1.536<sup>+0.226</sup><sub>-0.174</sub></td><td>1.598<sup>+0.219</sup><sub>-0.193</sub></td></tr><tr><td>a<sub>⊥</sub>/a<sub>snow</sub> [au]</td><td>1.757<sup>+0.066</sup><sub>-0.067</sub></td><td>1.726<sup>+0.064</sup><sub>-0.067</sub></td></tr><tr><td>D<sub>L</sub> [kpc]</td><td>1.008<sup>+0.167</sup><sub>-0.129</sub></td><td>1.034<sup>+0.153</sup><sub>-0.142</sub></td></tr><tr><td>μ<sup>geo</sup> [mas yr<sup>-1</sup>]</td><td>3.439<sup>+0.250</sup><sub>-0.190</sub></td><td>3.219<sup>+0.265</sup><sub>-0.204</sub></td></tr><tr><td>μ<sub>hel,N</sub> [mas yr<sup>-1</sup>]</td><td>-2.911<sup>+0.177</sup><sub>-0.215</sub></td><td>3.689<sup>+0.318</sup><sub>-0.249</sub></td></tr><tr><td>μ<sub>hel,E</sub> [mas yr<sup>-1</sup>]</td><td>1.344<sup>+0.317</sup><sub>-0.335</sub></td><td>1.535<sup>+0.342</sup><sub>-0.304</sub></td></tr><tr><td>ν<sub>lsr,l</sub> [km s<sup>-1</sup>]</td><td>1.782<sup>+2.166</sup><sub>-2.245</sub></td><td>34.370<sup>+2.358</sup><sub>-1.840</sub></td></tr><tr><td>ν<sub>lsr,b</sub> [km s<sup>-1</sup>]</td><td>-7.300<sup>+0.887</sup><sub>-1.199</sub></td><td>9.765<sup>+1.925</sup><sub>-1.424</sub></td></tr></tbody></table>

to a single star) would require that the star be much more massive than in the best solution, so much brighter, and therefore much farther in order to be consistent with the measured light from the blend.

Changing the mass and distance creates tension with the microlensing models because it requires changing either $\theta_E$ or $\pi_E$. The more plausible route is to decrease $\pi_E$ because (from Equation (10)) this both increases the mass and increases the distance as required by our hypothesis above. In addition, $\theta_E$ is much more robustly measured than $\pi_E$. In particular, the amplitude $\pi_E \equiv |\pi_E|$ is dominated by $\pi_{E,N}$, which has relatively large formal errors that are a direct result of the fact that it is more difficult to measure than $\pi_{E,E}$. One quickly finds that the lens must lie at least several kiloparsecs from us, so behind most or all of the dust. This implies $M \gtrsim 1 M_\odot$ to explain the color $(V-I)_{0,\text{base}} \simeq 0.56$. Then $\pi_{rel} = \theta_E^2/\kappa M \to 0.135 \text{ mas}(M/M_\odot)^{-1}$, where we have adopted $\theta_E = 1.05 \text{ mas}$ for reasons that we discuss below. Thus, for the case $M = 1 M_\odot$, $D_L = 3.8 \text{ kpc}$, and hence the measured blend flux would imply $M_I = 4.7$. This is somewhat dim for even a metal-poor turnoff star, but perhaps compatible given the uncertainties. The main remaining issue is that this lens mass requires $\pi_E = \theta_E/\kappa M \to 0.13$ compared to $\pi_E = 0.47 \pm 0.16$ for the lowest of the four parallax solutions in Table 2. However, given the relatively large errors, this “large” deviation may not be very seriously disfavored. Naively, it appears to be only a $2.1\sigma$ discrepancy. In order to test this possibility more rigorously, we force the “close” and “wide” ($u_0 < 0$) geometries by fixing $\pi_{E,N} = 0.12$. We find, indeed, that these solutions (Table 9) are only disfavored relative to the corresponding solutions in Table 2 by $\Delta\chi^2 \sim 3$, and therefore by only $\Delta\chi^2 \sim 7$ relative to the best solution. Hence, this resolution of the puzzle would have qualitatively similar tension to the one discussed above based on errors in the color and magnitude of the baseline object. Given that there are often systematic errors in microlensing at this level, these values do not in themselves disqualify this blend = lens solution. If future high-resolution imaging (see Section 6.4) shows that the blend is closely aligned with the lens, then this solution will become the most probable. If not, it will be ruled out.

Finally, another possibility is that the parallax has been underestimated by the fit, due to low-level systematics or contamination from unmodeled xallarap effects. For example, if the lens were 1.5 mag fainter than we have estimated from the best-fit microlensing model (i.e., using Equation (10)), then the lens would contribute only 11% to the I-band light, giving the “remaining light” a color of $(V-I)_{0,\text{remain}} \sim 0.46$. This is
<table><caption>**Table 8**

Physical Parameters for 2L2S Models with Source Orbital Motion</caption><thead><tr><th rowspan="2">Quantity</th><th colspan="2">Parallax</th></tr><tr><th>u<sub>0</sub> > 0</th><th>u<sub>0</sub> < 0</th></tr></thead><tbody><tr><td>M<sub>lens</sub> [M<sub>☉</sub>]</td><td>0.510<sup>+0.081</sup><sub>-0.053</sub></td><td>1.363<sup>+2.041</sup><sub>-0.631</sub></td></tr><tr><td>M<sub>planet</sub> [M<sub>☉</sub>]</td><td>9.434<sup>+1.748</sup><sub>-1.153</sub></td><td>34.924<sup>+52.062</sup><sub>-15.957</sub></td></tr><tr><td>a<sub>⊥</sub> [au]</td><td>2.121<sup>+0.303</sup><sub>-0.214</sub></td><td>4.587<sup>+2.589</sup><sub>-1.582</sub></td></tr><tr><td>a<sub>⊥</sub>/a<sub>snow</sub> [au]</td><td>1.533<sup>+0.039</sup><sub>-0.042</sub></td><td>1.245<sup>+0.278</sup><sub>-0.464</sub></td></tr><tr><td>D<sub>L</sub> [kpc]</td><td>1.263<sup>+0.212</sup><sub>-0.158</sub></td><td>3.062<sup>+1.732</sup><sub>-1.069</sub></td></tr><tr><td>μ<sub>geo</sub> [mas yr<sup>-1</sup>]</td><td>3.206<sup>+0.144</sup><sub>-0.150</sub></td><td>3.535<sup>+0.235</sup><sub>-0.121</sub></td></tr><tr><td>μ<sub>hel,N</sub> [mas yr<sup>-1</sup>]</td><td>-2.754<sup>+0.118</sup><sub>-0.103</sub></td><td>3.556<sup>+0.347</sup><sub>-6.585</sub></td></tr><tr><td>μ<sub>hel,E</sub> [mas yr<sup>-1</sup>]</td><td>0.748<sup>+0.263</sup><sub>-0.281</sub></td><td>-0.031<sup>+0.677</sup><sub>-1.312</sub></td></tr><tr><td>ν̃<sub>lsr,l</sub> [km s<sup>-1</sup>]</td><td>-2.439<sup>+2.670</sup><sub>-3.889</sub></td><td>64.693<sup>+44.091</sup><sub>-134.692</sub></td></tr><tr><td>ν̃<sub>lsr,b</sub> [km s<sup>-1</sup>]</td><td>-7.263<sup>+0.911</sup><sub>-1.442</sub></td><td>35.649<sup>+80.471</sup><sub>-47.020</sub></td></tr><tr><td>s<sub>z</sub> (source) [au]</td><td>0.274<sup>+0.018</sup><sub>-0.019</sub></td><td>0.330<sup>+0.027</sup><sub>-0.023</sub></td></tr><tr><td>period(source) [yr]</td><td>0.161<sup>+0.016</sup><sub>-0.017</sub></td><td>0.212<sup>+0.026</sup><sub>-0.021</sub></td></tr></tbody></table>

<table><caption>**Table 9**

3L1S Parallax-only Models for Fixed $\pi_{E,N}$</caption><thead><tr><th rowspan="2">Parameters</th><th>Close</th><th>Wide</th></tr><tr><th>$\pi_{E,N} = 0.12$</th><th>$\pi_{E,N} = 0.12$</th></tr></thead><tbody><tr><td>$\chi^2/\text{dof}$</td><td>9095.359/9120</td><td>9095.939/9120</td></tr><tr><td>$t_0$ (HJD')</td><td>8219.571 ± 0.007</td><td>8219.494 ± 0.018</td></tr><tr><td>$u_0$ (10<sup>-3</sup>)</td><td>-9.426 ± 0.586</td><td>-9.240 ± 0.488</td></tr><tr><td>$t_E$ (days)</td><td>121.167 ± 7.127</td><td>118.429 ± 5.344</td></tr><tr><td>$s_1$</td><td>1.014 ± 0.001</td><td>1.011 ± 0.001</td></tr><tr><td>$q_1$ (10<sup>-4</sup>)</td><td>1.371 ± 0.099</td><td>1.401 ± 0.083</td></tr><tr><td>$\alpha$ (rad)</td><td>0.462 ± 0.005</td><td>0.463 ± 0.005</td></tr><tr><td>$\rho$ (10<sup>-4</sup>)</td><td>3.618 ± 0.186</td><td>3.658 ± 0.174</td></tr><tr><td>$s_2$</td><td>0.425 ± 0.033</td><td>2.251 ± 0.636</td></tr><tr><td>$q_2$ (10<sup>-3</sup>)</td><td>2.179 ± 0.580</td><td>2.311 ± 0.638</td></tr><tr><td>$\psi$ (rad)</td><td>-0.160 ± 0.036</td><td>0.141 ± 0.034</td></tr><tr><td>$\pi_{E,E}$</td><td>-0.044 ± 0.015</td><td>-0.056 ± 0.015</td></tr><tr><td>$f_S$</td><td>0.0200 ± 0.0013</td><td>0.0205 ± 0.0011</td></tr><tr><td>$f_B$</td><td>0.7764 ± 0.0012</td><td>0.7759 ± 0.0010</td></tr><tr><td>$t_s$ (days)</td><td>0.044 ± 0.002</td><td>0.043 ± 0.002</td></tr></tbody></table>

still quite blue, but allowing for measurement errors, marginally acceptable, which would resolve the problems with supposing the blend is an unrelated star or a companion to the lens or source. The change in $\pi_{E,N}$ required to produce such a fainter lens would induce similar modest stress on the fit as the one proposed above. It is intrinsically less likely because there are fewer nearby stars in the observation cone than distant stars. However, if future high-resolution images show that the lens is displaced from the blend, then it will become more probable.

By contrast, the conflict between the predictions for the lens flux and the limits from the blended light is not so easily resolved for the 2L2S solutions. The problem is that the lens masses in these solutions are all substantially higher compared to the 3L1S solutions. We focus on the solutions with orbital motion (Table 8) because these are favored by $\Delta\chi^2 = 12$ and also because (as discussed in Section 3.6.1) Newton's laws imply that orbital motion should be significant for two bodies that are projected so close on the sky. However, we note that the issue is also very severe for the 2L2S models without orbital motion (Table 7).

For the $u_0 > 0$ solution, the best-fit host mass ($M = 0.51 M_☉$) at 1.2 kpc already implies $I_L \sim 17.3$. This substantially exceeds the estimate of the flux from the baseline object from Equation (6). This conflict is not substantially ameliorated by adopting the $2\sigma$ lower limit on the mass ($M \rightarrow 0.4 M_☉$), in which case $I_L \rightarrow 18$, which is still too bright. Moreover, even if one pushed the lens mass lower (so, in greater tension with the light-curve parameters), so as to be under the blended-light flux limit, the resulting lens would account for only a small fraction of the V-band light, so that the “remaining light” would have to be extraordinarily blue.

For the $u_0 < 0$ solution, the best-estimated lens mass is much larger, but so are the error bars. However, from a more detailed examination of the results of the Galactic model, we find that the lens mass is $M > 0.5 M_☉$ at $2\sigma$ and is $M > 0.4 M_☉$ at $3\sigma$. Moreover, these low-mass Galactic-model realizations also have small lens distances, $D_L \lesssim 1.5$ kpc. Hence, the same argument applies to this solution as to the $u_0 > 0$ solution.

In sum, if these arguments are taken at face value, we should regard the 2L2S solutions as excluded by the blended-light constraints assuming that the lens is a main-sequence (i.e., luminous) star. However, in view of the complexity of the modeling and the low level (scarcely visible by eye) of the extra-body (3L1S or 2L2S) effects, we adopt a more conservative viewpoint and consider the 2L2S solutions as “significantly disfavored” rather than “ruled out.”

## 6.4. Future High-resolution Imaging and Spectroscopy

Many puzzles remain for the OGLE-2018-BLG-0532L system, at least some of which can be clarified by high-resolution follow-up observations, either AO observations from the ground or possibly HST or James Webb Space Telescope observations from space.

The foremost question is the relation of the blended light to the event. Because the blended light is projected $\Delta\theta \lesssim 50$ mas from the lens, it is very likely to be associated with the event, i.e., either the lens itself, a companion to the lens, or a companion to the source. However, in principle, the blend could be an unrelated ambient star.

The most exciting possibility is that the blend is the lens. Recall that we concluded that this was against the “balance of evidence” but was still possible given the systematic astrometric errors induced by unresolved stars. A high-resolution image taken “immediately” (i.e., during 2019) could largely resolve this question. Recall that the astrometric position of the lens relative to the KMT template is known with an error of only 5 mas. The astrometric error of a high-resolution image would almost certainly be even smaller. Hence, the main uncertainty in identifying these two positions will simply be the displacement of the blend star relative to the field stars due to its unknown proper motion. Although this will already be on the order of 5 mas in 2019, and will continue to grow, in principle it can be corrected based on a subsequent high-resolution image. If this imaging shows a substantial offset, then the blend is definitely not the lens. And if the offset is within, say, 10 mas, then it is most likely the lens. In particular, such a small offset would decisively exclude the blend as a companion to the lens because such a close companion would have given rise to a huge microlensing signal near the peak of the event. And it would also virtually exclude the scenario that the blend is an ambient star because the probability for an $I < 19$ star to lie within 10 mas of a
predefined location is $<10^{-3}$. Hence, the only remaining
possible identifications for the blend would be the lens itself or
a companion to the source.

To better facilitate this, we provide a fits image of the KMTC template image together with an ascii file containing photometry and astrometry from this image. The source position (determined from difference images near peak) is (x, y) = (167.572, 151.009). The catalog can be directly cross-matched to another catalog in the relevant band to determine whether there is excess light at the position of the source.

Once the results of this early high-resolution imaging are known, then the next steps can be decided. Consequently, we do not attempt to chart these in excessive detail. However, we note that if the blended light proves to be the lens, then it will be possible to probe the planetary system using radial-velocity (RV) measurements with next generation telescopes. For example, in the ($s_2 < 1, u_0 < 0$) solution discussed in Section 6.3, the projected separation of the “Jupiter” would be ~1.7 au, while the mass ratio would be about $q_2 \sim 2.2 \times 10^{-3}$, implying a potential RV amplitude as high as ~50 m s$^{-1}$, with a period as short as $P = 2.2$ yr. Such measurements on an $I < 19$ star are plausibly feasible for 30 m telescopes. Of course, in the wide solutions, the maximum possible amplitude would be ~20 m s$^{-1}$, and the period would be on the order of 12 times longer. Nevertheless (again, if the blend were shown to be the lens), there would be at least some prospect for probing this system.

Work by A.G. was supported by AST-1516842 from the US NSF. I.G.S. and A.G. were supported by JPL grant 1500811. A.G. received support from the European Research Council under the European Union's Seventh Framework Programme (FP 7) ERC grant Agreement No. [321035] Work by C.H. was supported by the grant (2017R1A4A1015178) of National Research Foundation of Korea. 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 OGLE project has received funding from the National Science Centre, Poland, grant MAESTRO 2014/14/A/ST9/00121 to AU. This research used data obtained through the Telescope Access Program (TAP), which has been funded by the National Astronomical Observatories of China, the Chinese Academy of Sciences (the Strategic Priority Research Program “The Emergence of Cosmological Structures” grant No. XDB09000000), and the Special Fund for Astronomy from the Ministry of Finance. This work was partly supported by the National Science Foundation of China (grant No. 11333003, 11390372 to S.M.). Partly based on observations obtained with MegaPrime/MegaCam, a joint project of CFHT and CEA/DAPNIA, at the CanadaFranceHawaii Telescope (CFHT), which is operated by the National Research Council (NRC) of Canada, the Institut National des Science de l'Univers of the Centre National de la Recherche Scientifique (CNRS) of France, and the University of Hawaii. 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.

This work was performed in part under contract with the
California Institute of Technology (Caltech)/Jet Propulsion

Laboratory (JPL) funded by NASA through the Sagan Fellowship Program executed by the NASA Exoplanet Science Institute. M.T.P. was supported by NASA grants NNX14AF63G and NNG16PJ32C, as well as the Thomas Jefferson Chair for Discovery and Space Exploration.

Appendix
Decision to Search for 3L1S Solutions

As briefly telegraphed in Section 3.4, our decision to carry out the 3L1S investigation was prompted by “accidental” developments in the course of the 2L1S investigation, i.e., apparent “problems” in the 2L1S solution that were all eventually resolved. Such detours are normal in any relatively complex scientific study, but the “lab notebook details” describing them are usually omitted from journal papers describing the work, in order to avoid burdening the reader with the tribulations of the authors. In the present case, however, these “accidental” developments led the investigation in an unexpected new direction, which could have important implications for both the microlensing event described here and for future events in its class. For completeness, we therefore give a brief account of how we stumbled upon this direction.

As discussed in Section 3.1, parallax and orbital motion were introduced simultaneously, but the orbital motion only improved the fit by attaching a caustic feature to the baseline light curve. However, this was not recognized immediately, and the moderately significant differences between the “standard” and the “higher-order” solutions were attributed mainly to parallax. In particular, the biggest change was a much smaller value of $\rho$ (so larger $\theta_E$). When combined with the value of $\pi_E$ in that higher-order solution (about half the values in Table 1), this yielded a lens mass and distance (Equation (10)) that predicted a lens flux that was strongly inconsistent with the upper limits derived in Section 4.1.

We were led to investigate 3L1S solutions because they could solve this contradiction in one of two ways. First, if the host were actually a close binary, it could keep the same total mass but generate much less light because the total mass would be divided into two stars. Indeed this effect played a crucial role in the solution of OGLE-2007-BLG-349 (Bennett et al. 2016). Second, the tidal shear due to a wide binary (or corresponding quadrupole distortion of a close binary; Dominik 1999), could have affected the $\pi_{E,N}$ measurement, decreasing its amplitude relative to the true value, and so incorrectly raising the inferred mass.

However, the net result was mostly contrary to the expectations in the previous paragraph. Most importantly, the third body turned out to have $q \ll 1$, so in the close solution it did not play any role in reducing the total lens light. Second, the lens mass was actually reduced mostly by increasing $\rho$ (so decreasing $\theta_E$) back near the level of the “standard” 2L1S solution. In addition, the mass was further reduced by an increase in $|\pi_{E,N}|$, which was in accord with our naively reasoned expectations.

Because of these puzzling results (particularly the fact that the 3L1S “higher-order” solution looked more like the 2L1S “standard” solution than it looked like the 2L1S “higher-order” solution), we undertook an investigation that normally would have been made only at a later stage: examining cumulative $\Delta\chi^2$ plots between different solutions in order to locate the times of observations and the individual observatories that were contributing the most to distinguishing between solutions. It
was at this point that we discovered that the orbital-motion parameters were leading to spurious solutions through coupling to noise in the baseline. After removing these parameters, we obtained the solutions that are discussed in Section 3.1.

## ORCID iDs

Yoon-Hyun Ryu https://orcid.org/0000-0001-9823-2907

Andrzej Udalski https://orcid.org/0000-0001-5207-5619

Jennifer C. Yee https://orcid.org/0000-0001-9481-7123

Matthew T. Penny https://orcid.org/0000-0001-7506-5640

Weicheng Zang https://orcid.org/0000-0001-6000-3463

Michael D. Albrow https://orcid.org/0000-0003-3316-4012

Sun-Ju Chung https://orcid.org/0000-0001-6285-4528

Cheongho Han https://orcid.org/0000-0002-2641-9964

Kyu-Ha Hwang https://orcid.org/0000-0002-9241-4117

Youn Kil Jung https://orcid.org/0000-0002-0314-6000

In-Gu Shin https://orcid.org/0000-0002-4355-9838

Yossi Shvartzvald https://orcid.org/0000-0003-1525-5041

Hyoun-Woo Kim https://orcid.org/0000-0001-8263-1006

Seung-Lee Kim https://orcid.org/0000-0003-0562-5643

Chung-Uk Lee https://orcid.org/0000-0003-0043-3925

Richard W. Pogge https://orcid.org/0000-0003-1435-3053

Przemek Mróz https://orcid.org/0000-0001-7016-1692

Michał K. Szymański https://orcid.org/0000-0002-0548-8995

Jan Skowron https://orcid.org/0000-0002-2335-1730

Igor Soszyński https://orcid.org/0000-0002-7777-0842

Paweł Pietrukowicz https://orcid.org/0000-0002-2339-5899

Szymon Kozłowski https://orcid.org/0000-0003-4084-880X

Krzysztof Ulaczyk https://orcid.org/0000-0001-6364-408X

Patryk Iwanek https://orcid.org/0000-0002-6212-7221

Shude Mao https://orcid.org/0000-0001-8317-2788

## References

Alard, C., & Lupton, R. H. 1998, ApJ, 503, 325

Albrow, M. D., Horne, K., Bramich, D. M., et al. 2009, MNRAS, 397, 2099

Alonso-García, J., Mateo, M., Sen, B., et al. 2012, AJ, 143, 70

An, J. H., & Gould, A. 2001, ApJL, 563, L111

Bachelet, E., Shin, I.-G., Han, C., et al. 2012, ApJ, 754, 73

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

Batista, V., Gould, A., Dieters, S., et al. 2011, A&A, 529, 102

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., Rhie, S. H., Nikolaev, S., et al. 2010, ApJ, 713, 837

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

Bensby, T., Yee, J. C., Feltzing, S., et al. 2013, A&A, 549A, 147

Bessell, M. S., & Brett, J. M. 1988, PASP, 100, 1134

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

Claret, A. 2000, A&A, 363, 1081

Dominik, M. 1999, A&A, 349, 108

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., Bennett, D. P., Udalski, A., et al. 2008, Sci, 319, 927

Gaudi, B. S., Naber, R. M., & Sackett, P. 1998, ApJL, 502, L33

Gould, A. 1992, ApJ, 392, 442

Gould, A. 2000, ApJ, 542, 785

Gould, A. 2004, ApJL, 606, 319

Gould, A., Dong, S., Gaudi, B. S., et al. 2010, ApJ, 720, 1073

Gould, A., Udalski, A., An, D., et al. 2006, ApJL, 644, L37

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

Gould, A., Yee, J. C., Udalski, A., & Han, C. 2020, AcA, submitted

Griest, K., & Safizadeh, N. 1998, ApJ, 500, 37

Han, C., Chang, H.-Y., An, J. H., & Chang, K. 2001, MNRAS, 328, 986

Han, C., Udalski, A., & Gould, A. 2017, AJ, 154, 223

Han, C., Udalski, A., Lee Choi, J.-Y., et al. 2013, ApJ, 778, 38

Holtzman, J. A., Watson, A. M., Baum, W. A., et al. 1998, AJ, 115, 1946

Hwang, K.-H., Udalski, A., Bond, I. A., et al. 2018, AJ, 155, 259

Janczak, J., Fukui, A., Dong, S., et al. 2010, ApJ, 711, 731

Jung, Y. K., Gould, A., Zang, W., et al. 2019, AJ, 157, 72

Kervella, P., Thévenin, F., Di Folco, E., & Ségransan, D. 2004, A&A, 426, 297

Kim, D.-J., Kim, H.-W., Hwang, K.-H., et al. 2018, AJ, 155, 76

Kim, S.-L., Lee, C.-U., Park, B.-G., et al. 2016, JKAS, 49, 37

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

Mróz, P., Han, C., Udalski, A., et al. 2017, AJ, 153, 143

Muraki, Y., Han, C., Bennett, D. P., et al. 2011, ApJ, 741, 22

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

Paczyński, B. 1986, ApJ, 304, 1

Pál, A. 2012, MNRAS, 421, 1825

Park, B.-G., DePoy, D. L., Gaudi, B. S., et al. 2004, ApJ, 609, 166

Pascucci, I., Mulders, G. D., Gould, A., & Fernandes, R. 2018, ApJL, 856, L28

Poindexter, S., Afonso, C., Bennett, D. P., et al. 2005, ApJ, 633, 914

Poleski, R., Skowron, J., Udalski, A., et al. 2014, ApJ, 755, 42

Ranc, C., Bennett, D. P., Hirao, Y., et al. 2019, AJ, 157, 232

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

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

Skowron, J., Udalski, A., Gould, A., et al. 2011, ApJ, 738, 87

Smith, M., Mao, S., & Paczyński, B. 2003, MNRAS, 339, 925

Song, Y.-Y., Mao, S., & An, J. H. 2014, MNRAS, 437, 4006

Sumi, T., Bennett, D. P., Bond, I. A., et al. 2010, ApJ, 710, 1641

Sumi, T., Udalski, A., Bennett, D. P., et al. 2016, ApJ, 825, 112

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

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

Suzuki, D., Bennett, D. P., Udalski, A., et al. 2018b, AJ, 155, 263

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

Tomaney, A. B., & Crotts, A. P. S. 1996, AJ, 112, 2872

Udalski, A. 2003, AcA, 53, 291

Udalski, A., Jaroszyński, M., Paczyński, B., et al. 2005, ApJL, 628, L109

Udalski, A., Ryu, Y.-H., Sajadian, S., et al. 2018, AcA, 68, 1

Udalski, A., Szymanski, M., Kaluzny, J., et al. 1994, AcA, 44, 227

Udalski, A., Szymanski, M. K., Szymanski, G., et al. 2015, AcA, 65, 1

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

Woźniak, P. R. 2000, AcA, 50, 421

Wu, Y. 2019, ApJ, 874, 91

Yee, J. C., Han, C., Gould, A., et al. 2014, ApJ, 790, 14

Yee, J. C., Shvartzvald, Y., Gal-Yam, A., et al. 2012, ApJ, 755, 102

Yoo, J., DePoy, D. L., Gal-Yam, A., et al. 2004, ApJ, 603, 139

Zang, W., Penny, M. T., Zhu, W., et al. 2018, PASP, 130, 104401

Zhu, W., Gould, A., Penny, M., et al. 2014, ApJ, 794, 53