OGLE-2019-BLG-0304: Competing Interpretations between a Planet–binary Model and a Binary-source + Binary-lens Model

The search for planets belonging to binary and multiple stellar systems is important because these planets are expected to have gone through different formation and evolution processes than those of single stars, and thus they can provide a valuable test bed to better understand the formation and evolution processes of planets in general. The search for these planets is also important in estimating the global frequency of planets because a majority of stars form binary or multiple systems, and thus the planetary frequency of multiple systems can have a significant effect on the global planet frequency.

Planets in binary and multiple systems have been detected using various methods, including radial velocity (Correia et al. 2005), transit (Doyle et al. 2011), pulsar-timing (Thorsett et al. 1993), eclipsing-binary (Lee et al. 2009), and microlensing (Gould et al. 2014) methods. Among these methods, microlensing is useful for detecting some specific populations of planets, especially cold planets orbiting low-mass binary stars. The microlensing sensitivity to planets in low-mass binaries is due to the lensing property, in which detection does not depend on the luminosity of host stars. Another important advantage of the microlensing method is that it enables one to detect planets in both S- and P-type orbits, for which a planet of the former type orbits one of the the widely separated binary stars, and a planet of the latter type orbits both closely spaced binary stars. The microlensing sensitivity to both S- and P-type planets is due to the lensing characteristics, in which a binary companion, regardless of binary separation, induces a caustic in the central magnification region, near which a planet also induces a caustic, and thus the signatures of both the planet and stellar companion appear in the same peak region of high-magnification lensing light curves (Han 2008; Lee et al. 2008). Here caustics indicate the positions on the microlensing source plane at which the lensing magnification of a point source would become infinite. For this reason, high-magnification lensing events provide a major channel (central perturbation channel) for detecting planets in binaries. ¹⁴

Although high-magnification lensing events provide a channel to detect both S- and P-type planets, one is often confronted with cases in which the orbital type of the planet, S or P, cannot be uniquely determined. This happens because both the close- and wide-separation stellar companions can induce similar caustics in the central magnification region, producing companion signals of a similar shape. This so-called "close–wide degeneracy" (Griest & Safizadeh 1998; Dominik 1999) causes ambiguity not only in determining the binary separation but also in determining the star–planet separation. The close–wide degeneracy happened in three (OGLE-2006-BLG-284, OGLE-2016-BLG-0613L, and OGLE-2018-BLG-1700L) out of the total six known microlensing planetary systems in binaries, thus it is a deep-seated problem in determining the orbital type of planetary systems in binaries.

The close–wide degeneracy in determining the planet and binary separations can be lifted in some special lens-system configurations. The first such case is that the projected separation between the planet and host is similar to the angular Einstein radius θ_E of the lens system. In this case, the planet induces a resonant caustic, for which the planet-host separation is uniquely determined. The second case is that the source trajectory of a lensing event passes the region around a widely separated second stellar lens component well after (or well before) the main approach to the first lens component. In this case, the source approach to the second lens component produces an extra bump in addition to the main peak produced by the source approach close to the first lens component (Di Stefano & Scalzo 1999). Then, the existence of the second bump indicates that the planet is in an S-type orbit.

In this work, we present the analysis of the lensing event OGLE-2019-BLG-0304. The light curve of the event exhibits two distinctive features: a deviation in the peak region and a bump appearing long after the main peak. A 2L1S model approximately describes the overall feature of the light curve, but it leaves subtle residuals in the peak region. If the central residual is caused by the perturbation of a planetary companion, then the lens is a special case in which a planet belongs to a binary and its orbital type is identified by the second bump. In order to check this possibility, we conduct thorough model testing under various interpretations of the lens system and present the results.

For the presentation of the work, we organize the paper as follows. In Section 2, we describe the observations of the lensing event and the data obtained from the observations. We also mention the characteristics of the lensing event. In Section 3, we describe various lensing models that are tested to interpret the observed lensing light curve. The procedure of estimating the angular Einstein radius and relative lens-source proper motion is stated in Section 4. In Section 5, we describe the procedure of estimating the physical lens parameters using the measured observables of the lensing event. We discuss the reality of the signal in Section 6, and conclude in Section 7.

The source of the event OGLE-2019-BLG-0304 is located toward the Galactic bulge field at the equatorial coordinates (R.A., decl.)_J2000 = (17:36:06.41, −26:08:45.56). The corresponding Galactic coordinates are (l, b) = (1252, 3266). The brightness of the source had remained constant before the 2018 season with a baseline magnitude of I_base = 17.99. The lensing-induced brightening of the source flux started during the time between the end of the 2018 season and the beginning of the 2019 season. The lensing light curve peaked on 2019 March 1 (HJD' ≡ HJD−2,450,000 ∼ 8543.5) with a magnification of A_peak ∼ 60, and then gradually declined to the baseline. The source flux increased and peaked again on 2020 May 1 ($\mathrm{HJD}^{\prime} \sim 8605$), producing a second bump with a low magnification of A_bump ∼ 1.2.

The lensing event was first found by the Early Warning System of the Optical Gravitational Microlensing Experiment (OGLE, Udalski et al. 1994) survey on 2019 March 16. The OGLE survey utilized the 1.3 m telescope located at the Las Campanas Observatory in Chile and the telescope is equipped with a camera yielding a 1.4 deg² field of view. The event was rediscovered by the Korea Microlensing Telescope Network (KMTNet, Kim et al. 2016) survey in its post-season analysis (Kim et al. 2018), and it was dubbed KMT-2019-BLG-2583. The KMTNet survey employed three identical 1.6 m telescopes located at the Siding Spring Observatory (KMTA) in Australia, the Cerro Tololo Inter-American Observatory (KMTC) in South America, and the South African Astronomical Observatory (KMTS) in Africa. The camera mounted on each KMTNet telescope yields a 4 deg² field of view.

For both surveys, observations were conducted mainly in the I band and a fraction of images were taken in the V band for the source color measurement. We describe the detailed procedure of the source color measurement in Section 4. The OGLE survey of the 2019 season started on $\mathrm{HJD}^{\prime} =8547$, 4 days before the peak, but the peak was not covered because the source of the event lies in a low-cadence field (BLG667), which is very infrequently observed at the beginning of the season due to short observing nights. Fortunately, the peak was covered by the KMTNet survey, which commenced in the 2019 season on $\mathrm{HJD}^{\prime} =8534$, 8518, 8534 for the KMTA, KMTC, and KMTS telescopes, respectively.

Neither the OGLE nor the KMT survey issued an alert of the event prior to peak, and hence no follow-up observations were possible, although this was a high-magnification event, for which the sensitivity to planets is high. For example, follow-up observations of OGLE-2012-BLG-0026 (Han et al. 2013) revealed a two-planet system, despite the fact that it peaked at a similar calendar date to OGLE-2019-BLG-0304, so that observations were restricted by similarly short nights. The delay of the OGLE alert was caused by the sparse observations of this field, as mentioned above. KMT did not issue an alert at all because its AlertFinder (Kim et al. 2018) system began operation in 2019 on March 27, and it only searches for rising events. Hence, as mentioned above, the event was found in the post-season analysis.

Reduction and photometry of the data were done using the pipelines of the individual survey groups developed based on the difference imaging method (Tomaney & Crotts 1996; Alard & Lupton 1998): the Udalski (2003) code using the DIA technique of Woźniak (2000) for the OGLE survey and the Albrow et al. (2009) code for the KMTNet survey. An additional reduction was conducted for a subset of the KMTC data using the pyDIA code (Albrow 2017) to estimate the color of the source star. Error bars of the data estimated by the photometry pipelines were readjusted following the routine described in Yee et al. (2012).

Figure 1 shows the light curve of OGLE-2019-BLG-0304 with the combined data obtained from the OGLE and KMTNet surveys. The inset in the upper panel shows a zoom-in view of the peak region. Compared to the symmetric light curve of a single-lens, single-source (1L1S) event, the light curve of the event exhibits two distinctive features. First, the light curve is asymmetric due to the existence of the second bump, which is centered at $\mathrm{HJD}^{\prime} \sim 8605$. Second, the region around the main peak, centered at $\mathrm{HJD}^{\prime} \sim 8543.5$, appears to exhibit deviations caused by finite-source effects, which occur when the lens passes over the surface of the source star. The curve drawn over the data points is the 1L1S model obtained by fitting the observed data excluding the region of the second bump with the consideration of finite-source effects. The model parameters are (t₀, u₀, t_E, ρ) ∼ (8543.68, 0.04 × 10⁻³, 15.3 days, 4.1 × 10⁻²), where the individual parameters indicate the peak time (in $\mathrm{HJD}^{\prime}$), lens-source separation at that time (normalized to θ_E), event timescale, and source radius (also normalized to θ_E). Although the finite-source 1L1S model better describes the peak region than the point-source model does, it still leaves substantial residuals, indicating that a more sophisticated model is needed to explain not only the second bump but also the main peak.

Zoom In Zoom Out Reset image size

3.1. Three-body Interpretations

Considering the existence of the second bump, we first test two three-body (lens plus source) interpretations, in which one is that the source is a binary (1L2S model) and the other is that the lens is a binary (2L1S model). Besides the 1L1S parameters of (t₀, u₀, t_E, ρ), modeling the light curve under these interpretations requires one to include additional parameters. For the 1L2S model, these additional parameters are (t_0,2, u_0,2, q_F , ρ₂), which represent the time of the second bump, separation between the lens and the second source at t_0,2, the flux ratio between the two source stars (S₁ and S₂), and the normalized source radius of the second source, respectively (Hwang et al. 2013). For the 2L1S model, the additional parameters are (s, q, α), which represent the separation (scaled to θ_E) and mass ratio between the two lens components (M₁ and M₂) and the angle of the source trajectory as measured from the M₁–M₂ axis (source trajectory angle), respectively.

3.1.1. Binary-source (1L2S) Model

We start with the modeling under the 1L2S interpretation. In this modeling, the initial parameters of (t₀, u₀, t_E, ρ) are set as those obtained from the 1L1S modeling, and the initial parameters of (t_0,2, u_0,2, q_F , ρ₂) are set considering the time and height of the second bump in the light curve. From this modeling, we find a solution that can explain the second bump. The lensing parameters of the 1L2S model along with the χ² value of the fit are presented in Table 1. We note that ρ₂ is not presented because the second bump does not exhibit deviations caused by finite-source effects.

Table 1. Lensing Parameters of 1L2S, 2L1S, 2L2S, and 3L1S Models

Parameter	1L2S	2L1S	2L2S	3L1S
χ2	7359.5	5654.6	5614.3	5606.3
t0 ($\mathrm{HJD}^{\prime}$)	8543.677 ± 0.00266	8543.610 ± 0.004	8543.474 ± 0.030	8543.599 ± 0.005
u0 (10−3)	−0.24 ± 0.94	0.86 ± 0.27	0.07 ± 1.77	0.54 ± 0.36
t0,2 ($\mathrm{HJD}^{\prime}$)	8604.245 ± 0.148	⋯	8543.719 ± 0.037	⋯
u0,2 (10−3)	0.831 ± 0.026	⋯	1.80 ± 0.99	⋯
tE (days)	14.57 ± 0.10	17.47 ± 0.08	17.48 ± 0.09	17.76 ± 0.11
s2	⋯	3.784 ± 0.015	3.785 ± 0.016	3.733 ± 0.019
q2	⋯	0.391 ± 0.010	0.378 ± 0.012	0.363 ± 0.012
α (rad)	⋯	2.954 ± 0.002	2.958 ± 0.003	2.960 ± 0.003
s3	⋯	⋯		0.885 ± 0.005
q3 (10−3)	⋯	⋯		1.82 ± 0.26
ψ (rad)	⋯	⋯		2.403 ± 0.109
ρ1 (10−2)	4.37 ± 0.06	1.97 ± 0.03	1.83 ± 0.07	1.85 ± 0.04
ρ2 (10−2)	⋯	⋯	2.02 ± 0.09	⋯
qF	0.30 ± 0.02	⋯	1.277 ± 0.524	⋯

Note. HJD' ≡ HJD–2,450,000.

Download table as:  ASCIITypeset image

From inspecting the residuals of the model, presented in the bottom panels of Figures 2 (around the main peak) and 3 (around the second bump), it is found that the model still leaves substantial residuals, especially around the main peak. This indicates that the 1L2S model is not a correct interpretation of the observed light curve.

Zoom In Zoom Out Reset image size

3.1.2. Binary-lens (2L1S) Model

The 2L1S modeling is carried out in two steps. In the first step, the binary parameters s and q are searched for using a grid approach and the remaining parameters are searched for using a downhill approach based on the Markov Chain Monte Carlo algorithm. In the second step, the local solutions found from the first step are refined by allowing all parameters, including s and q, to vary. This two-step process is needed to investigate the existence of possible degenerate solutions.

From the 2L1S modeling, we find a unique model that can explain the overall features of the light curve by significantly reducing the residuals in both regions of the main peak and the second bump of the light curve. In Figures 2 and 3, we present the model curve and residuals of the 2L1S solution in the regions around the main peak and second bump, respectively. The best-fit lensing parameters of the solution are listed Table 1, in which we mark the binary parameters with a subscript "2," i.e., (s₂, q₂), to distinguish them from the parameters related to a possible tertiary lens component to be discussed in the following subsection. In Figure 4, we present the lens-system configuration of the 2L1S solution, showing the source trajectory with respect to the locations of the lens components and the resulting caustics. According to the model, the event is produced by a binary lens, in which the companion M₂ is separated in projection from the primary M₁ by s ∼ 3.8, and its mass ratio to the primary is q = M₂/M₁ ∼ 0.39. The anomaly around the main peak of the light curve is produced by the source star crossing over the small caustic located close to the primary of the lens and the second bump is explained by the approach of the source near the region around the second caustic located close to the companion. We note that despite the source star crossing over the caustic during the main peak, the usual sharp caustic crossing features do not appear in the light curve due to severe finite-source effects. It is found that the caustic crossing interpretation greatly reduces the residuals from the 1L1S model in the peak region of the light curve. This, together with the second bump, strongly indicates that the lens is accompanied by a binary companion.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

We note that the 2L1S solution would have been subject to the close–wide degeneracy, if it were not for the data around the second bump. To check this, we conduct an additional modeling by excluding the data lying in the region $8570\lt \mathrm{HJD}^{\prime} \lt 8635$. Figure 5 shows the Δχ² map in the s₂–q₂ parameter plane obtained from this modeling. The map shows two distinct locals, in which one with s₂ > 1.0 corresponds to the wide solution presented in Table 1, and the other with s₂ < 1.0 is the degenerate close solution. Therefore, the event is a special case, in which the close–wide degeneracy is clearly resolved by the existence of a second bump.

Zoom In Zoom Out Reset image size

Although the 2L1S model appears to depict the overall feature of the light curve, we find that it leaves subtle but noticeable residuals in the peak region. This can be seen in the third residual panel (labeled as "2L1S") of Figure 2, which shows that the data points in this region deviate from the model by ≲0.05 mag. The three dotted vertical lines (marked by t₁ = 8542.20, t₂ = 8542.85, and t₃ = 8543.85) indicate the three epochs of relatively large deviations. In order to check whether the deviations are caused by systematics in the data, we conduct multiple re-reductions of the data using different template images for difference imaging photometry. We find that the residuals persist regardless of the reduction, suggesting that the signal is real.

We further check the possibility that the residuals from the 2L1S model around the main peak are caused by the omission of higher-order effects in modeling. For this check, we conduct an additional modeling considering the microlens parallax (Gould 1992) and lens orbital (Dominik 1998) effects. The microlens-parallax effects are caused by the positional change of an observer due to the orbital motion of Earth around the Sun, and the lens orbital effects are caused by the change of the lens position due to the orbital motion of a binary lens. The modeling considering the microlens-parallax effect requires one to include two additional parameters of π_E,N and π_E,E , which denote the north and east components of the microlens parallax vector π _E = (π_rel/θ_E)( μ /μ), respectively. Here ${\pi }_{\mathrm{rel}}=\mathrm{au}({D}_{{\rm{L}}}^{-1}-{D}_{{\rm{S}}}^{-1})$ is the relative lens-source parallax, μ denotes the relative lens-source proper motion vector, and (D_L, D_S) are the distances to the lens and source, respectively. Considering the lens orbital effect also requires one to add two parameters of ds/dt and d α/dt, which represent the change rates of s and α, respectively. From this modeling, we find that the residuals from the model near the main peak of the light curve still persist, indicating that the cause of the residuals is not the higher-order effects.

3.2. Four-body Interpretations

To explain the residuals from the 2L1S model, we test two additional four-body interpretations, in which one is that both the lens and source are binaries (2L2S model) and the other is that the lens is a triple system (3L1S model).

3.2.1. Binary-source + Binary-lens (2L2S) Model

We conduct a 2L2S modeling considering the possibility that the discontinuous anomaly features at around t₁, t₂, and t₃ are caused by a companion to the source. The introduction of an additional source to a 2L1S model requires one to add extra parameters in modeling, including t_0,2, u_0,2, and ρ₂. Here ρ₂ denotes the normalized source radius of the companion (S₂) to the primary source star (S₁). In the modeling, we use the initial parameters related to S₁ as the ones obtained from the 2L1S modeling because the 2L1S model describes the overall feature of the observed data. The initial parameters related to S₂ (t_0,2, u_0,2, ρ₂, and q_F ) are chosen considering the times and magnitudes of the deviations from the 2L1S model.

We plot the model curve of the 2L2S solution and residuals from the model in Figures 2 and 3. The lensing parameters of the model are presented in Table 1 and the corresponding lens-system configuration is shown in Figure 6. The configuration is very similar to that of the 2L1S solution, except that there is an additional trajectory for S₂. According to the 2L2S solution, S₂, which is brighter than S₁ by ∼28% in the I-band flux, trails S₁ with a projected separation of ${\rm{\Delta }}u=\{{[{({t}_{\mathrm{0,2}}-{t}_{\mathrm{0,1}})/{t}_{{\rm{E}}}]}^{2}+({u}_{\mathrm{0,2}}-{u}_{\mathrm{0,1}})}^{2}\}=0.014$, and crosses the caustic. We mark the positions of S₁ and S₂ corresponding to the times of t₁, t₂, and t₃ as orange and brown circles, respectively.

Zoom In Zoom Out Reset image size

It is found that the model fit substantially improves with respect to the 2L1S model by introducing an additional source component. First, the anomaly feature at around t₃, which exhibits the largest deviation from the 2L1S model, is explained by the caustic crossing of S₂. Second, the anomaly feature at around t₂ is explained by the passage of S₂ through the positive deviation region extending from the left-side cusp of the caustic lying on the M₁–M₂ axis. However, the model fit around the anomaly at t₁ remains almost unchained. In combination, the 2L2S model results in a better fit than the 2L1S model by ${\rm{\Delta }}{\chi }^{2}={\chi }_{2{\rm{L}}1{\rm{S}}}^{2}-{\chi }_{2{\rm{L}}2{\rm{S}}}^{2}=40.3$. We reserve making conclusions about this model until we check an additional model to be discussed below.

3.2.2. Triple-lens (3L1S) Model

We check a 3L1S interpretation because the deviations from the 2L1S model appear in the peak region, at which an additional anomaly would occur if the lens has a tertiary component M₃. The 3L1S modeling is carried out under the assumption that the magnification pattern of a triple-lens system can be approximated by the superposition of the patterns produced by two (M₁–M₂ and M₁–M₃) binary pairs (Bozza 1999; Han et al. 2001). Under this approximation, we search for the parameters related to M₃ by fixing the lensing parameters related to M₂ as those of the 2L1S model. The parameters related to M₃ are (s₃, q₃, ψ), which represent the projected separation and mass ratio between M₁ and M₃, and the position angle of M₃ as measured from the M₁–M₂ axis, respectively (Han et al. 2013).

Under the superposition approximation, we first conduct thorough grid searches for (s₃, q₃, ψ). Figure 7 shows the Δχ² map on the s₃–q₃ plane obtained from the grid search. It shows a unique and distinct local at (s₃, q₃) ∼ (0.89, 1.8 × 10⁻³). We refine the solution by allowing all parameters to vary. In Table 1, we list the best-fit lensing parameters of the 3L1S model. The estimated mass ratio between M₁ and M₂ of q₂ = M₂/M₁ ∼ 0.36 is similar to the ratio estimated from the 2L1S modeling. On the other hand, the mass ratio between M₁ and M₃, q₃ = M₃/M₁ = (1.82 ± 0.26) × 10⁻³, is very low, indicating that the tertiary lens component is a planetary mass object according to the 3L1S interpretation of the event. The projected separation between M₁ and M₃, s₃ ∼ 0.89, is substantially smaller than the separation between M₁ and M₂, s₂ ∼ 3.7. This indicates that the planet is in an S-type orbit, in which the planet is orbiting the heavier member of a widely separated binary. From the additional modeling considering higher-order effects, it is found that determining the higher-order lensing parameters is difficult mainly due to the short timescale, t_E ∼ 18 days, of the event.

Zoom In Zoom Out Reset image size

Figure 8 shows the lens-system configuration of the 3L1S model. The main panel shows the enlarged view of the central magnification region and the inset shows a wide view including the positions of the individual lens components. The planet is located on the M₁ side with a positional angle of ψ ∼ 138° as measured from the M₁–M₂ axis centered at M₁ in a counterclockwise direction. From the comparison of the 2L1S caustic shown in Figure 4, it is found that the tertiary lens component M₃ induces an additional caustic in the central region. The additional caustic appears to be a resonant caustic induced by a planetary companion located at a projected separation that is very nearly equal ¹⁵ to the Einstein radius. The resulting central caustic appears to be the superposition of a caustic induced by a binary companion and another caustic induced by a planet companion, although some deviations from the superposition exist because of the twisting and intersection of the caustics caused by the interference between the two (Gaudi et al. 1998; Rhie 2002; Daněk & Heyrovský 2015, 2019).

Zoom In Zoom Out Reset image size

We note that the separation of the planet from the host, s₃, is uniquely determined without any ambiguity. In general, the central caustics induced by a pair of planetary companions with separations s and s⁻¹ appear to be similar to each other due to the invariance of the binary-lens equation with the inversion of s and s⁻¹. However, this invariance breaks down as the planet separation is similar to θ_E, i.e., s ∼ 1.0, and the lens system forms a single large resonant caustic (Bozza 1999; An 2005; Chung et al. 2005). For OGLE-2019-BLG-0304, the planet induces such a resonant caustic, and thus determining s₃ does not suffer from the close–wide degeneracy.

It is found that the 3L1S model explains all major anomaly features. To better show the region of fit improvement around the main peak, we present the cumulative function of ${\rm{\Delta }}{\chi }_{3{\rm{L}}1{\rm{S}}}^{2}={\chi }_{2{\rm{L}}1{\rm{S}}}^{2}-{\chi }_{3{\rm{L}}1{\rm{S}}}^{2}$ in Figure 9. For the comparison of the model fit with that of the 2L2S model, we also present the distribution of ${\rm{\Delta }}{\chi }_{2{\rm{L}}2{\rm{S}}}^{2}={\chi }_{2{\rm{L}}1{\rm{S}}}^{2}-{\chi }_{2{\rm{L}}2{\rm{S}}}^{2}$. The fit improvement with respect to the 2L1S model can also be seen by comparing the residuals of the models presented in Figure 2. Although the fit improves throughout the rising part of the light curve during $8537\lesssim \mathrm{HJD}^{\prime} \lesssim 8544$, a major improvement occurs at the three epochs of t₁, t₂, and t₃. To be noted is that the 3L1S model explains the anomaly at t₁ which could not be explained by the 2L2S model. As a result, the fit of the 3L1S model is better than that of the 2L2S model by Δχ² = 8.0. The fit improvement over the 2L1S model is Δχ² = 48.3. We mark the source positions at the three epochs of the major fit improvement, i.e., t₁, t₂, and t₃, as orange circles in Figure 8. From this, it is found that the epoch t₁ corresponds to the time at which the source passes the cusp of the planet-induced caustic and the other two epochs correspond to the times at which the source passes over the folds of the planet-induced caustic.

Zoom In Zoom Out Reset image size

It is found that the second bump in the lensing light curve is important for the clear detection of the planetary signal. We find this fact by comparing the fits of the two sets of 2L1S–3L1S solutions, in which one set of solutions is obtained from modeling with the use of all data, and the other set of solutions is obtained by conducting modeling with the exclusion of the data around the region of the second bump. The cumulative Δχ² distribution obtained with the exclusion of the bump data is presented as a dashed curve (with a label "w/o bump" in the legend) in Figure 9. The χ² difference between the 3L1S and 2L1S models without the bump data is Δχ² = 29.4, which is substantially smaller than Δχ² = 48.3 when the bump data are included.

For the origin of the anomalies from the 2L1S model, the explanation with the existence of a low-mass tertiary lens component (3L1S model) is more plausible than the explanation with the existence of a very close companion to the source (2L2S model) for several reasons. First, the fit of the 3L1S model is better than the fit of the 2L2S model. Although Δχ² = 8.0 is not very big, the 3L1S model explains all three major anomaly features (at t₁, t₂, and t₃), while the 2L2S model cannot describe one (t₁) of these features. Second, the signal of a planet according to the 3L1S model shows up in the expected region of an event, for which the chance to detect such a signal is high, that is, around the peak of a high-magnification event. Third, indirect evidence comes from the fact that the projected separation between S₁ and S₂ according to the 2L2S model is too close for the binary system to be physically stable. The dereddened colors and magnitudes of the source stars are ${(V-I,I)}_{0,{S}_{1}}=(0.888\pm 0.053,15.383\pm 0.026)$ and ${(V-I,I)}_{0,{S}_{2}}=(1.007\pm 0.049,15.141\pm 0.023)$ for S₁ and S₂, respectively. In Figure 10, we mark the positions of S₁ and S₂ in the color–magnitude diagram (CMD). These indicate that the source stars are G- and K-type giants with radii of ${R}_{{S}_{1}}\sim 6\,{R}_{\odot }$ and ${R}_{{S}_{2}}\sim 8\,{R}_{\odot }$, respectively. On the other hand, the physical projected separation between the source stars is a_S,⊥ = ΔuD_S θ_E ∼ 4.3 R_⊙, which is significantly smaller than ${R}_{{S}_{1}}+{R}_{{S}_{2}}\sim 14\,{R}_{\odot }$. Here Δu = 0.014 is the S₁–S₂ separation in units of θ_E, and θ_E is estimated based on the lensing parameters of the 2L2S solution. In principle, the source stars could avoid merging if S₂ were projected in front of or behind S₁. However, unless this difference were considerable, the ellipsoidal distortions generated by the mutual tides of the two stars would induce ellipsoidal variations in the baseline light curve, which are not seen. Fourth, in order for two source stars to have evolved to the same brightness on the giant branch, they would have to have nearly identical mass, which is also less likely.

Zoom In Zoom Out Reset image size

However, the suggested reasons for the preference of the 3L1S model over the 2L2S model are either weak or indirect, and thus do not collectively provide sufficient evidence to strongly favor the 3L1S model. For the reason based on the χ² difference, a Δχ² = 8.0 is not big enough for strong statistical support. This is particularly the case considering that small photometric variations can arise from the use of multiple data sets processed using different pipelines. Considering that microlensing search algorithms are set up to systematically exclude light curves with periodically varying baselines, the evidence of no detectable periodic baseline variation is not strong either to strongly support the 3L1S model. In order to strongly support the argument based on the low prior probability of observing a nearly equal-luminosity pairing of giant branch stars, one needs to know the relative probability of nearly equal-luminosity binary giant stars versus planets in binary systems detectable to microlensing. Without the information on the frequency of cool planets in binaries, this argument does not strongly support the 3L1S interpretation of the event. We, therefore, consider the 2L2S model as a viable solution.

The physical lens parameters can be constrained by measuring observables related to the parameters. The first of such observables is the event timescale, which is related to the lens parameters of the mass M and distance D_L by

$$\begin{eqnarray}\begin{array}{rcl}{t}_{{\rm{E}}} & = & \displaystyle \frac{{\theta }_{{\rm{E}}}}{\mu };\ {\theta }_{{\rm{E}}}={\left(\kappa M{\pi }_{\mathrm{rel}}\right)}^{1/2};\\ {\pi }_{\mathrm{rel}} & = & \mathrm{au}\left(\displaystyle \frac{1}{{D}_{{\rm{L}}}}-\displaystyle \frac{1}{{D}_{{\rm{S}}}}\right).\end{array}\end{eqnarray} \tag{ 1 }$$

The timescale is related to the three physical parameters M, μ and D_L, and thus the constraint on the lens parameters is weak. The lens parameters can be more tightly constrained with the measurement of the angular Einstein radius because θ_E is related to two parameters of M and D_L. In this section, we determine the angular Einstein radius for the use of constraining the physical lens parameters. In Section 5, we discuss the procedure of determining M and D_L in detail. We note that the lens parameters can be uniquely determined by additionally measuring the microlens parallax by

$$\begin{eqnarray}&&M=\displaystyle \frac{{\theta }_{{\rm{E}}}}{\kappa {\pi }_{{\rm{E}}}};\qquad {D}_{{\rm{L}}}=\displaystyle \frac{\mathrm{au}}{{\pi }_{{\rm{E}}}{\theta }_{{\rm{E}}}+{\pi }_{{\rm{S}}}}.\end{eqnarray} \tag{ 2 }$$

However, π_E cannot be measured for OGLE-2019-BLG-0304.

The angular Einstein radius is measured from the combination of the angular source radius θ_* and the normalized source radius ρ by

$$\begin{eqnarray}&&{\theta }_{{\rm{E}}}=\displaystyle \frac{{\theta }_{* }}{\rho }.\end{eqnarray} \tag{ 3 }$$

The value of ρ is measured by analyzing the part of the lensing light curve affected by finite-source effects, and it is presented in Table 1. The value of θ_* is estimated from the color and brightness of the source star. For the θ_* measurement, we use the standard method of Yoo et al. (2004). According to this method, we first calibrate the color and magnitude of the source using those of the red giant clump (RGC) centroid as a reference, and then estimate θ_* using the relation between the color and θ_*.

Figure 10 shows the locations of the source and RGC centroid in the CMD of stars lying in the vicinity of the source. The source position estimated based on the 3L1S solution is marked by a red dot and the positions of the individual binary source stars estimated based on the 2L2S solution are marked by black and gray dots. The CMD is constructed based on the pyDIA photometry of the KMTC data set. The measured instrumental colors and magnitudes of the source according to the 3L1S solution and RGC centroid are (V − I, I) = (3.681 ± 0.009, 17.950 ± 0.001) and (V − I, I)_RGC = (3.785, 17.964), respectively. From the offsets in color Δ(V − I) and magnitude ΔI between the source and RGC centroid together with the known extinction and reddening corrected values of the RGC centroid, (V − I, I)_RGC,0 = (1.060, 14.564), toward the field from Bensby et al. (2013) and Nataf et al. (2013), the dereddened color and magnitude of the source are determined as

$$\begin{eqnarray}\begin{array}{rcl}{(V-I,I)}_{0} & = & {(V-I,I)}_{\mathrm{RGC},0}+{\rm{\Delta }}(V-I,I)\\ & = & \,(0.956\pm 0.009,14.550\pm 0.001).\end{array}\end{eqnarray} \tag{ 4 }$$

According to the 3L1S solution, the source is an RGC star with an early-K spectral type. As mentioned in the previous section, the stellar types of S₁ and S₂ according to the 2L2S solution are G- and K-type giants, respectively.

The angular source radius is deduced from the measured source color and magnitude. For this, we first convert the measured V − I color into a V − K color using the color–color relation of Bessell & Brett (1988). We then estimate θ_* from V − K using the θ_*–(V − K) relation of Kervella et al. (2004). The estimated source radius from this procedure is

$$\begin{eqnarray}&&{\theta }_{* }=5.26\pm 0.37\,\mu \mathrm{as},\end{eqnarray} \tag{ 5 }$$

for the 3L1S solution. The angular Einstein radius estimated using the relation in Equation (3) is

$$\begin{eqnarray}&&{\theta }_{{\rm{E}}}=0.285\pm 0.020\,\mathrm{mas}.\end{eqnarray} \tag{ 6 }$$

This together with the event timescale yields the relative lens-source proper motion of

$$\begin{eqnarray}&&\mu =\displaystyle \frac{{\theta }_{{\rm{E}}}}{{t}_{{\rm{E}}}}=5.86\pm 0.41\,\mathrm{mas}\,{\mathrm{yr}}^{-1}.\end{eqnarray} \tag{ 7 }$$

The Einstein radius and proper motion based on the 2L2S solution result in different values from those estimated based on the 3L1S solution. This is because the flux from the source is divided into two roughly equal-luminosity stars and thus the source radii of the individual source stars are smaller than that of a single source for the 3L1S solution. The estimated angular radius of S₁, the Einstein radius, and the proper motion are ${\theta }_{* ,{S}_{1}}=3.24\pm 0.28$, ${\theta }_{{\rm{E}}}={\theta }_{* ,{S}_{1}}/{\rho }_{1}=0.177\pm 0.016$ mas, and μ = 3.69 ± 0.32 mas yr⁻¹, respectively. We note that the values of θ_E and μ are smaller than those estimated from the 3L1S solution due to the smaller source radius.

Although it is tough to uniquely determine M and D_L due to the difficulty of measuring π_E, the lens parameters can still be constrained with the measured observables of t_E and θ_E. For this, we conduct a Bayesian analysis of the event based on the priors of the lens mass function and Galactic model of the physical and dynamic distributions.

In the Bayesian analysis, we conduct a Monte Carlo simulation to produce a large number of artificial lensing events. The mass function used in the simulation is adopted from Zhang et al. (2020) for stellar and brown-dwarf lenses and Gould (2000) for remnant lenses, including white dwarfs, neutron stars, and black holes. The lens objects are physically distributed based on the modified Han & Gould (2003) model, in which the distribution of disk matter is modified from the original version using the model of Bennett et al. (2014). The motion of the lens is assigned using the dynamical model of Han & Gould (1995). The proper motion of the source, μ _S(R. A.,decl.) = (−6.737 ± 1.721, −9.277 ± 1.284) mas yr⁻¹, is known from the measurement by Gaia (Gaia Collaboration et al. 2018), and thus we consider the measured source proper motion in the computation of the relative lens-source motion. With these priors, we produce 4 × 10⁷ artificial events and then construct probability distributions of M and D_L for events with t_E and θ_E values lying within the ranges of the measured observables. Then, the median values of the probability distributions are presented as representative values of the physical lens parameters and the uncertainties are estimated as the 1σ range of the probability distributions, that is, 16% for the lower limit and 84% for the upper limit. Considering that the 3L1S and 2L2S solutions result in similar fits to the observed data, we conduct two sets of analysis based on the two solutions.

The shaded curves in Figure 11 show the posterior distributions of the M₁ (upper panel) and D_L (lower panel) based on the 3L1S solution. The estimated masses of the individual lens components are

$$\begin{eqnarray}\begin{array}{rcl}{M}_{1} & = & {0.27}_{-0.12}^{+0.27}\,{M}_{\odot },\\ {M}_{2} & = & {0.10}_{-0.04}^{+0.10}\,{M}_{\odot },\\ {M}_{3} & = & {0.51}_{-0.23}^{+0.51}\,{M}_{{\rm{J}}},\end{array}\end{eqnarray} \tag{ 8 }$$

indicating that the lens is a planetary system in which a giant plant belongs to a binary composed of a mid-M dwarf and a late-M dwarf. The estimated distance to the lens is

$$\begin{eqnarray}&&{D}_{{\rm{L}}}={6.98}_{-0.90}^{+0.91}\,\mathrm{kpc}.\end{eqnarray} \tag{ 9 }$$

The projected physical separations of the stellar (M₂) and planetary (M₃) companions from the primary (M₁) are

$$\begin{eqnarray}\begin{array}{rcl}{a}_{\perp ,2} & = & {s}_{2}{D}_{{\rm{L}}}{\theta }_{{\rm{E}}}={5.19}_{-4.52}^{+5.87}\,\mathrm{au},\\ {a}_{\perp ,3} & = & {s}_{3}{D}_{{\rm{L}}}{\theta }_{{\rm{E}}}={1.23}_{-1.07}^{+1.39}\,\mathrm{au},\end{array}\end{eqnarray} \tag{ 10 }$$

respectively. Considering that a_⊥,3 is substantially smaller than a_⊥,2, the planet is very likely to be in an S-type orbit around the heavier star of the binary. In Table 2, we summarize the estimated physical parameters of M₁, M₂, M₃, D_L, a_⊥,2, and a_⊥,3. We find that the Gaia measurement of the vector source proper motion, μ _S(l, b) = ( − 11.4, + 0.7) mas yr⁻¹, when combined with the scalar lens-source relative proper motion from Equation (7), strongly constrains the lens to lie in the bulge. Here μ _S(l, b) denotes the source proper motion vector in the Galactic coordinates, and it is related to μ _S(R. A.,decl.) by

$$\begin{eqnarray}&&\left(\begin{array}{c}{\mu }_{{\rm{S}},l}\\ {\mu }_{{\rm{S}},b}\end{array}\right)=\left(\begin{array}{cc}\cos \eta & \sin \eta \\ -\cos \eta & \cos \eta \end{array}\right)\left(\begin{array}{c}{\mu }_{{\rm{S}},{\rm{R}}.{\rm{A}}.}\\ {\mu }_{{\rm{S}},\mathrm{decl}.}\end{array}\right),\end{eqnarray} \tag{ 11 }$$

where η ∼ 576 is the tilt angle of the Galactic plane with respect to the celestial equator. This is because the Galactic model contains relatively few disk lenses with μ_L(l) ≲ − 5.8 mas yr⁻¹. Specifically, we find that when the Gaia measurement is ignored, the probability of a disk lens is 28%. However, after including the Gaia measurement, the probability of a disk lens is only ∼1%.

Zoom In Zoom Out Reset image size

Table 2. Physical Lens Parameters

Parameter	3L1S	2L2S
M1 (M⊙)	${0.27}_{-0.12}^{+0.27}$	${0.12}_{-0.05}^{+0.12}$
M2 (M⊙)	${0.10}_{-0.04}^{+0.10}$	${0.045}_{-.019}^{+0.045}$
M3 (MJ)	${0.51}_{-0.23}^{+0.51}$	⋯
DL (kpc)	${6.98}_{-0.90}^{+0.91}$	${6.97}_{-0.86}^{+0.70}$
a⊥,2 (au)	${5.19}_{-4.52}^{+5.87}$	${3.52}_{-3.09}^{+3.87}$
a⊥,3 (au)	${1.23}_{-1.07}^{+1.39}$	⋯

Download table as:  ASCIITypeset image

The Bayesian posterior distributions obtained with the 2L2S solution are also presented in Figure 11 (red unshaded curves). The estimated masses of the lens components and the distance to the lens system are

$$\begin{eqnarray}\begin{array}{rcl}{M}_{1} & = & {0.12}_{-0.05}^{+0.12}\,{M}_{\odot },\\ {M}_{2} & = & {0.045}_{-.019}^{+0.045}\,{M}_{\odot },\end{array}\end{eqnarray} \tag{ 12 }$$

and

$$\begin{eqnarray}&&{D}_{{\rm{L}}}={6.97}_{-0.86}^{+0.70}\,\mathrm{kpc},\end{eqnarray} \tag{ 13 }$$

respectively. Then, the lens is a binary composed of an M dwarf and a brown dwarf located in the bulge. The lens mass from the 2L2S solution, M_2L2S = M₁ + M₂, is smaller than the value estimated from the 3L1S solution, M_3L1S, because the Bayesian input value of the Einstein radius, θ_E,2L2S ∼ 0.18 mas, is smaller than the value of the 3L1S solution, θ_E,3L1S ∼ 0.29 mas. Considering that the lens distances for the two solutions are similar to each other, it is found that the lens masses estimated from the two solutions are approximately in the relation ${M}_{2{\rm{L}}2{\rm{S}}}\sim {({\theta }_{{\rm{E}},2{\rm{L}}2{\rm{S}}}/{\theta }_{{\rm{E}},3{\rm{L}}1{\rm{S}}})}^{2}{M}_{3{\rm{L}}1{\rm{S}}}$. The physical lens parameters of the 2L2S solution are also summarized in Table 2.

If the 3L1S model is correct, OGLE-2019-BLG-0304LAbB has the second lowest formal significance, Δχ² = 48.3, of any reported microlensing planet, with the lowest being Δχ² = 47 for OGLE-2018-BLG-0677 (Herrera-Martín et al. 2020) and the previous second lowest being Δχ² = 170 for KMT-2018-BLG-1025 (Han et al. 2021). Moreover, OGLE-2019-BLG-0304LAbB is a more complex (3L1S) system, whereas the other low formal-significance systems were 2L1S. This distinction is important because, of the seven previously claimed binary+planet microlensing systems, two were subsequently shown to be spurious, namely MACHO-97-BLG-41 (Bennett et al. 1999; Albrow et al. 2000; Jung et al. 2013) and OGLE-2013-BLG-0723 (Udalski et al. 2015; Han et al. 2016). In both cases, the original models did not incorporate orbital motion. The additional features in the light curve, which were inconsistent with static 2L1S models (and so were attributed to a third body), were found to be explained by the motion of a caustic (due to the underlying binary motion) within the context of orbital-motion models. This was the major motivation for our check in Section 3.1 of 2L1S orbital-motion models, which could not explain the additional "planetary" anomalies. However, it is also important to understand at a deeper level why orbital motion is strongly constrained to the point that it cannot explain these anomalies.

Table 3 shows the results of the fits that exclude data from the second bump. For the 2L1S model, the parameters (s₂, α₂) = (3.900 ± 0.040, 16988 ± 044) essentially predict the location of the second bump, without any direct light curve information about the bump, which peaks 61 days later. In the full-data (static) model, this bump location fixes these parameters with much greater precision: (s₂, α₂) = (3.784 ± 0.015, 16925 ± 011). We will comment on the 3σ discrepancy in s₂ below, but for the moment, the main point is that the overall agreement between this bump-free "prediction" and the actual location of the bump 61 days later, implies that the orbital motion on the 1 day timescales of the main peak (where the planetary anomaly appears) is extremely well constrained.

Table 3. 2L1S and 3L1S Models without the Second Bump

Parameter	2L1S	3L1S
χ2	4728.8	4699.5
t0 ($\mathrm{HJD}^{\prime}$)	8543.617 ± 0.005	8543.606 ± 0.006
u0 (10−3)	1.84 ± 0.46	0.98 ± 0.55
tE (days)	17.43 ± 0.15	17.62 ± 0.19
s2	3.900 ± 0.040	3.837 ± 0.063
q2	0.436 ± 0.015	0.396 ± 0.022
α (rad)	2.965 ± 0.008	2.965 ± 0.009
s3	⋯	0.879 ± 0.008
q3 (10−3)	⋯	1.67 ± 0.26
ψ (rad)	⋯	2.514 ± 0.109
ρ (10−2)	1.97 ± 0.04	1.88 ± 0.05

Note. HJD' ≡ HJD–2,450,000.

Download table as:  ASCIITypeset image

As shown by the Figure 2 residuals, the only really pronounced deviations from the 2L1S model come from the three KMTC points on each of HJD' = 8542.xx and 8543.xx. These contribute Δχ² = 7 and Δχ² = 18 to the total Δχ² = 48. Moreover, they could contribute indirectly to the Δχ² = 10 from the early light curve (because the KMTS data, with similar coverage during this period, do not show any significant Δχ²). We therefore focus particular attention on these two nights.

We note that the crescent moon passed within 61 of the event at HJD' = 8542.37, and during the three observations on 8542.xx, the moon was separated by 84–88, which gave rise to background counts of about 1900 per pixel (compared to ∼500 for dark time and ∼25,000 when the full moon is in the bulge). However, on the next night, which yielded a much larger "planetary signal," the Moon was 191–196 from the event, and the background was only about 800. Hence, it is unlikely that these relatively normal (especially on the second night) observing conditions could be responsible for the planetary anomaly. To double check, we investigate the effect of the moon by additionally (1) conducting visual inspection of images, (2) probing the dependence of the photometry on the moon phase and distance, and (3) comparing our photometric results with those of other comparison stars. From this, we find nothing irregular about them.

Finally, the Δχ² = 19 improvement when the bump is included (see Figure 9) lends added weight to the reality of the signal. This improvement arises from the greater consistency between the "predicted" and actual forms of the late-time bump when the planet is included in the model. That is, regardless of the origin of the deviations in the KMTC light curve (i.e., whether a planet or some random systematics), the 2L1S and 3L1S fits will try to adjust their parameters to accommodate these deviations. If the "planet" is just the result of accommodating random systematics in the 3L1S fit, then the 3L1S prediction should, on average, be no better or worse than the 2L1S prediction. Using the parameters presented in Tables 1 and 3, we find that the differences in the parameters between the two sets of models obtained with the full data and the partial data without the second bump are Δ(s₂, q₂, α) = (0.116, 0.044, 063) ± (0.040, 0.015, 044) for the 2L1S model and Δ(s₂, q₂, α) = (0.104, 0.032, 029) ± (0.063, 0.022, 050) for the 3L1S model. In all cases, there is improved agreement. And, by adding the planet, the χ² drops from 19.1 to 5.3 for three degrees of freedom. We therefore conclude that the residual from the 2L1S model cannot be ascribed to the lens orbital motion, although its 3L1S or 2L2S origin cannot be firmly distinguished.

We analyzed the microlensing event OGLE-2019-BLG-0304. The light curve of the event showed two distinctive features, in which the main peak appeared to exhibit deviations caused by finite-source effects, and a second bump appearing ∼61 days after the main peak. The fit with a finite-source single-lens model excluding the second bump left substantial residuals in the peak region, indicating that a more sophisticated model was needed to explain not only the second bump but also the deviations in the peak region. A 2L1S model could explain the overall features of the light curve by significantly reducing the residuals in both regions of the main peak and the second bump, but the model still left subtle but noticeable residuals in the peak region. We found that the residuals could be explained by the presence of either a planetary companion located close to the primary of the lens or an additional close companion to the source. Although the 3L1S model is favored over the 2L2S model, firm resolution of the degeneracy between the two models was difficult with the photometric data. Therefore, the event well illustrated the need for thorough model testing in the interpretation of lensing events with complex features.

We estimated the physical lens parameters expected from the two degenerate solutions by conducting Bayesian analysis. According to the 3L1S interpretation, the lens is a planetary system in which a planet with a mass ${0.51}_{-0.23}^{+0.51}\,{M}_{{\rm{J}}}$ is in an S-type orbit around a binary composed of stars with masses ${0.27}_{-0.12}^{+0.27}\,{M}_{\odot }$ and ${0.10}_{-0.04}^{+0.10}\,{M}_{\odot }$. According to the 2L2S interpretation, the source is composed of G- and K-type giant stars, and the lens is composed of a low-mass M dwarf and a brown dwarf with masses ${0.12}_{-0.05}^{+0.12}\,{M}_{\odot }$ and ${0.045}_{-.019}^{+0.045}\,{M}_{\odot }$, respectively.

Work by C.H. was supported by grants from the National Research Foundation of Korea (2019 R1A2C2085965 and 2020 R1A4A2002885). Work by A.G. was supported by JPL grant 1500811. This research has made use of the KMTNet system operated by the Korea Astronomy and Space Science Institute (KASI) and the data were obtained at three host sites of CTIO in Chile, SAAO in South Africa, and SSO in Australia. The OGLE project has received funding from the National Science Centre, Poland, grant MAESTRO 2014/14/A/ST9/00121 to AU.