© ESO 2022

1 Introduction

Studies based on radial velocity (RV) observations have shown that beyond ~1 AU about 10% of stars have giant planets (Cumming et al. 2008; Fulton et al. 2021). One question to ask is how common it is for such cold giant planets to have massive planetary-mass companions. On one hand, cold giant planets have eccentric orbits that are generally attributed to dynamical interactions with additional massive companions, suggesting that perhaps the majority of them have such companions, at least at birth (Jurić & Tremaine 2008; Chatterjee et al. 2008). On the other hand, after over two decades of searches, RV surveys have only been able to identify the presence of massive companions to ~20–30% of known cold giants (Wright et al. 2009; Rosenthal et al. 2021). This discrepancy can potentially be reduced with more discoveries of giant planet systems, but unfortunately RV becomes extremely inefficient in detecting planets with relatively low masses and/or long orbital periods.

Being most sensitive to cold planets located around or beyond the water snow line, gravitational microlensing can play an important role in completing the demographic census of exoplanets (Gaudi 2012; Zhu & Dong 2021). In terms of multi-planet systems, microlensing has so far detected three reliable two-planet systems, OGLE-2006-BLG-109L (Gaudi et al. 2008; Bennett et al. 2010), OGLE-2012-BLG-0026L (Han et al. 2013; Beaulieu et al. 2016), and OGLE-2018-BLG-1011L (Han et al. 2019), in addition to two candidate two-planet systems, OGLE-2014-BLG-1722L (Suzuki et al. 2018) and KMT-2019-BLG-1953L (Han et al. 2020a). Compared to the total number of over 100 planetary systems found with microlensing, the number of multi-planet systems is small. However, given the relatively low efficiency of detecting multi-planet systems with microlensing (Zhu et al. 2014), these numbers already somewhat indicate that a perhaps substantial fraction of microlensing planets have additional planetary-mass companions (Madsen & Zhu 2019).

Given the potential of microlensing in studying the multiplicity distribution of cold planets and thus the architecture of planetary systems in the cold region, it is important to detect more, secure multi-planet microlensing systems. This requires high-cadence observations over a large number of microlensing events as well as detailed light curve modelings of all anomalous events. The signal produced by multiple planets differs from that produced by a single planet because the individual planets induce their own caustics and these caustics often result in a complex pattern due to the interference between the caustics (Daněk & Heyrovský 2015a,b, 2019). As a result, these signals, in most cases, cannot be described by the usual lensing models based on the binary-lens (2L1S) or binary-source (1L2S) interpretations. This implies that some anomalous events with signals produced by multiple planets are probably left unanalyzed.

The amount of microlensing data was dramatically decreased during the COVID-19 pandemic because many of the major survey telescopes were shut down. In order to make the best use of this downtime, we conducted a project in which previous microlensing data collected by the Korea Microlensing Telescope Network (KMTNet) survey in and before the 2019 season were systematically reinvestigated. The aim of the project is to find events of scientific importance among those with no presented analyses. One group of events for this investigation are those with weak anomalies. Investigating such events led to the discoveries of 16 microlensing planets: KMT-2018-BLG-1025Lb (Han et al. 2021e), KMT-2016-BLG-2364Lb, KMT-2016-BLG-2397Lb, OGLE-2017-BLG-0604Lb, OGLE-2017-BLG-1375Lb (Han et al. 2020e), KMT-2018-BLG-0748Lb (Han et al. 2020d), KMT-2019-BLG-1339Lb (Han et al. 2020b), KMT-2018-BLG-1976, KMT-2018-BLG-1996, OGLE-2019-BLG-0954 (Han et al. 2021d), OGLE-2018-BLG-0977Lb, OGLE-2018-BLG-0506Lb, OGLE-2018-BLG-0516Lb, OGLE-2019-BLG-1492Lb, KMT-2019-BLG-0253 (Hwang et al. 2022), and OGLE-2019-BLG-1053 (Zang et al. 2021).

Another group of targeted events are those with known anomalies but for which interpretations of the anomalies have not been presented. From the investigation of such events, it was found that KMT-2019-BLG-1715 was a planetary event involved with three lens masses and two source stars (Han et al. 2021c), KMT-2019-BLG-0797 was a binary-lensing event occurring on a binary stellar system, a 2L2S event (Han et al. 2021b), and KMT-2019-BLG-1953 was a strong candidate planetary event with a lens composed of two planets and the host (Han et al. 2020a). The events in this group have the common characteristic that interpreting the lensing light curves of the events requires one to add extra source or lens components into the modeling.

In this work we present the results found from the reanalysis of the lensing event OGLE-2019-BLG-0468. The light curve of the event was previously investigated with 2L1S and 1L2S interpretations, but no plausible solution was suggested. From the reanalysis of the event based on more sophisticated models, we find that the event was produced by a triple-lens (3L1S) system, in which the lens is composed of two giant planets and their host star.

We present the analysis of the event according to the following organization. In Sect. 2 we describe observations of the lensing event and the characteristics of the observed light curve. In Sect. 3 we depict various models tested to explain the observed light curve, including 2L1S, 1L2S, 2L2S, and 3L1S models. In Sect. 4 we characterize the source star and estimate the angular Einstein radius. In Sect. 5 we estimate the physical lens parameters using the available observables of the event. In Sect. 6 we discuss the possibility that the baseline object is the lens. We then summarize our results and conclude in Sect. 7.

2 Observations and data

The source of the lensing event OGLE-2019-BLG-0468 lies in the Galactic bulge field at the equatorial coordinates ${(\text{RA},\text{Dec})}_{\text{J2000}}=($17:45:37.44, − 24:26:50.2), which correspond to the Galactic coordinates (l, b) = (3°.834,  2°.336). The apparent baseline I-band magnitude of the source is I_base = 17.8 according to the Optical Gravitational Lensing Experiment IV (OGLE-IV) photometry system. As we show in Sect. 4, the source is much fainter than the baseline magnitude, and the baseline flux comes mostly from a blend.

Figure 1 shows the light curve of OGLE-2019-BLG-0468. The rising of the source flux induced by lensing was first found by the OGLE-IV (Udalski et al. 2015) survey in the early part of the 2019 season, on April 13, 2019 (HJD′≡HJD − 2 450 000 ~ 8586). The OGLE teamutilizes the 1.3 m telescope at the Las Campanas Observatory in Chile, which is equipped with a camera that yields a 1.4 deg² field of view. The source flux increased until t₁ ~ 8586.3 and then declined during 8586 ≲HJD′≲ 8590, producing a weak bump at around t₁. The flux suddenly increased at t₂ ~ 8592.9, suggesting that the source crossed a caustic induced by the multiplicity of the lens. The detailed structure of the light curve for the three nights during the period 8596 ≲HJD′≲ 8598 could not bedelineated because the OGLE observation of the event was not conducted for that time. When the event was observed by the OGLE survey againon HJD′~ 8599, the source flux continued to decline until it reached the baseline.

The event was also located in the field covered by the KMTNet survey (Kim et al. 2016). The KMTNet group utilizes three identical telescopes, each of which has a 1.6 m aperture and is mounted with a camera that yields a 4 deg² field of view. For the continuous coverage of lensing events, the telescopes are distributed over the three continents of the Southern Hemisphere, at the Siding Spring Observatory in Australia (KMTA), the Cerro Tololo Inter-American Observatory in Chile (KMTC), and the South African Astronomical Observatory in South Africa (KMTS). For both the OGLE and KMTNet surveys, observations of the event were done mainly in the I band, and a fraction of V -band images wereacquired for the measurement of the source color. The event was identified by the KMTNet survey from the post-season inspection of the 2019 season data, and it was labeled as KMT-2019-BLG-2696. Hereafter we use only the OGLE event number, according to the order of discovery, to designate the event. The light curve constructed with the additional KMTNet data revealed that there existed an additional peak at t₃ ~ 8595.7 that was not covered by the OGLE data. Figure 2 shows the zoomed-in view of the light curve around the epochs t₂ and t₃.

The light curve of the event was constructed by conducting photometry of the source using the pipelines of the individual survey groups: Woźniak (2000) for OGLE and Albrow et al. (2009) for KMTNet. In order to consider the scatter of the data and to make χ² per degree of freedom for each data set unity, the error bars estimated by the pipelines were readjusted according to the prescription depicted in Yee et al. (2012).

In addition to the photometric data, we also obtained two spectra, with a 1000 s exposure for each, of the baseline object on the night of June 3, 2021(HJD′ = 9398), which is ~2.2 yr after the event, using the Robert Stobie Spectrograph (Burgh et al. 2003) mounted on the South African Large Telescope (SALT; Buckley et al. 2006). The spectroscopic data were reduced using a custom pipeline based on the PySALT package (Crawford et al. 2010), which accounts for basic charge-coupled device characteristics, removal of cosmic rays, wavelength calibration, and relative flux calibration. To estimate the stellar parameters (T_eff, log g, and [Fe/H]), we interpolated the observed spectra in an empirical grid (Du et al. 2019), which was constructed with massive spectra collected by the Large Sky Area Multi-Object Fibre Spectroscopic Telescope (LAMOST). Unfortunately, it was difficult to securely estimate the stellar parameters due to the low signal-to-noise ratios of the spectra.

Fig. 1 Light curve of the microlensing event OGLE-2019-BLG-0468. The arrows labeled t1 (8586.3), t2 (8592.9), and t3 (8595.7) indicate the three epochs of major anomalies. The curves drawn over the data points represent the 1L1S (dotted curve) and 3L1S (close-close model, solid curve) models. The zoomed-in view of the anomaly region around t2 and t3 is shown in Fig. 2.

Fig. 2 Zoomed-in view of the anomaly region around the epochs t2 and t3. Notations are the same as in Fig. 1.

3 Light curve interpretation

3.1 2L1S, 1L2S, and 2L2S models

We first modeled the observed light curve under the assumption of lens (2L1S model) or source (1L2S model) binarity, which is the most common cause of lensing light curve anomalies. The lensing light curve of a single-lens single-source (1L1S) event is characterized by the three parameters (t₀, u₀, t_E), which represent the time of the closest lens-source approach, the lens-source separation at t₀, and the event timescale, respectively. Adding an extra lens or source component into the modeling requires one to include extra parameters. For a 2L1S event, these extra parameters are (s, q, α), which denote the binary separation, the mass ratio between the lens components, and the angle between the direction of the source motion and the binary axis (source trajectory angle), respectively. For a 1L2S event, the extra parameters are (t_0,2, u_0,2, q_F), the first two of which are the closest approach time and separation of the source companion, and the last parameter indicates the flux ratio between the companion (S₂) and primary (S₁) source stars. In all cases of the tested models, we included an additional parameter, ρ, which denotes the ratio of the angular source radius, θ_*, to the angular Einstein radius, θ_E, that is, ρ = θ_*∕θ_E (normalized source radius), to account for possible finite-source effects in the lensing light curve. To distinguish parameters related to S₁ and S₂, we use the notations (t_0,1, u_0,1, ρ₁) and (t_0,2, u_0,2, ρ₂), respectively.

From the modeling of the observed light curve with the 2L1S and 1L2S interpretations, it was found that the data cannot be explained by these models. In Fig. 3 we plot the model curves and residuals of the 2L1S (dashed curve) and 1L2S (dot-dashed curve) models. The lensing parameters of these models are listed in Table 1 together with the χ² values of the fits.

We made an additional check with a 2L2S model, in which both the lens and source are binaries; such a system is exemplified by MOA-2010-BLG-117 (Bennett et al. 2018), OGLE-2016-BLG-1003 (Jung et al. 2017), KMT-2019-BLG-0797 (Han et al. 2021b), and KMT-2018-BLG-1743 (Han et al. 2021a). The model curve and residual from the 2L2S solution are presented in Fig. 3, and the lensing parameters of the solution are listed inTable 1. This model provides a better fit than the 2L1S and 1L2S models by Δχ² = 295.2 and 2069.7, respectively. However, the model still leaves substantial residuals, indicating that a new interpretation of the light curve is needed.

3.2 3L1S model

We then tested a model in which the lens is a triple system (3L1S model). In the first step of this modeling, we checked whether a 2L1S model can describe part of the anomalies, although it turned out that the model could not explain all anomaly features. We did this check because lensing light curves with three lens components (M₁, M₂, and M₃) can, in many cases, be approximated by the superposition of two 2L1S light curves produced by M₁ –M₂ and M₁ –M₃ lens pairs (Bozza 1999; Han et al. 2001). This is exemplified by the light curve of OGLE-2012-BLG-0026 (Han et al. 2013), which is generated by a lens composed of two planets and a host star, and that of OGLE-2018-BLG-1700 (Han et al. 2020c), which is produced by a lens composed of a planet and binary stars. From the modeling conducted with the partial data, excluding those after t₂, it is found that the light curve is well approximated by a 2L1S model with binary parameters (s, q, α) ~ (0.86, 3 × 10⁻³, 154°). This suggests the possibility that the event may be produced by a triple-lens (3L) system

Modeling with the inclusion of a third lens component requires one to add three extra parameters to those of the 2L1S model. These parameters are (s₃, q₃, ψ), which represent the M₁ –M₃ separation, the M₃∕M₁ mass ratio, and the orientation angle of M₃ as measured from the M₁ –M₂ axis with the center at the position of M₁, respectively.In the first round of the 3L1S modeling, we searched for the lensing parameters related to M₃, that is, s₃, q₃, and ψ, with a grid approach by fixing the parameters related to M₂, that is, s₂, q₂, and α, as the ones obtained from the preliminary 2L1S modeling conducted with the use of the partial data. Figure 4 shows the Δχ² distribution in the log s₃– logq₃ parameter plane obtained from this first-round modeling. We note that there exist two locals in the Δχ² map, indicating that there are degenerate solutions: at (logs₃, logq₃) ~ (+0.15, −2.0) and ~(−0.15, −2.0). We discuss this degeneracy more below. In the second round, we refined the local solutions found from the first-round modeling using a downhill approach based on the Markov chain Monte Carlo (MCMC) method. In this process, we released all parameters as free parameters.

We find that the anomaly features of the observed light curve are well explained by a 3L1S model. The model curve is plotted over the data points in Figs. 1 and 2, and the residual from the model is compared to those of the other models (1L2S, 2L1S, and 2L2S models) in Fig. 3. We find two solutions, in which s₂ < 1.0 and s₃ < 1.0 for one solution and s₂ < 1.0 and s₃ > 1.0 for the other solution. We refer to the individual solutions as “close-close” and “close-wide” solutions, respectively.

The lensing parameters of the two 3L1S solutions are listed in Table 2. Also listed in the table are the χ² values of the fits and the flux parameters of the source, f_s, and the blend, f_b, as measured on the OGLE flux scale. We note that the estimated source flux, f_s ~ 0.03, is much smaller than the blend flux, f_b ~ 1.16, indicating that the observed flux is heavily blended. The estimated mass ratios of the companions to the primary lens are q₂ ~ (3.3−3.5) × 10⁻³ and q₃ ~ 10.5 × 10⁻³, both of which correspond to the ratio between a giant planet and a star. This indicates that the lens is a planetary system composed of two giant planets. We note that the designation of q₂ and q₃ is not arranged according to the mass, and it turns out that q₂ < q₃. We find that the fit of the 3L1S model is better than those of the 1L2S, 2L1S, and 2L2S models by Δχ² = 2750.0, 974.5, and 680.3, respectively. The degeneracy between the close-close and close-wide solutions is very severe, and the close-close solution is preferred over the close-wide solution by merely Δχ² = 0.6. The fact that the M₁ –M₃ separations of the two degenerate solutions are in the relation of s_{3,close-close} × s_3,close-wide ~ 1 indicates that the degeneracy is caused by the close-wide degeneracy in s₃ (Griest & Safizadeh 1998; Dominik 1999).

Figure 5 shows the lens system configurations of the two 3L1S models. For each model, the caustic configuration (red figures)appears to be the superposition of the two sets of caustics induced by two low-mass companions (blue dots labeled M₂ and M₃ in the inset of each panel) of the primary (M₁). We mark the positions of the source at the three epochs of t₁, t₂, and t₃ (small empty circles drawn in magenta). The source position at t₁ corresponds to the source passage through the positive deviation region that extends from one of the sharp cusps of the caustic induced by M₂, and this produced a weak bump at around t₁. The source then crossed the tip of the caustic, producing a sharp caustic-crossing feature, and the two data points at t₂, one from OGLE and the other from KMTC observations, correspond to the time of the caustic entrance. Another caustic crossing occurred when the source passed the tip of the caustic induced by M₃, and this produced a sharp peak around t₃. The last peak was covered by the KMTC data, both on the rising and falling sides of the peak.

The insets of the panels in Fig. 5 show that the source passes close to one of the triangular planetary caustics induced by M₂. If the separation between the source and the caustic were very close, this source approach would induce a low bump at the time of the caustic approach. We checked this possibility by inspecting the magnification pattern around the planetary caustic. The magnification map around the caustics (for the close-wide 3L1S model) is shown in Fig. 6, in which the source location at the time of the caustic approach is marked by an empty magenta circle. The source size, ρ ~ 0.5 × 10⁻³, is too small to be seen in the presented scale, and thus we arbitrarily set the circle size. The lower panel shows the light curve around the time of the caustic approach at t_approach ~ 8581.7. Both the magnification map and the light curve show that the deviation at around this approach is too weak to emerge from the baseline.

The fact that the caustic configuration of the 3L1S solutions appears to be the superposition of two 2L1S caustics, together with the fact that M₃ is in the planetary-mass regime, suggests that there may be additional degenerate solutions caused by the close-wide degeneracy in s₂. We, therefore, searched for additional solutions resulting from this degeneracy: solutions with s₂ > 1.0, s₃ < 1.0 (“wide-close” solution), and s₂ > 1.0, s₃ > 1.0 (“wide-wide” solution). These degenerate solutions were found using the initial parameters of s_2,wide,xx = 1∕s_2,close-xx. The lensing parameters of these solutions are listed in Table 2. We find that, although the wide-xx solutions approximately describe the observed light curve, their fits are worse than the close-xx solutions by Δχ² ~ 14; the main contributions to this χ² difference come from the three data points (one from OGLE and two from KMTC) taken at HJD′~ 8594.7 and the seven data points (from KMTC) covering the anomaly at t₃. The superiority of the close-xx solutions over the wide-xx solutions is found from the comparison of the residuals of the models, presented in Fig. 7, as well as the cumulative distributions of the χ² difference from the best-fit solution (close-close solution), presented in Fig. 8. The degeneracy in s₂ is resolvable because the separation (s₂ ~ 0.85 for the close solution and ~ 1.1 for the wide solution) is close to unity. In this case, the M₁ –M₂ pair forms aresonant caustic, for which the close-wide degeneracy is less severe (Chung et al. 2005).

We checked whether the parameters of the normalized source radius, ρ, and the microlens-parallax, π_E, can be measured. Measurements of these parameters are important because the mass, M, and distance to the lens, D_L, can be uniquely determined with these parameters by $$M=\frac{{\theta }_{\text{E}}}{\kappa {\pi }_{\text{E}}};D_{\text{L}}=\frac{\text{AU}}{{\pi }_{\text{E}}{\theta }_{\text{E}}+{\pi }_{\text{S}}};{\theta }_{\text{E}}=\frac{{\theta }_{*}}{\rho },$$(1)

where κ = 4G∕(c²AU) and π_S denotes the annual parallax of the source (Gould 2000). The value of the normalized source radius, ρ ~ 0.5 × 10⁻³, is firmly measured from the data points involved with caustic crossings at around t₂ and t₃. On the other hand, modeling that considers the microlens-parallax effect results in a π_E value with aconsiderable uncertainty due to the fact that the photometry quality of the data is mediocre because of the faintness of the source, together with the fact that all the main features occur in an interval of 12 days, which is too short to appreciate any deviations due to parallax effects. Figure 9 shows the scatter plot of points in the MCMC chain on the π_E,E–π_E,N plane, where π_E,E and π_E,N are the east and north components of the microlens-parallax vector π_E, respectively. The scatter plot shows that the microlens-parallax is consistent with a zero-parallax model within the 2σ level, and the uncertainty of π_E,N is big, although the uncertainty of π_E,E is relatively small. Considering that the anomaly features in the light curve occurred within a short time interval, the lens-orbital motion would not invoke any significant effects. Nevertheless, we conducted an additional modeling considering the lens-orbital model and found that the fit improvement by the lens-orbital effect is negligible and that the orbital motion was poorly constrained.

Fig. 3 Comparison of the four tested models, 1L2S, 2L1S, 2S2L, and 3L1S (close-close solution). The lower panels show the residuals from the individual models.

Fig. 4 Distribution of Δχ2 in the logs3– logq3 parameter plane obtained from the preliminary 3L1S grid searches for the parameters related to the third lens component. Red, yellow, green, cyan, and blue are used to denote points with Δχ2 < n(12), < n(22), < n(32), < n(42), and < n(52), respectively, where n = 10.

Fig. 5 Lens system configurations of the close-close (left panel) and close-wide (right panel) 3L1L models. The inset in each panel shows the locations of the three lens components, denoted by blue dots labeled M1, M2, and M3. The dashed circle in the inset represents the Einstein ring. The line with an arrow represents the source trajectory. The three small empty dots on the source trajectory represent the source positions at the epochs of t1, t2, and t3 that are marked in Fig. 1. The size of the dots is not scaled to the source size. The red figure represents the caustic, and the gray curves around the caustic represent equi-magnification contours.

Fig. 6 Magnification pattern around the caustics (for the close-wide 3L1S model) and the light curve around the time of the source approach close to the planetary caustic induced by M2. The small empty circle on the source trajectory represents the source position at the time of the source approach, at tapproach ~ 8581.7 in HJD′. The size of the circle is arbitrarily set and is not scaled to the source size.

Fig. 7 Comparison of the residuals from the four degenerate 3L1S solutions: close-close, close-wide, wide-close, and wide-wide solutions. As a representative model, the model curve of the close-close solution is drawn over the data points in the top panel.

Fig. 8 Cumulative distributions of the Δχ2 difference from the best-fit solution (close-close model). The light curve in the upper panel is inserted to show the locations of the fit difference.

Fig. 9 Scatter plot of points in the MCMC chain on the πE,E–πE,N parameter plane for the close-close 3L1S solution. The color coding is the same as in Fig. 4, except n = 1.

4 Source star and angular Einstein radius

For the determination of the angular Einstein radius, we needed to estimate the angular radius of the source. We deduced θ_* from the color and brightness of the source. To estimate the extinction- and reddening-corrected (de-reddened) values, ${(V-I,I)}_0$, from the instrumental values, (V − I, I), we used the Yoo et al. (2004) method, in which the centroid of a red giant clump (RGC), ${(V-I,I)}_{\text{RGC}}$, in the color-magnitude diagram (CMD) is used for the calibration of the source color and magnitude.

Figure 10 shows the locations of the source (empty blue dot with error bars) and the RGC centroid (solid red dot) in the instrumental CMD constructed using the KMTC data processed using the pyDIA software (Albrow 2017). The source position in the CMD was determined via the regression of the I- and V -band pyDIA data with the variation of the lensing magnification. Also marked in the CMD is the blend position (solid green dot), which lies on the main-sequence branch of foreground disk stars. We show in Sect. 6 that the blended light is due to the host star, a companion to the host star, or a combination of the two. The measured instrumental color and brightness are (V − I, I) = (3.378 ± 0.099, 22.018 ± 0.004) for the source and ${(V-I,I)}_{\text{RGC}}=(3.606,16.856)$ for the RGC centroid. By measuring the offsets in color, Δ(V − I), and magnitude, ΔI, between thesource and RGC centroid, together with the known de-reddened values of the RGC centroid, ${(V-I,I)}_{\text{RGC},\text{0}}=(1.060,14.268)$, from Bensby et al. (2013) and Nataf et al. (2013), respectively, we estimate that the de-reddened source color and magnitude are $$\begin{array}{lll}{(V-I,I)}_0\hfill & =\hfill & {(V-I,I)}_{\text{RGC},0}+\Delta (V-I,I)\hfill \\ \hfill & =\hfill & (0.832\pm 0.099,19.431\pm 0.004).\hfill \end{array}$$(2)

This indicates that the source is a late G-type main-sequence star located in the bulge.

For the θ_* estimation from ${(V-I,I)}_0$, we first converted the measured V − I color into V − K color using the color-color relation of Bessell & Brett (1988) and second deduced the source radius from the (V − K)–θ_* relation of Kervella et al. (2004). The estimated value is $${\theta }_{*}=0.47\pm 0.06 \mu \text{as}.$$(3)

Together with the ρ and t_E values measured from the modeling, the angular Einstein radius and the relative lens-source proper motion are estimated as $${\theta }_{\text{E}}=0.91\pm 0.13 \text{mas}$$(4)

and $$\mu =4.41\pm 0.63 \text{mas} {\text{yr}}^{-1},$$(5)

respectively.

Fig. 10 Locations of the source (empty blue dot with error bars), RGC centroid (solid red dot), and blend (solid green dot) in the instrumental CMD of stars around the source constructed using the pyDIA photometry of the KMTC data set.

5 Physical lens parameters

Although the mass and distance to the lens cannot be uniquely determined using the relation in Eq. (1) due to the insecure measurement of the microlens parallax, the physical parameters can still be constrained using the measured observables of t_E and θ_E, which are related to M and D_L by $$t_{\text{E}}=\frac{{\theta }_{\text{E}}}{\mu };{\theta }_{\text{E}}={(\kappa M{\pi }_{\text{rel}})}^{1/2}.$$(6)

Here, ${\pi }_{\text{rel}}=\text{AU}(D_{\text{L}}^{-1}-D_{\text{S}}^{-1})$ represents the relative lens-source parallax and D_S denotes the distance to the source. For this constraint, we conducted a Bayesian analysis using a prior Galactic model. Although the uncertainty of the measured π_E is big, we considered the measured π_E by computing the parallax-ellipse covariance matrix in the Bayesian estimation of the physical lens parameters.

The Galactic model defines the physical and dynamical distributions and the mass function of Galactic objects. We adopted the Galacticmodel described in Jung et al. (2021). This model uses the Robin et al. (2003) and Han & Gould (2003) physical matter distributions to specify the locations of disk and bulge objects, respectively. It also uses the Han & Gould (1995) and Jung et al. (2021) dynamical distributions to describe the motions of disk and bulge objects, respectively. The mass function adopted in the model is described in Jung et al. (2018), and it includes stellar remnants and brown dwarfs.

In the first step of the Bayesian analysis, we conducted a Monte Carlo simulation to produce many (10⁷) artificial lensing events, for which lens and source locations, transverse lens-source speeds, and lens masses are allocated following the Galactic model. Then, the microlensing observables (t_E, θ_E, π_E) of the individual events were computed using the relations in Eqs. (1) and (6). In the second step, we constructed the distributions of M and D_L of events for which their observables are consistent with the observed values.

Figure 11 shows the posterior distributions of M and D_L obtained from the Bayesian analysis. In Table 3, we list the estimated physical parameters of the masses of the individual lens components (M₁, M₂, M₃), the distance, and the projected physical separations (a_⊥,2, a_⊥,3) between M₁ –M₂ and M₁ –M₃ pairs for the close-close and close-wide solutions. We adopted the median values of the distributions as representative parameters, and the uncertainties of the parameters were estimated as 16 and 84% of the distributions. The two solutions result in similar physical parameters, except for a_⊥,3. The estimated mass of the host, $$M_1={0.92}_{-0.32}^{+0.49} M_{\odot },$$(7)

indicates that the host is a G-type star, and the distance, $$D_{\text{L}}={4.40}_{-1.02}^{+1.15} \text{kpc},$$(8)

indicates that the lens is likely to be in the disk. The two planets have masses $$(M_2,M_3)=\{\begin{array}{ll}({3.43}_{-1.17}^{+1.83} M_{\text{J}},{10.22}_{-3.50}^{+5.46} M_{\text{J}}),\hfill & \text{for}\text{close-close},\hfill \\ ({3.21}_{-1.10}^{+1.71} M_{\text{J}},{10.14}_{-3.47}^{+5.42} M_{\text{J}}),\hfill & \text{for}\text{close-wide},\hfill \end{array}$$(9)

and thus both planets are heavier than Jupiter. Considering that the snow line distance from the host is d_sl = 2.7(M∕M_⊙) ~ 2.5 AU, the estimated separations, $$(a_{\perp ,2},a_{\perp ,3})=\{\begin{array}{ll}({3.29}_{-2.52}^{+4.15} \text{AU},{2.77}_{-2.12}^{+3.49} \text{AU}),\hfill & \text{for}\text{close-close},\hfill \\ ({3.29}_{-2.53}^{+4.16} \text{AU},{5.32}_{-4.08}^{+6.71} \text{AU}),\hfill & \text{for}\text{close-wide},\hfill \end{array}$$(10)

indicate that one planet lies near the snow line, and the other planet lies beyond the snow line.

Fig. 11 Bayesian posterior distributions of the mass of (upper panel) and distance to (lower panel) the lens. In each panel, the solid vertical line indicates the median of the distribution, and the two dotted vertical lines indicate the 1σ range of the distribution. The blue and red curves represent the contributions by disk and bulge lenses, respectively.

6 Baseline object

From the position of the blend on the CMD (Fig. 10), it is found that the flux of the baseline object is dominated by the light from a main-sequence star in the foreground disk. Logically, there are only four possibilities for this star: (1) the host of the planets, (2) a companion to the host, (3) a companion to the source, or (4) an ambient star that is unrelated to the event (or possibly a combination of two or more of these possibilities). The position of the blend on the CMD already rules out possibility (3) because it is inconsistent with a bulge star.

Because θ_E ~ 0.9 mas, the lens must also be in the foreground disk (unless it is a dark remnant). That is, if the lens were in the bulge, with π_rel ≲0.02 mas, then its mass would be $M={\theta }_{\text{E}}^2/\kappa {\pi }_{\text{rel}}\gtrsim 5M_{\odot }$, implying that it would be easily seen. This suggests that the blended light may be due primarily to the host of the planets.

To test this conjecture, we first measured the astrometric offset between the “baseline object” and the source (which has the same position as the host at the time of the event). The position of the source is measured from difference images formed by subtracting the reference image from a series of images taken at high magnifications. This difference image essentially consists of an isolated star on a blank background. Hence, its position can be accurately measured. We find a scatter among 15 images of just (0.026, 0.037) pixels in the (east, north) directions, yielding a standard error of the mean of (0.007, 0.010) pixels, that is, (3, 4) mas, given the 400 mas pixel size. These errors are far below the other errors in the problem, which are discussed below. We therefore ignore them in what follows.

By contrast, the problem of measuring the position of the baseline object is far subtler. Its position is returned by the DoPhot photometry package (Schechter et al. 1993), which simultaneously fits for the positions and fluxes from possibly overlapping stellar images. Fortunately, the baseline object appears isolated in these images, so the statistical errors of this measurement are not strongly impacted by the algorithm’s procedure for separating stars. We estimated this error to be 23 mas in each direction from a comparison of two completely independent reductions, that is, by reductions determined by different people. We find the offset between the source and the baseline object to be $$\Delta \theta (E,N)=(36\pm 23,73\pm 23) \text{mas}.$$(11)

If we could ignore systematic effects and take the error distributions to be a two-dimensional Gaussian, this would rule out the identification of the host with baseline light at p = 0.002 probability.

We note that this close alignment rules out possibility (4), that is, that the blended light is due to an ambient star. We simply counted the number of stars on the foreground main sequence in the 2 arcmin² area of the OGLE finding chart that are brighter than the baseline object, finding N = 64. This translates to a surface density of n = 64∕120² = 4.4 × 10⁻³ arcsec⁻². Thus, the probability of such an alignment is $p=\pi {(\Delta \theta )}^{2}n={10}^{-4}$. Hence, the light from the baseline object must be due to either the planet host or a companion to the host (or a combination of the two).

Before considering these two remaining possibilities, we examined the role of the key systematic effect that could corrupt the position measurement of the baseline object at the $\Delta \theta ={[\Delta \theta {(N)}^2+\Delta \theta {(E)}^2]}^{1/2}~80$ mas level: namely, the possibility that an ambient star lies within the point-spread function of the baseline object. That is, if a fraction, f, of the blended light were due to an ambient star with a separation from the source of δθ, then the centroid of light (measured by DoPhot) would be displaced by Δθ = fδθ. For instance, for f = 0.2 and δθ = 400 mas, Δθ = 80 mas. In this example, the faint ambient starwith a separation of 0.4′′ would not be identified as a separate star by DoPhot, so the centroid would be shifted.

We have just argued that ambient stars similar to or brighter than the baseline object are rare. However, ambient stars that could corrupt the astrometric measurement at this level are not rare, for three reasons. First, such stars are not restricted to the foreground main sequence, so the much larger population of bulge stars is available. Second, there are more faint stars than bright stars. Third, ambient stars can corrupt the astrometric measurement with separations of up to δθ ~ 1′′, whereas ambient stars that would explain the baseline object light must lie within Δθ = 80 mas. These three features also place limits on which ambient stars can play this role. First, if the ambient star were too bright (and came from the more populous, but redder, bulge population), then it would also drive the combined color of the resulting baseline object to the red, and so off the foreground main-sequence feature, contrary to its actual location. Second, if the ambient star were too faint, its separation would have to be δθ≳1′′, at which point its light would no longer enter into the DoPhot centroid. Third, as δθ is increased,ambient stars of sufficient brightness will be separately resolved by DoPhot before the δθ = 1′′ limit is reached.

To make our evaluation, we first focused our attention on the 2349 stars in the 2 arcmin² OGLE-IV finding chart that satisfy 0.7 < (I − I_base) < 2.7, where I_base = 17.8 is evaluatedin the OGLE-IV system. The faint limit is set by the requirement δθ < 1′′, while the bright limit is set by the requirement that inclusion of the (usually) red ambient light does not drive the baseline object off the foreground main sequence. Then, for each of these 2349 stars we found the annulus of positions such that the star would induce an astrometric error of > 80 mas. The inner radius of the annulus is set by the flux ratio: δθ_inner = Δθ × 10^{0.4(I−I_base)}. We set δθ_outer = 0.2′′ [(I − I_base) + 2.3]. We then find a total area subtended by these 2349 annuli to be 2800 arcsec² (i.e., a fraction p = 19% of the 2 arcmin²) That is, if the host were primarily responsible for the blended light, then there is a p = 19% probability that an ambient star would corrupt the astrometric measurement by enough to account for the observed Δθ ~ 80 mas offset.

On the other hand, it is also possible that the baseline object has a fainter companion that generated the main microlensing event and so serves as the host for the planets. The main constraint on this scenario is that the blend, which is likely a G dwarf at several kiloparsecs, should have a widely separated companion at log (P∕day) ~ 4, with a mass ratio of q_base≳0.5. According to Table 7 of Duquennoy & Mayor (1991), about 25% of G dwarfs have companions with log (P∕day) > 4 and a mass ratio q ≥ 0.5. However, if the baseline object did have a companion in this parameter range, it is about equally likely that the baseline object served as the lens host for the event while the companion generated the flux required to corrupt the astrometric measurement.

In brief, there are two channels for the baseline object to be the host, with the astrometric measurement corrupted either by an ambient star (p = 0.19) or by a less luminous companion widely separated from the planet host (p = 0.25∕2 ~ 0.12) for a total of p = 0.31. By contrast, there is one channel for the host to be a companion of the baseline object, with probability p = 0.25∕2 ~ 0.12. Therefore, the host of the planet is likely responsible for the light of the baseline object, but the host could also plausibly be a companion to the baseline object. Considering the inconclusive nature of the blend together with the low quality of the spectra, we do not impose the constraint from the spectra on the estimation of the physical parameters presented in Table 3.

7 Summary and conclusion

We have reported a multiple planetary system discovered from the analysis of the lensing event OGLE-2019-BLG-0468, which was reinvestigated as part of a project of reviewing the microlensing data collected in and before the 2019 season by the KMTNet survey. The light curve of the event, which consists of three distinctive anomaly features, could not be explained by the usual 2L1S or 1L2S interpretations. We found a solution explaining all anomaly features with a 3L1S interpretation, in which the lens is composed of two planets and their host, making the lens the fourth multi-planet system securely found by microlensing. The lensing solution is subject to two-fold degeneracies caused by the ambiguity in estimating the separations of the planets from the host. One of the degeneracies is very severe, but the other was resolvable due to the resonant nature of the caustic induced by the second planet. The two planets have masses of ~ 3.4 M_J and ~ 10.2 M_J, and they are orbiting around a G-type star with a mass of ~0.9 M_⊙ and a distance of ~4.4 kpc. It was found that the planet host was most likely responsible for the light of the baseline object, although the possibility of the host being acompanion to the baseline object could not be ruled out.

Acknowledgements

Work by C.H. was supported by the grants of National Research Foundation of Korea (2020R1A4A2002885 and 2019R1A2C2085965). S.D. acknowledges the science research grants from the China Manned Space Project with NO. CMS-CSST-2021-A11. 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 observations using the SALT telescope were conducted under the transients follow up programme 2018-2-LSP-001 (PI: D.B.), which is supported by Poland under grant no. MNiSW DIR/WK/2016/07. D.B. also acknowledges research support from the National Research Foundation. M.G. is supported by the EU Horizon 2020 research and innovation programme under grant agreement No. 101004719.

References

All Tables

All Figures

	Fig. 1 Light curve of the microlensing event OGLE-2019-BLG-0468. The arrows labeled t1 (8586.3), t2 (8592.9), and t3 (8595.7) indicate the three epochs of major anomalies. The curves drawn over the data points represent the 1L1S (dotted curve) and 3L1S (close-close model, solid curve) models. The zoomed-in view of the anomaly region around t2 and t3 is shown in Fig. 2.
In the text

	Fig. 2 Zoomed-in view of the anomaly region around the epochs t2 and t3. Notations are the same as in Fig. 1.
In the text

	Fig. 3 Comparison of the four tested models, 1L2S, 2L1S, 2S2L, and 3L1S (close-close solution). The lower panels show the residuals from the individual models.
In the text

	Fig. 4 Distribution of Δχ2 in the logs3– logq3 parameter plane obtained from the preliminary 3L1S grid searches for the parameters related to the third lens component. Red, yellow, green, cyan, and blue are used to denote points with Δχ2 < n(12), < n(22), < n(32), < n(42), and < n(52), respectively, where n = 10.
In the text

Fig. 5 Lens system configurations of the close-close (left panel) and close-wide (right panel) 3L1L models. The inset in each panel shows the locations of the three lens components, denoted by blue dots labeled M1, M2, and M3. The dashed circle in the inset represents the Einstein ring. The line with an arrow represents the source trajectory. The three small empty dots on the source trajectory represent the source positions at the epochs of t1, t2, and t3 that are marked in Fig. 1. The size of the dots is not scaled to the source size. The red figure represents the caustic, and the gray curves around the caustic represent equi-magnification contours.

In the text

Fig. 6 Magnification pattern around the caustics (for the close-wide 3L1S model) and the light curve around the time of the source approach close to the planetary caustic induced by M2. The small empty circle on the source trajectory represents the source position at the time of the source approach, at tapproach ~ 8581.7 in HJD′. The size of the circle is arbitrarily set and is not scaled to the source size.

In the text

	Fig. 7 Comparison of the residuals from the four degenerate 3L1S solutions: close-close, close-wide, wide-close, and wide-wide solutions. As a representative model, the model curve of the close-close solution is drawn over the data points in the top panel.
In the text

	Fig. 8 Cumulative distributions of the Δχ2 difference from the best-fit solution (close-close model). The light curve in the upper panel is inserted to show the locations of the fit difference.
In the text

	Fig. 9 Scatter plot of points in the MCMC chain on the πE,E–πE,N parameter plane for the close-close 3L1S solution. The color coding is the same as in Fig. 4, except n = 1.
In the text

	Fig. 10 Locations of the source (empty blue dot with error bars), RGC centroid (solid red dot), and blend (solid green dot) in the instrumental CMD of stars around the source constructed using the pyDIA photometry of the KMTC data set.
In the text

	Fig. 11 Bayesian posterior distributions of the mass of (upper panel) and distance to (lower panel) the lens. In each panel, the solid vertical line indicates the median of the distribution, and the two dotted vertical lines indicate the 1σ range of the distribution. The blue and red curves represent the contributions by disk and bulge lenses, respectively.
In the text