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

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, ≳100 days, and the peak of the main event is relatively sharp. Together, these characteristics would normally indicate a high-magnification event, A_max ≫ 1. However, the flux at the peak is only a factor of F_peak/F_base ∼ 3.6 above baseline. These two indicators can only be reconciled if the event is heavily blended, f_b ≫ 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 .

Zoom In Zoom Out Reset image size

The major anomaly, which is shown in the upper (zoomed) panel of Figure 1, is comprised of a trough that lasts ∼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 ≫ 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}_{{\rm{E}}},\rho ,\alpha )$ are allowed to vary in a Markov Chain Monte Carlo (MCMC). The first three (Paczyński 1986) parameters (t₀, u₀, t_E) are seeded at values derived from a 1L1S fit to the light curve with the anomaly removed. The angle α between the binary axis and lens-source relative proper motion ${{\boldsymbol{\mu }}}_{\mathrm{rel}}$ is seeded at six equally spaced points around the unit circle, while $\rho ={\theta }_{* }/{\theta }_{{\rm{E}}}$, i.e., the source radius normalized to the Einstein radius, is seeded at ρ = 5 × 10⁻⁴. 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}_{{\rm{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 linear limb-darkening coefficient Γ_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  ∼ 40. That is, despite the fact that the baseline appears relatively bright, ${I}_{\mathrm{base}}=18.24$ in the OGLE-IV system (see Section 4.1.2), the source is extremely faint, I_s  ∼ 22.25.

Table 1. Best-fit Solutions for 2L1S Models

		Parallax Models
Parameters	Standard	u0 > 0	u0 < 0
${\chi }^{2}/\mathrm{dof}$	9240.291/9124	9169.300/9122	9169.877/9122
t0 $({\mathrm{HJD}}^{{\prime} })$	8219.579 ± 0.004	8219.579 ± 0.004	8219.579 ± 0.004
u0 (10−3)	9.341 ± 0.371	8.295 ± 0.353	−7.873 ± 0.369
tE (days)	118.327 ± 4.709	133.269 ± 5.849	140.226 ± 5.962
s	1.014 ± 0.0005	1.012 ± 0.001	1.013 ± 0.001
q (10−4)	1.389 ± 0.053	1.015 ± 0.057	0.962 ± 0.056
α $(\mathrm{rad})$	−0.431 ± 0.002	−0.465 ± 0.005	0.465 ± 0.005
ρ (10−4)	3.154 ± 0.099	2.734 ± 0.105	2.680 ± 0.101
${\pi }_{{\rm{E}},N}$	⋯	−0.785 ± 0.105	0.765 ± 0.097
${\pi }_{{\rm{E}},E}$	⋯	−0.104 ± 0.018	−0.054 ± 0.015
fS	0.0201 ± 0.0008	0.0181 ± 0.0008	0.0171 ± 0.0008
fB	0.7764 ± 0.0008	0.7787 ± 0.0008	0.7797 ± 0.0008
t* (days)	0.037 ± 0.001	0.036 ± 0.001	0.038 ± 0.002

Note. ${t}_{* }\equiv \rho {t}_{{\rm{E}}}$ is not an independent quantity.

Download table as:  ASCIITypeset image

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

$$\begin{eqnarray}&&\begin{array}{l}{(V-I)}_{\mathrm{calib}}=0.93[{\left(V-I\right)}_{\mathrm{OGLE} \mbox{-} \mathrm{IV}}+0.309];\\ {I}_{\mathrm{calib}}={I}_{\mathrm{OGLE} \mbox{-} \mathrm{IV}}-0.005-0.008{\left(V-I\right)}_{\mathrm{calib}}.\end{array}\end{eqnarray} \tag{ 1 }$$

Thus, in the standard Johnson-Cousins system, ${I}_{\mathrm{base},\mathrm{standard}}\,=18.23$ and ${I}_{s,\mathrm{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.

Zoom In Zoom Out Reset image size

In spite of the faintness of the source (hence, high fractional errors $\sigma (F)/{f}_{s}$), the long Einstein timescale ${t}_{{\rm{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}_{{\rm{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 }_{{\rm{E}},N},{\pi }_{{\rm{E}},E})$, i.e., the components in equatorial coordinates of the "vector microlens parallax" ${{\boldsymbol{\pi }}}_{{\rm{E}}}$ (Gould 1992, 2000),

$$\begin{eqnarray}&&{{\boldsymbol{\pi }}}_{{\rm{E}}}=\displaystyle \frac{{\pi }_{\mathrm{rel}}}{{\theta }_{{\rm{E}}}}\,\displaystyle \frac{{{\boldsymbol{\mu }}}_{\mathrm{rel}}}{{\mu }_{\mathrm{rel}}},\end{eqnarray} \tag{ 2 }$$

where ${\theta }_{{\rm{E}}}$ is the Einstein radius,

$$\begin{eqnarray}&&{\theta }_{{\rm{E}}}=\sqrt{\kappa M{\pi }_{\mathrm{rel}}};\qquad \kappa \equiv \displaystyle \frac{4G}{{c}^{2}\,\mathrm{au}}\simeq 8.14\,\displaystyle \frac{\mathrm{mas}}{{M}_{\odot }},\end{eqnarray} \tag{ 3 }$$

M is the total mass of the lens, and ${\pi }_{\mathrm{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 }_{{\rm{E}}},\alpha ,d\alpha /{dt})\to -({u}_{0},{\pi }_{{\rm{E}}},\alpha ,d\alpha /{dt})$. We indeed find two such solutions with very similar χ², which are both shown in Table 1. We also search for jerk-parallax solutions (Gould 2004), which would approximately take ${\pi }_{{\rm{E}},N}\to -{\pi }_{{\rm{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 Δχ² = 71 is highly significant, and corresponding to this, the parallax amplitude ${\pi }_{{\rm{E}}}$ is measured to about 13%. By far the largest change, apart from the parallax, is the mass ratio, which drops by roughly 25% in the parallax solutions.

One point of possible concern is that ${\pi }_{{\rm{E}},E}$ is very close to zero, with very small measurement errors, while ${\pi }_{{\rm{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, ${{\boldsymbol{\pi }}}_{{\rm{E}}}\sim 0$, then this would be robustly reflected in ${\pi }_{{\rm{E}},\parallel }\simeq {\pi }_{{\rm{E}},E}$ but might be masked by subtle long-term systematics in the light curve that gave rise to a spurious ${\pi }_{{\rm{E}},\perp }\simeq {\pi }_{{\rm{E}},N}$.

However, first, it will be straightforward to show that the parallax must obey ${\pi }_{{\rm{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 }_{{\rm{E}},E}|$ is small compared to $| {\pi }_{{\rm{E}},N}|$, it is strongly inconsistent with zero, and quite consistent with the lower limit on ${\pi }_{{\rm{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 σ lower than the values reported in Table 1.

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 }_{{\rm{E}},{\rm{N}}}$ and ${\xi }_{{\rm{E}},{\rm{E}}}$), the orbital period P, and the ecliptic longitude (${l}_{\mathrm{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 χ² minimum near P ∼ 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.

Zoom In Zoom Out Reset image size

As a second test, we examine the best-fit coordinates for P = 1 yr. These coordinates are offset from the true coordinates by ∼27° and improve the fit by Δχ² = 19. In principle, a xallarap signal could have any set of coordinates over π sr, taking account of the four-fold symmetry of xallarap solutions. Thus, there is only a ∼20% probability of lying within 27° 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 measurement in Section 6.3.

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.

We begin by considering static 3L1S models with 10 geometric parameters, i.e., the seven "standard" 2L1S parameters (with (q,s) renamed (q₁, s₁)), plus three additional parameters, (q₂, s₂, ψ). 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₂ = 0.3 and s₂ = 0.05 and at several different values of ψ and then using MCMC to search for a local minimum. Such s₂ ≪  1 models are in the extreme Chang–Refsdal regime, in which the pseudo-shear ¹⁸ ${\gamma }_{2}^{c}\equiv {q}_{2}{s}_{2}^{2}$ is an approximate invariant. Therefore, one expects that (if the choice of ψ 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₂ and s₂ 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₂, s₂) minimum, probably due to the weakness of the χ² gradient. We therefore substituted ${s}_{2}\to {\gamma }_{2}^{c}$ as an MCMC variable, after which the minimum was found relatively quickly.

This search yields solutions with q₂ ≪ 1 and s₂ < 0.5 similar model light curve from the corresponding wide solution $({q}_{2},{s}_{2})\to ({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 }_{{\rm{E}},N},{\pi }_{{\rm{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 χ² 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, u₀ > 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 Δχ² = 81. This is formally very significant, but it would be by far the lowest Δχ² to justify an "officially accepted" planet. Therefore, some discussion is required to assess the reality of this putative planet. See Section 6.

Table 2. Best-fit Solutions for 3L1S Models

	Close			Wide
Parameters	u0 > 0	u0 < 0		u0 > 0	${u}_{0}\lt 0$
χ2/dof	9088.321/9119	9092.882/9119		9090.243/9119	9093.705/9119
t0 (HJD')	8219.590 ± 0.009	8219.571 ± 0.008		8219.479 ± 0.028	8219.493 ± 0.028
u0 (10−3)	8.228 ± 0.613	−8.893 ± 0.645		7.133 ± 0.503	−7.428 ± 0.709
tE (days)	139.317 ± 9.735	128.530 ± 9.344		151.039 ± 8.858	145.361 ± 11.011
s1	1.013 ± 0.001	1.013 ± 0.001		1.011 ± 0.001	1.011 ± 0.001
q1 (10−4)	0.975 ± 0.120	1.187 ± 0.133		0.927 ± 0.100	0.963 ± 0.168
α (rad)	−0.477 ± 0.007	0.475 ± 0.007		−0.478 ± 0.007	0.474 ± 0.009
ρ (10−4)	2.793 ± 0.240	3.310 ± 0.255		2.905 ± 0.198	2.884 ± 0.308
s2	0.364 ± 0.043	0.406 ± 0.041		2.656 ± 0.375	2.634 ± 0.396
q2 (10−3)	3.081 ± 0.871	2.231 ± 0.712		2.456 ± 1.059	2.208 ± 1.022
ψ (rad)	−0.036 ± 0.047	0.144 ± 0.042		−0.086 ± 0.047	0.084 ± 0.055
${\pi }_{{\rm{E}},N}$	−0.675 ± 0.129	0.473 ± 0.157		−0.630 ± 0.100	0.625 ± 0.227
${\pi }_{{\rm{E}},E}$	−0.105 ± 0.018	−0.041 ± 0.016		−0.087 ± 0.016	−0.038 ± 0.015
fS	0.0175 ± 0.0013	0.0189 ± 0.0013		0.0161 ± 0.0011	0.0168 ± 0.0015
fB	0.7790 ± 0.0013	0.7776 ± 0.0013		0.7804 ± 0.0011	0.7799 ± 0.0015
t* (days)	0.039 ± 0.002	0.043 ± 0.002		0.044 ± 0.002	0.042 ± 0.002

Note. ${t}_{* }\equiv \rho {t}_{{\rm{E}}}$ is not an independent quantity.

Download table as:  ASCIITypeset image

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). Whenever 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})\to ({t}_{\mathrm{0,1}},{u}_{\mathrm{0,1}})$, introduce a second pair $({t}_{\mathrm{0,2}},{u}_{\mathrm{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 Δχ² = 124 relative to the standard 2L1S solution. We then add the two parallax parameters, $({\pi }_{{\rm{E}},N},{\pi }_{{\rm{E}},E})$, which yields an additional ${\rm{\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 Δχ² = 16 relative to the 2L2S solution.

Table 3. Best-fit Solutions for 2L2S Models

		Parallax
Parameters	Standard	u0 > 0	u0 < 0
χ2/dof	9116.315/9121	9104.215/9119	9106.616/9119
${t}_{\mathrm{0,1}}$ $({\mathrm{HJD}}^{{\prime} })$	8219.598 ± 0.005	8219.594 ± 0.005	8219.594 ± 0.005
${u}_{\mathrm{0,1}}$ (10−3)	7.832 ± 0.429	5.926 ± 0.392	−5.557 ± 0.438
${t}_{{\rm{E}}}$ (days)	131.607 ± 6.148	173.809 ± 8.555	182.493 ± 11.562
s	1.014 ± 0.0006	1.013 ± 0.0006	1.014 ± 0.0006
q (10−4)	1.110 ± 0.060	0.705 ± 0.056	0.660 ± 0.064
α $(\mathrm{rad})$	−0.403 ± 0.004	−0.436 ± 0.007	0.431 ± 0.007
ρ1 (10−4)	2.747 ± 0.112	2.149 ± 0.122	1.986 ± 0.135
${t}_{\mathrm{0,2}}$ (HJD')	8217.655 ± 0.598	8219.763 ± 0.575	8220.006 ± 0.342
${u}_{\mathrm{0,2}}$	0.048 ± 0.007	0.026 ± 0.006	−0.023 ± 0.004
${qF},I$	0.137 ± 0.017	0.119 ± 0.018	0.130 ± 0.017
${\pi }_{{\rm{E}},N}$	⋯	−0.557 ± 0.095	0.527 ± 0.098
${\pi }_{{\rm{E}},E}$	⋯	−0.072 ± 0.020	−0.041 ± 0.015
fS	0.019 ± 0.0011	0.014 ± 0.0009	0.013 ± 0.0010
fB	0.777 ± 0.0010	0.782 ± 0.0009	0.783 ± 0.0010
t* (days)	0.036 ± 0.002	0.037 ± 0.002	0.036 ± 0.002

Note. ${t}_{* }\equiv \rho {t}_{{\rm{E}}}$ is not an independent quantity.

Download table as:  ASCIITypeset image

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 $\sim 0.02{\theta }_{{\rm{E}}}$, 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 Δχ² = 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_⊙ and 0.3 M_⊙ for the primary and secondary, respectively. We adopt a source distance of $8.0\,\mathrm{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}_{\mathrm{0,1}},{u}_{\mathrm{0,1}},{t}_{\mathrm{0,2}},{u}_{\mathrm{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, (α_s , δ_s ) describe the orientation of the source-orbit angular momentum vector, s_s is the source (3D) separation in units of ${\theta }_{{\rm{E}}}$, and ψ 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, 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 Δχ² = 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 Δχ² = 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 χ² difference is on the order of or smaller than typical systematic errors in microlensing. Hence the two classes of solutions cannot be distinguished based on goodness-of-fit criteria alone.

Table 4. Best-fit Solutions for 2L2S Models with Source Orbital Motion

Parameters	u0 > 0	u0 < 0
χ2/dof	9092.077/9117	9094.720/9117
${t}_{0,c}$ (HJD')	8219.683 ± 0.065	8219.461 ± 0.076
u0, c (10−2)	1.353 ± 0.065	−1.704 ± 0.012
${t}_{{\rm{E}}}$ (days)	185.023 ± 3.914	155.299 ± 5.269
s	1.014 ± 0.001	1.014 ± 0.001
q (10−4)	0.553 ± 0.018	0.750 ± 0.034
α (rad)	−0.487 ± 0.007	0.408 ± 0.017
ρ1 (10−4)	1.813 ± 0.060	2.234 ± 0.074
${s}_{{\rm{s}}}$	0.024 ± 0.002	0.029 ± 0.003
ψ (rad)	−0.362 ± 0.089	3.511 ± 0.059
${qF},I$	0.153 ± 0.013	0.154 ± 0.021
${\pi }_{{\rm{E}},N}$	−0.393 ± 0.063	0.140 ± 0.146
${\pi }_{{\rm{E}},E}$	−0.090 ± 0.018	−0.021 ± 0.017
αs (deg)	107.198 ± 5.614	84.458 ± 4.336
δs (deg)	40.817 ± 12.182	−29.558 ± 8.312
fS	0.012 ± 0.001	0.015 ± 0.001
fB	0.784 ± 0.001	0.781 ± 0.001
t* (days)	0.034 ± 0.002	0.035 ± 0.001

Note. ${t}_{* }\equiv \rho {t}_{{\rm{E}}}$ is not an independent quantity.

Download table as:  ASCIITypeset image

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}_{{\rm{E}}}$ and ${{qt}}_{{\rm{E}}}$ are approximate invariants (Yee et al. 2012), and the value of ${t}_{{\rm{E}}}$ is substantially higher in the 2L2S solution. Note in particular that the value of ${{\boldsymbol{\pi }}}_{{\rm{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.

We first note that the various light-curve models described above all obey ${F}_{\mathrm{base}}/{f}_{s}\gt 40$, where ${F}_{\mathrm{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 ≃ 045) 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.

Zoom In Zoom Out Reset image size

4.1.1. Astrometry

Here our goal is to find the astrometric offset between the microlensed 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 microlensed 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 $\sim 8/\sqrt{50}\sim 1\,$ mas. This error is negligible compared to the other errors, and so we ignore it. Hence, finally we derive

$$\begin{eqnarray}&&{\rm{\Delta }}{\theta }_{\mathrm{KMT}}(E,N)={{\boldsymbol{\theta }}}_{\mathrm{base}}-{{\boldsymbol{\theta }}}_{\mathrm{source}}=(-43,+20)\pm (6,4)\mathrm{mas}.\end{eqnarray} \tag{ 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

$$\begin{eqnarray}&&{\rm{\Delta }}{\theta }_{\mathrm{OGLE}}(E,N)=(-61,+31)\pm (5,5)\mathrm{mas}.\end{eqnarray} \tag{ 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 $+{{\boldsymbol{\mu }}}_{\mathrm{base}}{\rm{\Delta }}t/\eta$ relative to the true position of the source, where ${{\boldsymbol{\mu }}}_{\mathrm{base}}$ 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 ${\mu }_{\mathrm{base},{\rm{E}}}=+3.4\pm 1.5\,\mathrm{mas}\,{\mathrm{yr}}^{-1}$. 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 $| {\rm{\Delta }}{\boldsymbol{\theta }}| \lt 50\,\mathrm{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\,\mathrm{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)}_{\mathrm{base}}=(19.992,18.596)\pm (0.014,0.023)$ yields

$$\begin{eqnarray}&&{[(V-I),I]}_{\mathrm{base}}=(1.362,18.596)\pm (0.019,0.024).\end{eqnarray} \tag{ 6 }$$

4.1.3. Role of Surface-brightness Fluctuations

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})\gt 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}\lt 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 ± 0.05. To relate this measurement in K with the predicted limit in I, we write

$$\begin{eqnarray}&&{I}_{L}={K}_{L}+{(I-K)}_{0,L}+E{(I-K)}_{L}.\end{eqnarray} \tag{ 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_⊙, 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).

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 θ_* (Yoo et al. 2004). This in turn will yield the angular Einstein radius via the relation ${\theta }_{{\rm{E}}}={\theta }_{* }/\rho$, where ρ 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. ¹⁹ 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 θ_* 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,\mathrm{fid}}=0.020$, which corresponds to ${I}_{s,\mathrm{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,\mathrm{fid}}=22.23$, these imply an offset from the clump of ${\rm{\Delta }}[(V-I),I]=(0.66,6.95)\pm (0.11,0.08)$.

Zoom In Zoom Out Reset image size

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

We adopt ${[(V-I),I]}_{0,\mathrm{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 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

$$\begin{eqnarray}&&{\theta }_{* ,\mathrm{fid}}=0.388\pm 0.019\,\mu \mathrm{as}.\end{eqnarray} \tag{ 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 ≃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 θ_*. Including this additional 5%, the general formula for ${\theta }_{* }$ becomes

$$\begin{eqnarray}&&{\theta }_{* }=\sqrt{\displaystyle \frac{{f}_{s}}{0.020}}(0.388\pm 0.027\,\mu \mathrm{as}).\end{eqnarray} \tag{ 9 }$$

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

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 × 10⁻⁴. They had argued for a sharp break in the mass-ratio function at q_brk ∼ 0.56 × 10⁻⁴ 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 × 10⁻⁴. 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 ≲ 1 × 10⁻⁴ 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.

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 Δχ² ∼ 50–142 in three of these cases, which are similar to the Δχ² = 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 Δχ² = 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 Δχ² = 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 Δχ² = 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 Δχ² 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 ${\rm{\Delta }}{\chi }^{2}(t)={\chi }_{2{\rm{L}}1{\rm{S}}}^{2}(t)\,-{\chi }_{3{\rm{L}}1{\rm{S}}}^{2}(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 Δχ². 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. ²⁰

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

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 central caustic. At the time of the Shin et al. (2015) study, there were a total of nine "high-magnification" (u₀ < 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.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 Δχ² as OGLE-2018-BLG-0532. Further, one other had evidence at a lower level (Δχ² = 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 Δχ² = 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 (Δχ² < 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 and add these to models with the same q₂ but different s₂ and ψ compared to the best-fit model. The seven standard parameters $({t}_{0},{u}_{0},{t}_{{\rm{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 Δχ² 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 ${{\boldsymbol{\pi }}}_{{\rm{E}}}=0$. Conceivably, this might alter the ${\rm{\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 Δχ² 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 $\mathrm{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.65. In each case, we examine 10 values of ${\psi }^{{\prime} }={\psi }_{\mathrm{best}}-2\pi n/10$ with $n=0,1,\,\ldots 9$. For each simulation, we show the caustic geometry on the left and ${\rm{\Delta }}{\chi }^{2}$ on the right. We see that at lns = 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.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

By contrast, for lns₂ = 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 lns₂ = 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 lns₂ = 0.5 case, nearly all the examples have very obvious residuals that would almost certainly prompt investigation for an additional body or bodies in the system. We expect that for second planets even closer to the Einstein ring, 0 < lns₂ < 0.5, essentially all cases would yield obvious residuals that would have to be investigated. Further, due to the s ↔ s⁻¹ 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 lns (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, 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 Δχ² = 3.8. Given the level of systematics that is typical of microlensing light curves, this small χ² 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 ΔI ≃  1.9 mag fainter than the primary. Hence, one would expect that it would also be redder by ${\rm{\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 ${\rm{\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 ${\rm{\Delta }}(V-I)\,=2.5\mathrm{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 ${\rm{\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 Δχ² = 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 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.

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 $\lesssim 50\,\mathrm{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]}_{\mathrm{base}}=(1.36,18.60)$. The extinction toward the clump is ${[E(V-I),{A}_{I}]}_{\mathrm{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,\mathrm{base}}=(0.56,17.63)$. In particular, if the blend is in the bulge, then ${[{(V-I)}_{0},{M}_{I}]}_{\mathrm{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_⊙ and ${D}_{L}\simeq 1\,\mathrm{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)}_{\mathrm{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}]}_{\mathrm{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}]}_{\mathrm{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σ, the lens could be M = 0.17 M_⊙ (so, ${({(V-I)}_{0},{I}_{0})}_{L}\sim (3.5,20.1)$), and at 1σ the baseline object (assumed to be behind the dust) could have dereddened values ${[(V-I),I]}_{0,\mathrm{base}}=(0.62,17.5)$. Then, the same reasoning given above would lead to ${[{(V-I)}_{0},{M}_{I}]}_{\mathrm{remain}}\sim (0.53,3.1)$, which is consistent with a metal-poor turnoff star.

This means that there is only 2σ 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σ 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 }_{{\rm{E}},\perp }$ (the component of ${{\boldsymbol{\pi }}}_{{\rm{E}}}$ perpendicular to Earth's acceleration) differed from the measured value by 2σ 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σ tension, and this was resolved by a 2σ 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\,\mathrm{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 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 }_{{\rm{E}}}$ or ${\pi }_{{\rm{E}}}$. The more plausible route is to decrease ${\pi }_{{\rm{E}}}$ because (from Equation (10)) this both increases the mass and increases the distance as required by our hypothesis above. In addition, ${\theta }_{{\rm{E}}}$ is much more robustly measured than ${\pi }_{{\rm{E}}}$. In particular, the amplitude ${\pi }_{{\rm{E}}}\equiv | {{\boldsymbol{\pi }}}_{{\rm{E}}}|$ is dominated by ${\pi }_{{\rm{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 }_{{\rm{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 ≳ 1 M_⊙ to explain the color ${(V-I)}_{0,\mathrm{base}}\simeq 0.56$. Then ${\pi }_{\mathrm{rel}}={\theta }_{{\rm{E}}}^{2}/\kappa M\to 0.135\,\mathrm{mas}{(M/{M}_{\odot })}^{-1}$, where we have adopted ${\theta }_{{\rm{E}}}=1.05\,\mathrm{mas}$ for reasons that we discuss below. Thus, for the case M = 1 M_⊙, D_L  = 3.8 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 }_{{\rm{E}}}={\theta }_{{\rm{E}}}/\kappa M\to 0.13$ compared to ${\pi }_{{\rm{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) geometries by fixing ${\pi }_{{\rm{E}},N}=0.12$. We find, indeed, that these solutions (Table 9) are only disfavored relative to the corresponding solutions in Table 2 by Δχ² ∼ 3, and therefore by only Δχ² ∼ 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.

Table 9. 3L1S Parallax-only Models for Fixed ${\pi }_{{\rm{E}},N}$

	Close	Wide
Parameters	${\pi }_{{\rm{E}},N}=0.12$	${\pi }_{{\rm{E}},N}=0.12$
${\chi }^{2}/\mathrm{dof}$	9095.359/9120	9095.939/9120
t0 $({\mathrm{HJD}}^{{\prime} })$	8219.571 ± 0.007	8219.494 ± 0.018
u0 $({10}^{-3})$	−9.426 ± 0.586	−9.240 ± 0.488
${t}_{{\rm{E}}}$ $(\mathrm{days})$	121.167 ± 7.127	118.429 ± 5.344
s1	1.014 ± 0.001	1.011 ± 0.001
q1 $({10}^{-4})$	1.371 ± 0.099	1.401 ± 0.083
α $(\mathrm{rad})$	0.462 ± 0.005	0.463 ± 0.005
ρ $({10}^{-4})$	3.618 ± 0.186	3.658 ± 0.174
s2	0.425 ± 0.033	2.251 ± 0.636
q2 $({10}^{-3})$	2.179 ± 0.580	2.311 ± 0.638
ψ $(\mathrm{rad})$	−0.160 ± 0.036	0.141 ± 0.034
${\pi }_{{\rm{E}},E}$	−0.044 ± 0.015	−0.056 ± 0.015
fS	0.0200 ± 0.0013	0.0205 ± 0.0011
fB	0.7764 ± 0.0012	0.7759 ± 0.0010
${t}_{* }$ $(\mathrm{days})$	0.044 ± 0.002	0.043 ± 0.002

Download table as:  ASCIITypeset image

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,\mathrm{remain}}\sim 0.46$. This is 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 }_{{\rm{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 Δχ² = 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 solution, the best-fit host mass (M =0.51 M_⊙) at $1.2\,\mathrm{kpc}$ already implies I_L  ∼ 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σ lower limit on the mass ($M\to 0.4\,{M}_{\odot }$), in which case ${I}_{L}\to 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 the that the "remaining light" would have to be extraordinarily blue.

For the u₀ < 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σ and is M > 0.4 M_⊙ at 3σ. Moreover, these low-mass Galactic-model realizations also have small lens distances, D_L  ≲  1.5 kpc. Hence, the same argument applies to this solution as to the u₀ > 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."

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 Δθ ≲  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⁻³. 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₂ < 1, u₀ < 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.2 × 10⁻³, implying a potential RV amplitude as high as $\sim 50\,{\rm{m}}\,{{\rm{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 $\sim 20\,{\rm{m}}\,{{\rm{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.